PRO check_phangs_ssps_against_uv_fast,show_map=show_map,from_restore=from_restore

;check_phangs_ssps_against_uv_fast

;===== 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

window,1,xsize=700,ysize=700

dustem_define_la_common

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

IF keyword_set(from_restore) THEN goto,from_restore

;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)

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
;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_plugin_phangs_stellar_continuum_1']         ;stellar continuum amplitude
;iv =   [1.6, [2.2e-3, 2.2e-3],5.e-5] 
;Npar=n_elements(pd)

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

;pd_filter_names=['NIRCAM11','MIRI2']
;pd_filter_names=['SDSS4','MIRI2']       ;This normalizes SSPs to Muse filter SDSS4 at ~0.8 mic
pd_filter_names=['NIRCAM11','MIRI2']       ;This normalizes SSPs to NIRCAM12 filter  at ~2 mic 

pd_filter_wav=dustem_filter2wav(pd_filter_names)
pd_filter_index=lonarr(n_elements(pd_filter_names))
FOR i=0L,n_elements(pd_filter_names)-1 DO pd_filter_index[i]=where(all_filters EQ pd_filter_names[i])

;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots
;stop
!quiet=1

;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)

actual_seds=fltarr(Nvor,Nfilters)
predicted_seds=fltarr(Nvor,Nfilters)

all_predicted_specs=ptrarr(Nvor)
all_predicted_specs_normalized=ptrarr(Nvor)

;predicted_maps=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),Nfilters) ;predicted maps
;actual_maps=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),Nfilters)    ;actual maps

diff_maps=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),Nfilters)   ;difference maps

;predicted_maps[*]=la_undef()
;actual_maps[*]=la_undef()

diff_maps[*]=la_undef()
;imrange=[-1.,1.]   ;MJy/sr
imrange=[-0.3,0.2]   ;MJy/sr

;===== Compute seds (pointer)
all_filters=[astrosat_filters,jwst_filters,muse_filters]
Nfilters=n_elements(all_filters)
seds_ptr=ptrarr(Nvor)
nperc=5./100.
FOR ii=0L,Nvor-1 DO BEGIN
	seds_ptr[ii]=ptr_new([*astrosat_seds[ii],*muse_seds[ii],*jwst_seds[ii]])
    ;=== check for 0 variances
	ind=where((*seds_ptr[ii]).sigmaII EQ 0,count)
    IF count NE 0 THEN BEGIN
    	this_sed=*seds_ptr[ii]
    	Ivalues=((*seds_ptr[ii]).stokesI)[ind]
    	sigmaii_values=(Ivalues*nperc)^2
    	this_sed[ind].sigmaii=sigmaii_values
    	seds_ptr[ii]=ptr_new(this_sed)
    ENDIF
ENDFOR

t1=systime(0,/sec)
first_vid=0L

;first_vid=12640L   ;for tests

;do_normalize=1
do_normalize=0

lambir=dustem_get_wavelengths()

;show_each=10
show_each=100

NHvor=fltarr(Nvor)
FOR vid=first_vid,Nvor-1 DO BEGIN
	IF vid mod show_each 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

	ENDIF

	;===== get the SSP weights for the given voronoi bins
	weights=phangs_binid2weights(st_muse_weights,vid,st_templates,age_values,metalicity_values,reddening=reddening)

	sed=*seds_ptr[vid]

	;===== get combined Astro_sat+Muse+JWST SED for that voronoi bin

	ind=where(sed.stokesI EQ la_undef(),count)
	IF count NE 0 THEN goto,do_the_next

	;===== set the fixed parameters and initial values

	amplitude=1.
	reddening=0.     ;This sets Muse reddening to 0. Will be contsrained later
	fpd=phangs_stellar_continuum_plugin_weight2params(weights,parameter_values=fiv,redenning=reddening,/force_include_reddening,amplitude=amplitude,/force_include_amplitude)

   NH_value=la_mul(la_mean(use_NHmap[*all_seds_indices[vid]]),10.) ;NH value in 1.e20 H/cm2
   NHvor[vid]=NH_value
	IF do_normalize EQ 1 THEN BEGIN
		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

	;print,fpd,fiv
	;stop

   Nparams=n_elements(fiv)
   key=intarr(Nparams)
   FOR i=0L,Nparams-1 DO BEGIN
	    toto=dustem_parameter_description2type(fpd[i],string_name=string_name,key=one_key)
	    key[i]=one_key
 	ENDFOR
   spec=dustem_plugin_phangs_stellar_continuum(key=key,val=fiv)

   ;sed_ref=interpol(spec[*,0],lambir,pd_filter_wav[0])

   dustem_predicted_sed=interpol(spec[*,0],lambir,sed.wave)

   ;Those are reference SEDs used for normalization

   sed_ref=dustem_predicted_sed[pd_filter_index[0]]
   sed_goal=(sed.stokesI)[pd_filter_index[0]]

   ;Normalize

   dustem_predicted_sed=dustem_predicted_sed/sed_ref*sed_goal

   all_predicted_specs[vid]=ptr_new(spec[*,0])
   all_predicted_specs_normalized[vid]=ptr_new(spec[*,0]/sed_ref*sed_goal)

	IF vid mod show_each EQ 0 THEN BEGIN
		wset,1
		cgplot,sed.wave,sed.stokesI,/xlog,/ylog,psym=-4
		cgoplot,sed.wave,dustem_predicted_sed,color='red'
 	ENDIF

   actual_seds[vid,*]=sed.stokesI
   predicted_seds[vid,*]=dustem_predicted_sed

	do_the_next:
