PRO check_phangs_ssps_against_uv,show_map=show_map

;check_phangs_ssps_against_uv

;===== This computes a map of prediction by dustem_plugin_phangs_stellar_continuum
;===== constrained on JWST NIR data to astrosat data in the UV for NGC0628.
;===== uses the plugin in a mode where the fortran is not run, but yet, the output files of the fortran are read at each pixel,
;===== which is not efficient (supposed to run in 46hrs .... See  check_phangs_ssps_against_uv_fast for an alternative use.

window,0,xsize=900,ysize=1000
pdp_define_la_common
obp=[1.1,0,1.15,1]
source_name='ngc0628'
data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'

;dustem_init,show_plots=show_plots
;needed only for NHCO
restore,data_dir+'ngc0628_jwst_images.sav',/verb
;needed for stellar parameters and voronoi bin info
restore,data_dir+'ngc0628_muse_images.sav',/verb
;restore,data_dir+'ngc0628_muse_NH.sav',/verb
;restore,data_dir+'ngc0628_astrosat_data.sav',/verb

;stop

;=== restore SEDs in voronoi bins
restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb
jwst_seds=all_seds
jwst_filters=(*jwst_seds[0]).filter
restore,data_dir+'ngc0628_muse_seds_muse_pixels.sav',/verb
muse_seds=all_seds
muse_filters=(*muse_seds[0]).filter
restore,data_dir+'ngc0628_astrosat_seds_muse_pixels.sav',/verb
astrosat_seds=all_seds
astrosat_filters=(*astrosat_seds[0]).filter
all_filters=[astrosat_filters,muse_filters,jwst_filters]

Nvor=max(voronoi_id)
Nfilters=n_elements(jwst_filters)+n_elements(astrosat_filters)+n_elements(muse_filters)

use_NHmap=NHCO       ;used NH map in 1.e21 (from CO)

;galaxy_distance=10.  ;MPc

use_model='DBP90'    ;Example with default keywords uses the DBP90 model
;use_model='MC10'    ;Example with default keywords uses the DBP90 model
;use_model='DL07'    ;Example with default keywords uses the DBP90 model
use_polarization=0   ; initialize Dustemwrap in no polarization mode 
use_window=2          ; default graphics window number to use for plotting the results
use_verbose=0
use_Nitermax=5        ; maximum number of iterations for the fit
dustem_define_la_common
	;parameters to fit
pd = [ $
     '(*!dustem_params).G0', $                           ;G0
     '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
     '(*!dustem_params).grains(1).mdust_o_mh',$          ;PAH1 mass fraction
     ;'(*!dustem_params).grains(2).mdust_o_mh', $         ;VSGs
     ;'(*!dustem_params).grains(3).mdust_o_mh', $         ;VSGs
     ;'(*!dustem_params).grains(4).mdust_o_mh', $         ;BGs
     'dustem_plugin_phangs_stellar_continuum_1']         ;stellar continuum amplitude
;iv =   [1.6, [2.2e-3, 2.2e-3, 5.7e-3, 5.7e-3],2.e-3] 
iv =   [1.6, [2.2e-3, 2.2e-3],5.e-5] 
;iv =   [1.6, [2.2e-3, 2.2e-3,2e-3, 2e-3],5.e-5] 
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
;llims=[1.e-2,[1.e-4,1.e-4,1.e-4,1.e-4],1.e-10]
llims=[1.e-2,[1.e-4,1.e-4],1.e-10]
;llims=[1.e-2,[1.e-4,1.e-4,1.e-4,1.e-4],1.e-10]

pd_to_update=['dustem_plugin_phangs_stellar_continuum_1','(*!dustem_params).grains(0).mdust_o_mh']
pd_filter_names=['NIRCAM11','MIRI2']

;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots
;stop
!quiet=1
!dustem_nocatch=1
!dustem_verbose=use_verbose
IF keyword_set(noobj) THEN !dustem_noobj=1
;!EXCEPT=2 ; for debugging (not sure hat this is)
;=== INFORMATION TO MAKE THE PLOT
;yr=[1.00e-2,10.0] ; y-axis limits
yr=[1.00e-3,100.0] ; y-axis limits
;xr=[0.5,1.00e2] ; x-axis limits
xr=[0.1,1.00e2] ; x-axis limits
tit='FIT INTENSITY EXAMPLE' ; plot title
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title
xtit=textoidl('\lambda (\mum)') ; x-axis title

;voronoi bins for test. Last series are randomly picked
one_sed=(*jwst_seds[0])[0]
Nfilters=n_elements(astrosat_filters)+n_elements(muse_filters)+n_elements(filters)
seds=replicate(one_sed,Nfilters)
diff_seds=fltarr(Nvor,Nfilters)
diff_maps=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),Nfilters)   ;difference maps
diff_maps[*]=la_undef()
;imrange=[-1.,1.]   ;MJy/sr
imrange=[-0.3,0.2]   ;MJy/sr

