loadct,39

;l=mrdfits('UKIDSS_GCS_SDSS_z_Z.fits',1)
l=mrdfits('UKIDSS_GCS_SDSS_z_Z_DR9.fits',1)

;mags='MODEL'
mags='PETRO'

l=l[where(l.specClass eq 2)]
n_obj=n_elements(l)

tags_ph=tag_names(l)
tag_u = where(tags_ph eq mags+'MAG_U')
tag_g = where(tags_ph eq mags+'MAG_G')
tag_r = where(tags_ph eq mags+'MAG_R')
tag_i = where(tags_ph eq mags+'MAG_I')
tag_z = where(tags_ph eq mags+'MAG_Z')   
tag_uerr = where(tags_ph eq mags+'MAGERR_U')
tag_gerr = where(tags_ph eq mags+'MAGERR_G')
tag_rerr = where(tags_ph eq mags+'MAGERR_R')
tag_ierr = where(tags_ph eq mags+'MAGERR_I')
tag_zerr = where(tags_ph eq mags+'MAGERR_Z')

mags=dblarr(9,n_obj)
mags[0,*]=l.(tag_u)-0*l.extinction_u - 0.927
mags[1,*]=l.(tag_g)-0*l.extinction_g + 0.103
mags[2,*]=l.(tag_r)-0*l.extinction_r - 0.146
mags[3,*]=l.(tag_i)-0*l.extinction_i - 0.366
mags[4,*]=l.(tag_z)-0*l.extinction_z - 0.533
;;;mags[5,*]=l.ZPETROMAG + 0.634*0
mags[5,*]=l.YPETROMAG + 0.634*0
mags[6,*]=l.JPETROMAG + 0.938*0
mags[7,*]=l.HPETROMAG + 1.379*0
mags[8,*]=l.K_1PETROMAG + 1.900*0

;;;diffzZ=mags[4,*]-(l.ZPETROMAG+0.015-0.103*(l.ZPETROMAG-l.YPETROMAG-0.106))
diffzZ=mags[4,*]-(l.ZPETROMAG+0.015-0.103*(l.ZPETROMAG-l.YPETROMAG))
diffzY=dblarr(n_obj)

for i=0,n_obj-1 do begin
    mags_new=nir_col_corr(mags[*,i], diffzY=d,ndeg=3)
    diffzY[i]=d
endfor

o=where(finite(diffzy-diffzz) eq 1 and abs(diffzy) le 2 and abs(diffzz) le 2)
g=where(((diffzy[o]+0.584) gt -0.0) and ((diffzy[o]+0.584) lt 0.65) $
    and (diffzz[o] gt -0.4) and (diffzz[o] lt 0.8) $
    and l[o].petromagerr_z le 0.05 and l[o].zpetromagerr le 0.05 $
    and l[o].petror50_r gt 0.4 and l[o].petror50_r lt 30.0 $
    ,cg)


tt=poly_fit(diffzy[o[g]]+0.584,diffzz[o[g]],2)
var2=alog10(l.petror50_r)
;;var2=l.petromag_r+5*alog10(l.petror50_r)-17.0

d_arr=dblarr(3,cg)
d_arr[0,*]=transpose(diffzy[o[g]]+0.584)
d_arr[1,*]=transpose(var2[o[g]])
d_arr[2,*]=transpose(diffzz[o[g]])
cf=sfit_2deg(d_arr,2,2,/irr,kx=kx)

diffzz_mod_1d=poly(diffzy+0.584,tt)
diffzz_mod=poly2d(diffzy+0.584,var2,kx,/irr,deg1=2,deg2=2)

plot,diffzy+0.584,diffzz,psym=4,xr=[-0.4,1.2],yr=[-0.8,0.8],/iso
oplot,diffzy[o[g]]+0.584,diffzz[o[g]],psym=4,col=128
;;oplot,diffzy[g]+0.584,poly(diffzy[g]+0.584,tt),col=254,psym=7
oplot,diffzy[o]+0.584,diffzz[o]-diffzz_mod[o]-0.5,col=254,psym=7,syms=0.4

print,robust_sigma(diffzz[o[g]]-diffzz_mod[o[g]]),$
      stdev(diffzz[o[g]]-diffzz_mod[o[g]])
print,robust_sigma(diffzz[o[g]]-diffzz_mod_1d[o[g]]),$
      stdev(diffzz[o[g]]-diffzz_mod_1d[o[g]])


end