ENDFOR

;stop

;compute E(B-V) needed to bring prediction to observed astrosat flux
first_vid=0L
ebv=fltarr(Nvor)
astrosat_filters=dustem_filter_names2filters(['F148W','F154W','F169M','F172M','N219M'])
wave_refs=dustem_filter2wav(astrosat_filters)*1.e4 ;AA
which_ebv_filter=1
FOR vid=first_vid,Nvor-1 DO BEGIN
	;flux_ratio=predicted_seds[vid,which_ebv_filter]/actual_seds[vid,which_ebv_filter]  ;This uses only one astrosat filter to derive E(B-V)
	flux_ratio=la_mean(predicted_seds[vid,0:4])/la_mean(actual_seds[vid,0:4]) ;This is an attempt to increase SNR on derived E(B-V)
   ebv[vid]=calz_guess_ebv(wave_refs[which_ebv_filter],flux_ratio,R_v=R_v)
ENDFOR

stop

save,actual_seds,predicted_seds,href,all_seds_indices,do_normalize,all_filters,NHvor,ebv,file=data_dir+'ngc0628_astrosat_voronoi_prediction_fast.sav',/verb

from_restore:

message,'see plot_phangs_ssps_against_uv_fast instead'

restore,data_dir+'ngc0628_astrosat_minus_prediction_fast.sav',/verb
restore,data_dir+'ngc0628_astrosat_voronoi_prediction_fast.sav',/verb

;needed only for NHCO
restore,data_dir+'ngc0628_jwst_images.sav',/verb
use_NHmap=NHCO

;stop

;===== plot histograms

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

xtit='(data-model) [MJy/sr]'
IF do_normalize THEN xtit='(data-model) normalized to 1e21 H/cm2 [MJy/sr]'
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=51

;vmin=-5. & vmax=5. & Nbins=100
;vmin=-10. & vmax=10. & Nbins=100
;im=diff_maps[*,*,0]
ims=actual_seds-predicted_seds

im=ims[*,0]

ind=where(im NE la_undef())
xtit='(data-model) normalized to 1e21 H/cm2 [MJy/sr]'
res=histogram(im[ind],locations=x,min=vmin,max=vmax,Nbins=Nbins)
yrange=[max([min(res),1.]),1.e4]
colors=['violet','blue','yellow','brown','red']

cgplot,x,res,/ylog,xtit=xtit,ytit='N',psym=10,yrange=yrange,/ysty,color='violet',tit=tit,/nodata

FOR k=0L,4 DO BEGIN
	im=ims[*,k]
	ind=where(im NE la_undef())

	res1=histogram(im[ind],locations=x1,min=vmin,max=vmax,Nbins=Nbins)

   cgoplot,x1,res1,color=colors[k],psym=10

ENDFOR

im=ims[*,0]
;ind=where(im NE la_undef() AND use_NHmap GT NHthreshold)
ind=where(im NE la_undef())
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=[max([min(res0),1.]),1.e4],/ysty,color='violet',tit=tit
im=ims[*,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=ims[*,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=ims[*,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=ims[*,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

;same in percent
window,0,xsize=800,ysize=800
!p.multi=[0,1,2]

xtit='(data-model) [%]'
IF do_normalize THEN xtit='(data-model) [%]'
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=-500. & vmax=500. & Nbins=100
;im=diff_maps[*,*,0]
ims=(actual_seds-predicted_seds)/predicted_seds*100.
;ind=where(im NE la_undef() AND use_NHmap GT NHthreshold)
im=ims[*,1]
ind=where(im NE la_undef())
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=[max([min(res0),1.]),1.e7],/ysty,color='violet',tit=tit
im=ims[*,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=ims[*,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=ims[*,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=ims[*,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