t1=systime(0,/sec)
first_vid=0L
;first_vid=5960L
;first_vid=12631L
;first_vid=12640L
do_normalize=1
FOR vid=first_vid,Nvor-1 DO BEGIN
	IF vid mod 10 EQ 0 THEN BEGIN
		t2=systime(0,/sec)
		delta_t_hr=(t2-t1)/60.^2 ;elapsed time [hr]
		frac_perc=vid/Nvor*100
		delta_t_remain=delta_t_hr/frac_perc*100.-delta_t_hr  ;remaining time [hr]
		message,'done '+strtrim(frac_perc)+' % in '+strtrim(delta_t_hr,2)+' hr remaining '+strtrim(delta_t_remain,2)+' hr',/continue
		ind=where(diff_maps[*,*,0] NE la_undef(),count)
		print,la_min(diff_maps[*,*,0]),la_max(diff_maps[*,*,0]),count
		image_cont20,diff_maps[*,*,0],href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrange
		;stop
	ENDIF
    ;index=vid

	;===== get the SSP weights for the given voronoi bins
	weights=phangs_binid2weights(st_muse_weights,vid,st_templates,age_values,metalicity_values,reddening=reddening)
	;===== extract SED
	Nf=0L & Nt=Nf+n_elements(astrosat_filters)-1
	seds[Nf:Nt]=*astrosat_seds[vid]
	Nf=Nt+1 & Nt=Nf+n_elements(muse_filters)-1
	seds[Nf:Nt]=*muse_seds[vid]
	Nf=Nt+1 & Nt=Nf+n_elements(jwst_filters)-1
	seds[Nf:Nt]=*jwst_seds[vid]
	sed=seds
	ind=where(sed.stokesI EQ la_undef(),count)
	IF count NE 0 THEN goto,do_the_next
	;needed to add this because astrosat had 0 variances
	sed.sigmaii=(5./100.*sed.stokesI)^2
	;===== set the fixed parameters and initial values
	fpd=['dustem_plugin_phangs_stellar_continuum_2']   ;This is E(B-V)
	fiv=[reddening] ;This is E(B-V) value
	Nage=n_elements(age_values)
	NZ=n_elements(metalicity_values)
	IF do_normalize EQ 1 THEN BEGIN
		NH_value=la_mul(la_mean(use_NHmap[*all_seds_indices[vid]]),10.) ;NH value in 1.e20 H/cm2
		sed.stokesI=sed.stokesI/NH_value   ;normalize SED to 1.e20 H/cm2
		sed.sigmaII=sed.sigmaII/NH_value^2   ;normalize SED to 1.e20 H/cm2
	ENDIF
	FOR i=0L,Nage-1 DO BEGIN
		FOR j=0L,NZ-1 DO BEGIN
			IF weights[i,j] NE 0 AND bins[i,j] NE la_undef() THEN BEGIN
				ij=lonarr(1,2) & ij[0,0]=i & ij[0,1]=j
				ind=ij2index(ij,[Nage,Nz])+1
				fpd=[fpd,'dustem_plugin_phangs_stellar_continuum_'+strtrim(ind[0]+2,2)]   ;+2 is because amplitude,E(B-V) is parameter 0 of pluggin and bins start at zero
				fiv=[fiv,weights[i,j]]   
			ENDIF
		ENDFOR
	ENDFOR

    ;all_params=[pd,fpd]
    ;all_values=[iv,fiv]
    ;try=dustem_parameter_description2type(all_params[0],string_name=string_name)
    ;stop

    ;key=lindgen(n_elements(fpd))
    ;val=fiv
	;toto=dustem_plugin_phangs_stellar_continuum(key=key,val=val)
	;== SET THE OBSERVATIONAL STRUCTURE
	dustem_set_data,m_fit=sed,m_show=sed
	;stop
    IF vid EQ first_vid THEN BEGIN
       ;input_dustem_st=0
	  use_previous_fortran = 0
	  ;==== This is necessary to intialise the phangs pluggin (read the SSPs in particular)
	  ;CAUTION: not sure about the no_reset_plugin_structure keyword ...
	  dustem_init_plugins,pd,fpd=fpd
	  new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $
    						use_previous_fortran=use_previous_fortran,output_dustem_st=output_dustem_st,input_dustem_st=input_dustem_st,dustem_predicted_sed=dustem_predicted_sed,/no_reset_plugin_structure)       
    ENDIF ELSE BEGIN
       use_previous_fortran = 1
       ;Now, This computes the SED from the plugin without runing the Fortran or reading its output (!!)
    	new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $
    						use_previous_fortran=use_previous_fortran,input_dustem_st=output_dustem_st,dustem_predicted_sed=dustem_predicted_sed,/no_reset_plugin_structure)       
    ENDELSE
    ;new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $
   ; 						use_previous_fortran=use_previous_fortran,output_dustem_st=output_dustem_st,input_dustem_st=input_dustem_st,dustem_predicted_sed=dustem_predicted_sed)       
    ;print,vid
    ;print,iv
    ;print,new_iv
    ;print,dustem_predicted_sed
    ;     0.632216     0.641847     0.652599     0.659949     0.602021     0.911603     0.857628     0.873598     0.867334     0.468007     0.432245     0.358588     0.244511     0.179884     0.653411     0.207659
    ;=== get a first guess at parameters (must be done after dustem_set_data)
    ;== compute predicted sed
    ;toto=dustem_compute_sed(new_iv)
    ;print,toto
    ;     0.919570     0.933578     0.949217     0.959909     0.875650      1.32596      1.24750      1.27080      1.26204     0.681323     0.772087     0.554530     0.839953     0.624881      3.04538     0.661734
    ;toto=dustem_predicted_sed*new_iv[3]/iv[3]
    ;print,toto
    ;     0.919570     0.933578     0.949217     0.959909     0.875650      1.32594      1.24743      1.27065      1.26149     0.680648     0.610574     0.517396     0.293335     0.212426     0.680191     0.217372

    ;here !dustem_plugin is not set anymore
    ;as input_dustem_st is set, dustem_activate_plugins is not run
    ;new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $
   ; 						use_previous_fortran=use_previous_fortran,input_dustem_st=output_dustem_st,dustem_predicted_sed=dustem_predicted_sed)       
    ;print,dustem_predicted_sed
    dustem_predicted_sed=dustem_predicted_sed*new_iv[3]/iv[3] ;This takes care of rescaling
    ;stop
    ;cgplot,sed.wave,sed.stokesI,/xlog,/ylog,psym=-4
    ;cgoplot,sed.wave,dustem_predicted_sed,color='red'
    ;stop
    ;== compute difference between predicted and observed SED 
    diff_sed=la_sub(sed.stokesI,dustem_predicted_sed)
    diff_seds[vid,*]=diff_sed
    ;=== update difference maps
    FOR j=0L,Nfilters-1 DO BEGIN
    	map=diff_maps[*,*,j]
    	map[*all_seds_indices[vid]]=diff_seds[vid,j]
    	diff_maps[*,*,j]=map
	ENDFOR
	do_the_next:
