PRO plot_phangs_ssps_against_uv_fast

;plot_phangs_ssps_against_uv_fast

dustem_define_la_common

source_name='ngc0628'
data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'
obp=[1.1,0,1.15,1]

;restore,data_dir+'ngc0628_astrosat_minus_prediction_fast.sav',/verb
restore,data_dir+'ngc0628_astrosat_voronoi_prediction_fast.sav',/verb
;% RESTORE: Restored variable: ACTUAL_SEDS.
;% RESTORE: Restored variable: PREDICTED_SEDS.
;% RESTORE: Restored variable: HREF.
;% RESTORE: Restored variable: ALL_SEDS_INDICES.
;% RESTORE: Restored variable: DO_NORMALIZE.
;% RESTORE: Restored variable: ALL_FILTERS.
;% RESTORE: Restored variable: NHVOR.
restore,data_dir+'ngc0628_astrosat_data.sav',/verb
;% RESTORE: Restored variable: ASTROSAT_IMAGES.
;% RESTORE: Restored variable: ASTROSAT_FILTERS.
;% RESTORE: Restored variable: HREF.
st_muse_weights=read_muse_phangs_weights(object='NGC0628',bin_numbers=bin_numbers)

Nvor=n_elements(st_muse_weights)
ebv_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))

FOR i=0L,Nvor-2 DO BEGIN    ;why -2 ??
	ebv_map[*all_seds_indices[i]]=st_muse_weights[i].reddening
ENDFOR

win=0L
window,win,xsize=800,ysize=800 & win=win+1

imrange=[-0.02,0.4]
image_cont20,ebv_map,href,/silent,/square,image_color_table='jpbloadct',axis_color_table=1,imrange=imrange,off_bar_pos=obp

stop

window,win,xsize=800,ysize=800 & win=win+1
k=0L
tit='Astrosat '+ astrosat_filters[k]
imrange=[-0.01,0.2]
image_cont20,astrosat_images[*,*,0,0],href,/silent,/square,image_color_table='jpbloadct',axis_color_table=1,imrange=imrange,off_bar_pos=obp

;stop

window,win,xsize=800,ysize=800 & win=win+1
!p.multi=[0,1,2]

;Nbins=11
Nbins=101
yrange=[1.,1.e5]
;==== plot the histogram of data-model in MJy/sr
xtit='(data-model) [MJy/sr]'
;vmin=-0.1 & vmax=0.1 & Nbins=50
;vmin=-1. & vmax=1. & Nbins=50
;vmin=-5. & vmax=5. & Nbins=100
vmin=-10. & vmax=10.
ims=la_sub(actual_seds,predicted_seds)
im=ims[*,0]
ind=where(im NE la_undef())
res=histogram(im,locations=x,min=vmin,max=vmax,Nbins=Nbins)
x=x+(vmax-vmin)/Nbins/2.
;yrange=[max([min(res),1.]),1.e4]
colors=['violet','blue','yellow','brown','red']
cgplot,x,res,/ylog,xtit=xtit,ytit='N',psym=10,color='violet',tit=tit,/nodata,yrange=yrange,/ysty
FOR k=0L,4 DO BEGIN
	;im=ims[*,k]
	;ind=where(im NE la_undef())
	;res=histogram(im[ind],locations=x,min=vmin,max=vmax,Nbins=Nbins)
	res=histogram(ims[*,k],locations=x,min=vmin,max=vmax,Nbins=Nbins)
	x=x+(vmax-vmin)/Nbins/2.
   	cgoplot,x,res,color=colors[k],psym=10
ENDFOR

;stop

