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
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
	;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)
	;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
	;
	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
       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=dustem_st,input_dustem_st=input_dustem_st,dustem_predicted_sed=dustem_predicted_sed)
       ;pb after this !dustem_plugin)[i].spec not set
    ENDIF ELSE BEGIN
       use_previous_fortran = 1
       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=dustem_st,dustem_predicted_sed=dustem_predicted_sed)       
    ENDELSE
    ;stop
    ;=== get a first guess at parameters (must be done after dustem_set_data)
    ;== compute predicted sed
    ;dustem_predicted_sed=dustem_compute_sed(new_iv,st=st)
    ;== 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
ENDFOR

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


END