ENDFOR

save,diff_maps,href,file=data_dir+'ngc0628_astrosat_minus_prediction.sav',/verb

;NHthreshold=0. & tit=''
;NHthreshold=0.1 & tit='NH > 1e20 H/cm2'
;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'
;vmin=-0.1 & vmax=0.1 & Nbins=50
vmin=-1. & vmax=1. & Nbins=50
;vmin=-5. & vmax=5. & Nbins=100
;vmin=-10. & vmax=10. & Nbins=100
im=diff_maps[*,*,0]
ind=where(im NE la_undef() AND use_NHmap GT NHthreshold)
res0=histogram(im[ind],locations=x0,min=vmin,max=vmax,Nbins=Nbins)
cgplot,x0,res0,/ylog,xtit='(data-model) normalized to 1e21 H/cm2 [MJy/sr]',ytit='N',psym=10,yrange=[min(res0),1.e7],/ysty,color='violet',tit=tit
im=diff_maps[*,*,1]
ind=where(im NE la_undef() AND use_NHmap GT NHthreshold)
res1=histogram(im[ind],locations=x1,min=vmin,max=vmax,Nbins=Nbins)
cgoplot,x1,res1,color='blue',psym=10
im=diff_maps[*,*,2]
ind=where(im NE la_undef() AND use_NHmap GT NHthreshold)
res2=histogram(im[ind],locations=x2,min=vmin,max=vmax,Nbins=Nbins)
cgoplot,x2,res2,color='yellow',psym=10
im=diff_maps[*,*,3]
ind=where(im NE la_undef() AND use_NHmap GT NHthreshold)
res3=histogram(im[ind],locations=x3,min=vmin,max=vmax,Nbins=Nbins)
cgoplot,x3,res3,color='brown',psym=10
im=diff_maps[*,*,4]
ind=where(im NE la_undef() AND use_NHmap GT NHthreshold)
res4=histogram(im[ind],locations=x4,min=vmin,max=vmax,Nbins=Nbins)
cgoplot,x4,res4,color='red',psym=10

cgplot,use_NHmap,diff_maps[*,*,0],/xlog,yr=[-1.,1.],/ysty,psym=3,xr=[1.e-3,1.e2],xtit='NH/1e21',ytit='(data-model) normalized to 1e21 H/cm2 [MJy/sr]'
cgplot,use_NHmap,diff_maps[*,*,0],/xlog,yr=[-10.,10.],/ysty,psym=3,xr=[1.e-3,1.e2],xtit='NH/1e21',ytit='(data-model) normalized to 1e21 H/cm2 [MJy/sr]'

message,'========== Finished computing difference SEDs',/continue
stop


END