;==== plot the histogram of data-model in %
;imsp=la_mul(la_div(la_sub(actual_seds,predicted_seds),predicted_seds),100.)
imsp=la_mul(la_div(ims,predicted_seds),100.)
;imsp=ims/predicted_seds*100.
vmin=-500. & vmax=500.
xtit='(data-model)/model [%]'
im=imsp[*,0]
ind=where(im NE la_undef())
res=histogram(im,locations=x,min=vmin,max=vmax,Nbins=Nbins)
x=x+(vmax-vmin)/Nbins/2.
;yrange=[max([min(res),1.]),1.e4]
colors=['violet','blue','yellow','brown','red']
cgplot,x,res,/ylog,xtit=xtit,ytit='N',psym=10,color='violet',tit=tit,/nodata,yrange=yrange,/ysty
FOR k=0L,4 DO BEGIN
   ;im=imsp[*,k]
   ;ind=where(im NE la_undef())
   ;res=histogram(im[ind],locations=x,min=vmin,max=vmax,Nbins=Nbins)
   res=histogram(imsp[*,k],locations=x,min=vmin,max=vmax,Nbins=Nbins)
   x=x+(vmax-vmin)/Nbins/2.
   cgoplot,x,res,color=colors[k],psym=10
ENDFOR

stop

;=== why is there so many positive differences in % ?

;==== plot the histogram of data-model above some NH
;NHthreshold=0.5 & tit='NH > 5e20 H/cm2'
;NHthreshold=1. & tit='NH > 1e21 H/cm2'
NHthreshold=5. & tit='NH > 5e21 H/cm2'
;NHthreshold=10. & tit='NH > 1e22 H/cm2'
;NHthreshold=50. & tit='NH > 5e22 H/cm2'
Nbins=101
yrange=[1.,1.e5]
;==== plot the histogram of data-model in MJy/sr
xtit='(data-model) [MJy/sr]'
;vmin=-0.1 & vmax=0.1 & Nbins=50
;vmin=-1. & vmax=1. & Nbins=50
;vmin=-5. & vmax=5. & Nbins=100
vmin=-10. & vmax=10.
ims=la_sub(actual_seds,predicted_seds)
im=ims[*,0]
ind=where(im NE la_undef())
res=histogram(im,locations=x,min=vmin,max=vmax,Nbins=Nbins)
x=x+(vmax-vmin)/Nbins
;yrange=[max([min(res),1.]),1.e4]
colors=['violet','blue','yellow','brown','red']
cgplot,x,res,/ylog,xtit=xtit,ytit='N',psym=10,color='violet',tit=tit,/nodata,yrange=yrange,/ysty
FOR k=0L,4 DO BEGIN
	;im=ims[*,k]
	;ind=where(im NE la_undef())
	;res=histogram(im[ind],locations=x,min=vmin,max=vmax,Nbins=Nbins)
	ind=where(NHvor GT NHthreshold)
	res=histogram(ims[ind,k],locations=x,min=vmin,max=vmax,Nbins=Nbins)
	x=x+(vmax-vmin)/Nbins
   	cgoplot,x,res,color=colors[k],psym=10
ENDFOR

;stop

;==== plot the histogram of data-model in %
;imsp=la_mul(la_div(la_sub(actual_seds,predicted_seds),predicted_seds),100.)
imsp=la_mul(la_div(ims,predicted_seds),100.)
;imsp=ims/predicted_seds*100.
vmin=-500. & vmax=500.
xtit='(data-model)/model [%]'
im=imsp[*,0]
ind=where(im NE la_undef())
res=histogram(im,locations=x,min=vmin,max=vmax,Nbins=Nbins)
x=x+(vmax-vmin)/Nbins
;yrange=[max([min(res),1.]),1.e4]
colors=['violet','blue','yellow','brown','red']
cgplot,x,res,/ylog,xtit=xtit,ytit='N',psym=10,color='violet',tit=tit,/nodata,yrange=yrange,/ysty
FOR k=0L,4 DO BEGIN
   ;im=imsp[*,k]
   ;ind=where(im NE la_undef())
   ;res=histogram(im[ind],locations=x,min=vmin,max=vmax,Nbins=Nbins)
   ind=where(NHvor GT NHthreshold)
   res=histogram(imsp[ind,k],locations=x,min=vmin,max=vmax,Nbins=Nbins)
   x=x+(vmax-vmin)/Nbins
   cgoplot,x,res,color=colors[k],psym=10
ENDFOR

stop

;==== make maps
window,0,xsize=800,ysize=800
window,1,xsize=800,ysize=800
window,2,xsize=800,ysize=800
window,3,xsize=800,ysize=800

Nvor=(size(actual_seds))[1]
diff_p=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),5)
diff=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),5)
diff_n=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),5)
NHmap=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
ebv_map_fit=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
NHmap[*,*]=la_undef()
diff_p[*,*,*]=la_undef()
diff[*,*,*]=la_undef()
diff_n[*,*,*]=la_undef()
imrange=[-0.1,0.1]
obp=[1.1,0,1.15,1]
imrangep=[-200,200.]

FOR i=0L,Nvor-1 DO BEGIN
	NHmap[*all_seds_indices[i]]=NHvor[i]
	ebv_map_fit[*all_seds_indices[i]]=ebv[i] ;This is only a correction, as the E(B-V) Muse was already applied when doinf the SEDs
ENDFOR
ebv_tot=la_add(ebv_map,ebv_map_fit)
ind=where(ebv_tot LE 0. and EBV_tot NE la_undef(),count)
ebv_tot[ind]=0.
ebv_tot_vor=la_add(ebv,st_muse_weights[0:47068-1].reddening)  ;This is E(B-V) total on voronoi bins
ind=where(ebv_tot_vor LT 0 AND ebv_tot_vor NE la_undef(),count)
print,(1.*count)/Nvor*100.,' %'
;      23.8570 %  when using average of 1 astrosat bands to set E(B-V) and normalizing at SDSS4
;      18.5774 %  when using average of 5 astrosat bands to set E(B-V) and normalizing at SDSS4
;      14.6087 %  when using average of 5 astrosat bands to set E(B-V) and normalizing at NIRCAM11
print,la_median(ebv_tot_vor[ind])
;    -0.099808343 when using average of 5 astrosat bands to set E(B-V) and normalizing at SDSS4
;    -0.088146150 when using average of 5 astrosat bands to set E(B-V) and normalizing at NIRCAM11
stop

wset,0
tit='E(B-V) correction SED fit []'
;imrange_ebv=[-0.1,1.]
imrange_ebv=[-0.02,0.4]
image_cont20,ebv_map_fit,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrange_ebv,title=tit

stop

tit='E(B-V) total (Muse + SED fit) []'
image_cont20,ebv_tot,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrange_ebv,title=tit

stop

cgplot,NHvor,ebv+st_muse_weights.reddening,psym=3,xtit='NH',ytit='E(B-V) tot'

stop

FOR k=0L,4 DO BEGIN
	mapp=diff_p[*,*,k]
	map=diff[*,*,k]
	FOR i=0L,Nvor-1 DO BEGIN
		IF i mod 1000 EQ 0 THEN BEGIN
			message,'Done '+strtrim(1.*i/Nvor*100,2)+' %',/cont
		ENDIF
		map[*all_seds_indices[i]]=la_sub(actual_seds[i,k],predicted_seds[i,k])
		mapp[*all_seds_indices[i]]=la_mul(la_div(map[*all_seds_indices[i]],predicted_seds[i,k]),100.)
	ENDFOR
	diff_p[*,*,k]=mapp
	diff[*,*,k]=map
	diff_n[*,*,k]=la_div(map,NHmap)

	;image_cont20,map,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrange
    wset,1
	tit='Data-model [%] '+all_filters[k]
	image_cont20,diff_p[*,*,k],href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrangep,title=tit
    wset,2
	tit='Data-model [MJy/sr] '+all_filters[k]
    image_cont20,diff[*,*,k],href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrange,title=tit
    wset,3
	tit='Data-model Normalised [MJy/sr] '+all_filters[k]
    image_cont20,diff_n[*,*,k],href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrange,title=tit
    filename='/tmp/'+source_name+'_diff_map_perc_astrosat'+strtrim(k+1,2)+'.fits'
    writefits,filename,diff_p[*,*,k],href
    message,'Wrote '+filename,/info
	stop
ENDFOR

stop

END