PRO fit_phangs_ngc0628_nir_continuum,from_restore=from_restore,brute_force=brute_force, $
									extract_seds=extract_seds,nostop=nostop,muse_stuff=muse_stuff,muse_data=muse_data,astrosat_data=astrosat_data, $
									show_plots=show_plots

;fit_phangs_ngc0628_nir_continuum,/nostop
;fit_phangs_ngc0628_nir_continuum,/muse_stuff
;fit_phangs_ngc0628_nir_continuum,/muse_data         ;do the Muse data "filter" save file
;fit_phangs_ngc0628_nir_continuum,/astrosat_data         ;do the Astrosat data save file
;fit_phangs_ngc0628_nir_continuum,/extract_seds
;fit_phangs_ngc0628_nir_continuum,/from_restore,/show_plots     ;fit seds of some voronoi bins
;fit_phangs_ngc0628_nir_continuum,/from_restore,/brute_force

;fit_phangs_ngc0628_nir_continuum,/muse_data

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

;stop

IF keyword_set(from_restore) THEN goto,from_restore
IF keyword_set(extract_seds) THEN goto,extract_seds
IF keyword_set(muse_stuff) THEN goto,muse_stuff
;IF keyword_set(extract_muse) THEN goto,NH_map
IF keyword_set(muse_data) THEN goto,muse_data
IF keyword_set(astrosat_data) THEN goto,astrosat_data

jwst_stuff:

;========== get JWST reference header
file=data_dir+source_name+'_nircam_lv3_f300m_i2d_anchor_atgauss1.fits' 
d=mrdfits(file,1,h)
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
print,sxpar(h,'CDELT1')
href=cd2astro_header(h)
sxaddpar,href,'EQUINOX',2000.
image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
;stop
Nx=sxpar(href,'NAXIS1')
Ny=sxpar(href,'NAXIS2')

filters_names=['F200W','F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W']
filters=dustem_filter_names2filters(filters_names)
Nfilters=n_elements(filters)
jwst_images=fltarr(Nx,Ny,2,Nfilters)

IF not keyword_set(nostop) THEN stop

;========== JWST stuff
i=0L
file=data_dir+source_name+'_nircam_lv3_f200w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
help,d
;D               DOUBLE    = Array[7289, 9476]
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name

image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

;stop

file=data_dir+source_name+'_nircam_lv3_f300m_i2d_anchor_atgauss1.fits' 
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
help,d
;D               DOUBLE    = Array[3515, 4576]
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_nircam_lv3_f335m_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_nircam_lv3_f360m_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.2,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_miri_lv3_f770w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.2,10],image_color_table='jpbloadct',/silent,tit=tit,off_bar=obp
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_miri_lv3_f1000w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.2,5],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_miri_lv3_f1130w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.3,20],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

file=data_dir+source_name+'_miri_lv3_f2100w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
image_cont20,d,href,/square,imrange=[-0.5,10],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop

;WCO map:
file='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/ngc0628_12m+7m+tp_co21_broad_mom0.fits'
d=readfits(file,h)
;sxaddpar,h,'CTYPE1','RA---TAN'
;sxaddpar,h,'CTYPE2','RA---TAN'
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN WCO=project2(h,d,href,/silent) ELSE WCO=d
fact=4.e20/1.e21
NHCO=la_mul(WCO,fact)   ;NH from CO in 1e21 H/cm2
tit=source_name+' '+'NHCO [1e21 H/cm2]'
image_cont20,NHCO,href,/square,imrange=[-0.5,10],image_color_table='jpbloadct',/silent,tit=tit
IF not keyword_set(nostop) THEN stop


;Invent variances (will have to do better)
;jwst_images[*,*,1,*]=(jwst_images[*,*,0,*]*5./100.)^2   ;assumed intensity variance
perc_error=5./100.
jwst_images[*,*,1,*]=la_power(la_mul(jwst_images[*,*,0,*],perc_error),2)   ;assumed intensity variance
;=== check for null variances
ind=where(jwst_images[*,*,1,*] EQ 0.,count)
IF count NE 0 THEN BEGIN
	message,'There are null variances ...',/continue
	stop
ENDIF

save,jwst_images,filters,href,NHCO,file=data_dir+'ngc0628_jwst_images.sav'

MUSE_data:

;=== this is just to get href
restore,data_dir+'ngc0628_jwst_images.sav',/verb

muse_data_dir='/Volumes/PILOT_FLIGHT1/PHANGS_MUSE/DR2p2/coopt/filterimages/NGC0628/'
;muse_filters=['COUSIN','DUPONT']
;muse_filters=['SDSS2','SDSS3','SDSS4']
muse_filters=dustem_filter_names2filters(['sdss_g','sdss_r','sdss_i'])
Nmuse_filter=n_elements(muse_filters)
muse_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nmuse_filter)
files=muse_data_dir+'NGC0628_'+['SDSS_gcopt','SDSS_rcopt','SDSS_icopt']+'.fits'

;imrange=[-0.5,200]/2000.
imrange=[-1.e-4,0.001]
imrange=[-1.e-2,1.]
muse_undefined=0.
FOR i=0L,Nmuse_filter-1 DO BEGIN
	file=files[i]
	;toto=mrdfits(file,0,hh)
	d=mrdfits(file,1,h)        ;here undefined values are at 0.
	hh=cd2astro_header(h)
	dv=mrdfits(file,2,hv)     ;is this looks like it's a variance
	mask=mrdfits(file,3,hm)
	ind=where(mask EQ 1,count)     ;this looks like a mask
	IF count NE 0 THEN BEGIN
		d[ind]=la_undef()
		dv[ind]=la_undef()
	ENDIF
	print,la_min(d),la_max(d)
	;     -7.37902      17186.9
	factor=dustem_muse_filter_conversion_factor(muse_filters[i],pivot_wav=pivot_wav)  ;goes from Muse filter units to muJy/pix
	omega_pix=(sxpar(hh,'CDELT2')/!radeg)^2    ;simeq 9e-13 sr
	;use_factor=factor/omega_pix*1.e-12         ;simeq 500
	;use_factor=factor/omega_pix*1.e-20         ;simeq 500
	use_factor=(1./factor)/omega_pix*1.e-12         ;simeq 500
	;jpb's estimate:
	cmic=2.99792458e14
	use_factor_jp=1.e-7*1.e4*pivot_wav*(pivot_wav/1.e4)/cmic/omega_pix
	print,factor,omega_pix,use_factor,use_factor_jp
	use_factor=use_factor_jp
	;stop
	;      492.255   9.4017732e-13   5.2357685e-06
    d=la_mul(d,use_factor)  ;Now supposidely in MJy/sr (but in what unit convention ??)
    dv=la_mul(dv,use_factor^2)  ;Now supposidely in (MJy/sr)^2
	print,la_min(d),la_max(d)
    ;-3.8634841e-05     0.089986754
    tit=source_name+' '+muse_filters[i]
	image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit,off_bar_pos=obp,axis_color_table=1
    ;stop
	IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN
		d=project2(h,d,href,/silent,save_it='/tmp/save_project2.sav')
		dv=project2(h,dv,href,/silent,restore_it='/tmp/save_project2.sav')
	ENDIF
	muse_images[*,*,0,i]=d
	muse_images[*,*,1,i]=dv
	;print,filter_name
	image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit,off_bar_pos=obp,axis_color_table=1
ENDFOR

;Invent variances (will have to do better)
;CAUTION: some variances were 0 because some fluxes are 0, had to add abs_variance
;perc_error=5./100.
;abs_variance=(1.e-12)^2
;muse_images[*,*,1,*]=la_add(la_power(la_mul(muse_images[*,*,0,*],perc_error),2),abs_variance)   ;assumed intensity variance
;=== check for null variances
ind=where(muse_images[*,*,1,*] EQ 0.,count)
IF count NE 0 THEN BEGIN
	message,'There are null variances ...',/continue
	stop
	toto=muse_images[*,*,1,*]
	toto[ind]=la_undef()
	muse_images[*,*,1,*]=toto
ENDIF

save,muse_images,muse_filters,href,file=data_dir+'ngc0628_muse_filters_data.sav'
message,'Saved '+data_dir+'ngc0628_muse_filters_data.sav',/info

stop

astrosat_data:
;BUNIT   = 'Angstrom-1 cm-2 erg s-1'
;in fact Angstrom-1 cm-2 erg s-1 PER PIXEL ?!

restore,data_dir+'ngc0628_jwst_images.sav',/verb

astrosat_data_dir='/Volumes/PILOT_FLIGHT1/PHANGS_ASTROSAT/v1p0/release/'
;muse_filters=['COUSIN','DUPONT']
astrosat_filters=dustem_filter_names2filters(['F148W','F154W','F169M','F172M','N219M'])
Nastrosat_filter=n_elements(astrosat_filters)
astrosat_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nastrosat_filter)
files=astrosat_data_dir+'NGC0628_'+['F148','F154','F169','F172','N219']+'_bkg_subtracted_mw_corrected.fits'
;NGC0628_F148_bkg_subtracted_mw_corrected.fits	NGC0628_F169_bkg_subtracted_mw_corrected.fits	NGC0628_N219_bkg_subtracted_mw_corrected.fits
;NGC0628_F154_bkg_subtracted_mw_corrected.fits	NGC0628_F172_bkg_subtracted_mw_corrected.fits
file_save_astrosat=data_dir+'ngc0628_astrosat_data.sav'
wavs_mic=dustem_filter2wav(astrosat_filters)
cmic=3.e14   ;mic/sec
facts=wavs_mic/cmic*(wavs_mic*1.e4)/4./!pi*1.e-7*1.e4*1.e20    ;from Flambda in ergs/cm2/s/AA to MJy/sr
imrange=[-0.5,2]*1e-17    ;erg/s/cm2/AA
imrange=[-0.5,1]   ;MJy/sr

FOR i=0L,Nastrosat_filter-1 DO BEGIN
	file=files[i]
	d=mrdfits(file,0,h)
	steradians=(sxpar(h,'CDELT2')/!radeg)^2   ;pixel in sr
	fact=facts[i]*4.*!pi/steradians
	ind=where(finite(d) NE 1,count)
	IF count NE 0 THEN d[ind]=la_undef()
	d=la_mul(d,fact)
	image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
	IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN d=project2(h,d,href,/silent)
	tit=source_name+' '+astrosat_filters[i]
	image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
	astrosat_images[*,*,0,i]=d
	stop
ENDFOR

;Invent variances (will have to do better)
perc_error=5./100.
astrosat_images[*,*,1,*]=la_power(la_mul(astrosat_images[*,*,0,*],perc_error),2)   ;assumed intensity variance
;=== check for null variances
ind=where(astrosat_images[*,*,1,*] EQ 0.,count)
IF count NE 0 THEN BEGIN
	message,'There are null variances ...',/continue
	stop
ENDIF

save,astrosat_images,astrosat_filters,href,file=file_save_astrosat
message,'Saved '+file_save_astrosat,/info

stop


muse_stuff:

restore,data_dir+'ngc0628_jwst_images.sav',/verb

;========== MUSE stuff
;Muse stellar bins definition
;st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,Nbins=Nbins,Nage=Nage,NZ=NZ)
;print,Nage,NZ,Nbins
;          13           6          78
;st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,/young,Nbins=Nbins,Nage=Nage,NZ=NZ,bins=bins)
;print,Nage,NZ,Nbins
;          17           7         119
;print,la_min(bins),la_max(bins)
;            0          67
st_templates=read_muse_templates_info(age_values=age_values,metalicity_values=metalicity_values,Nbins=Nbins,Nage=Nage,NZ=NZ,bins=bins)
print,Nage,NZ,Nbins
;                    13           6          78
print,la_min(bins),la_max(bins)
;           0          77

;Muse weights
;st_muse_weights=read_muse_phangs_weights(object='NGC0628',/young,bin_numbers=bin_numbers)
;print,n_elements(bin_numbers)
;          68
st_muse_weights=read_muse_phangs_weights(object='NGC0628',bin_numbers=bin_numbers)
print,n_elements(bin_numbers)
;          78

voronoi_id=read_muse_phangs_voronoi_bins('NGC0628',header_in=href,snr_bin=snr_bin,snr_flux=snr_flux,flux=muse_flux)   ;This is the voronoi bin map

image_cont20,voronoi_id,href,/square,/silent,image_color_table='jpbloadct',title='Voronoi Bins'

save,st_templates,st_muse_weights,voronoi_id,age_values,metalicity_values,bins,href,file=data_dir+'ngc0628_muse_images.sav'
message,'saved '+data_dir+'ngc0628_muse_images.sav',/info

stop

extract_seds:

restore,data_dir+'ngc0628_jwst_images.sav',/verb
restore,data_dir+'ngc0628_muse_images.sav',/verb
restore,data_dir+'ngc0628_astrosat_data.sav',/verb
restore,data_dir+'ngc0628_muse_filters_data.sav',/verb

stop

;=== extract and save observed seds in Muse Voronoi bins
extract_all_muse_phangs_seds,source_name,voronoi_id,jwst_images,filters,indices=all_seds_indices,file=data_dir+'ngc0628_jwst_seds_muse_pixels.sav',counts=counts
restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb  ;just to get all_seds_indices
extract_all_muse_phangs_seds,source_name,voronoi_id,muse_images,muse_filters,use_these_indices=all_seds_indices,file=data_dir+'ngc0628_muse_seds_muse_pixels.sav',counts=counts
extract_all_muse_phangs_seds,source_name,voronoi_id,astrosat_images,astrosat_filters,use_these_indices=all_seds_indices,file=data_dir+'ngc0628_astrosat_seds_muse_pixels.sav',counts=counts

;==== The following is not really useful, can be done through SED aggregation later
;==== However, not equivalent to the above, because of undefined pixels in some of the images, and the way this is handled by sed extractor
stop
use_filters=[astrosat_filters,muse_filters,filters]
use_Nfilters=n_elements(use_filters)
use_images=fltarr([sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,use_Nfilters])
use_images[*,*,*,0:n_elements(filters)-1]=jwst_images
use_images[*,*,*,n_elements(filters):n_elements(filters)+n_elements(muse_filters)-1]=muse_images
use_images[*,*,*,n_elements(filters)+n_elements(muse_filters):n_elements(filters)+n_elements(muse_filters)+n_elements(astrosat_filters)-1]=astrosat_images

extract_all_muse_phangs_seds,source_name,voronoi_id,use_images,use_filters,indices=all_seds_indices,file=data_dir+'ngc0628_all_seds_muse_pixels.sav',counts=counts

stop

NH_map:
restore,data_dir+'ngc0628_jwst_images.sav',/verb
restore,data_dir+'ngc0628_muse_images.sav',/verb
restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb

Nvor=max(voronoi_id)

ebv=st_muse_weights.reddening
Rv=3.1
Av_o_NH=1./2.  ;2 mag per 1e21 H/cm2
;Do an Av map
NH_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))  ;in 1e21 H/cm2
FOR vid=0LL,Nvor-1 DO BEGIN
	IF vid mod 100 EQ 0 THEN print,strtrim(vid/Nvor*100.,2)+' %'
	;ind=where(voronoi_id EQ vid)
    Nh_map[*all_seds_indices[vid]]=Rv*ebv[vid]/Av_o_NH
ENDFOR

save,NH_map,file=data_dir+'ngc0628_muse_NH.sav'

from_restore:

dustem_init,show_plots=show_plots
restore,data_dir+'ngc0628_jwst_images.sav',/verb
;% RESTORE: Restored variable: JWST_IMAGES.
;% RESTORE: Restored variable: FILTERS.
;% RESTORE: Restored variable: HREF.
;% RESTORE: Restored variable: NHCO.
restore,data_dir+'ngc0628_muse_images.sav',/verb
;% RESTORE: Restored variable: ST_TEMPLATES.
;% RESTORE: Restored variable: ST_MUSE_WEIGHTS.
;% RESTORE: Restored variable: VORONOI_ID.
;% RESTORE: Restored variable: AGE_VALUES.
;% RESTORE: Restored variable: METALICITY_VALUES.
;% RESTORE: Restored variable: BINS.
;% RESTORE: Restored variable: HREF.
restore,data_dir+'ngc0628_muse_NH.sav',/verb
;% RESTORE: Restored variable: NH_MAP.
restore,data_dir+'ngc0628_astrosat_data.sav',/verb
;% RESTORE: Restored variable: ASTROSAT_IMAGES.
;% RESTORE: Restored variable: ASTROSAT_FILTERS.
;% RESTORE: Restored variable: HREF.
restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb
;% RESTORE: Restored variable: ALL_SEDS.
;% RESTORE: Restored variable: ALL_SEDS_INDICES.
jwst_seds=all_seds
jwst_filters=(*jwst_seds[0]).filter
restore,data_dir+'ngc0628_muse_seds_muse_pixels.sav',/verb   ;contains the Muse SEDs in voronoi bins
;% RESTORE: Restored variable: ALL_SEDS.
;% RESTORE: Restored variable: ALL_SEDS_INDICES.
muse_seds=all_seds
muse_filters=(*muse_seds[0]).filter
restore,data_dir+'ngc0628_astrosat_seds_muse_pixels.sav',/verb
;% RESTORE: Restored variable: ALL_SEDS.
;% RESTORE: Restored variable: ALL_SEDS_INDICES.
astrosat_seds=all_seds
astrosat_filters=(*astrosat_seds[0]).filter

Nvor=max(voronoi_id)
all_filters=[jwst_filters,astrosat_filters,muse_filters]
Nfilters=n_elements(all_filters)
all_seds=fltarr(Nfilters,Nvor)

;===== Compute seds (pointer)
one_sed=(*jwst_seds[0])
all_filters=[jwst_filters,astrosat_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

;stop
;FOR i=0L,Nvor-1 DO BEGIN
;	(*muse_seds[i]).stokesI=((*muse_seds[i]).stokesI)/242315.*1.e8
;ENDFOR
FOR i=0L,Nvor-1 DO all_seds[*,i]=[(*astrosat_seds[i]).stokesI,((*muse_seds[i]).stokesI),(*jwst_seds[i]).stokesI]
;look for maximum astrosat value
ind=where(all_seds[2,*] EQ max(all_seds[2,*]))
print,ind
;       39665
print,all_seds[*,ind]
;      2.35202      2.75261      3.46727      3.68084      1.62485      1054.94      927.261      404.267      10.5454     -32768.0     -32768.0     -32768.0     -32768.0     -32768.0     -32768.0     -32768.0
ind=where(all_seds[2,*] GT 3.*la_mean(all_seds[2,*]))
print,ind[100]
;4699
print,all_seds[*,ind[100]]
;    0.0954950     0.126592     0.144323     0.120900     0.114967      68.6855      52.6373      38.7998      1.15606     0.638490     0.606459     0.529654      2.30200     0.782402      2.51761      1.86948
ind=where(all_seds[2,*]/all_seds[10,*] GT 1.)
print,ind[0]
;12306
print,all_seds[*,ind[0]]
;     0.323970     0.441198     0.591049     0.494629     0.182164      67.2941      50.9925      33.8923     0.967236     0.531406     0.498163     0.374723      1.39314     0.683776      1.90688      1.85371
print,ind[3]
print,all_seds[*,ind[3]]
;       12631
;     0.532320     0.654246     0.788243     0.643166     0.604815      85.0280      61.7483      41.8740      1.29023     0.754277     0.714327     0.569998      1.49024      1.07281      2.24375      1.84148
print,ind[4]
print,all_seds[*,ind[4]]
;       12640
;     0.421350     0.600934     0.639170     0.736419     0.392632      102.730      74.2435      50.7555      1.26155     0.704108     0.633040     0.518208      1.23460     0.844332      1.82023      2.09209
ind=where(all_seds[2,*]/all_seds[10,*] GT 0.6 and all_seds[2,*]/all_seds[10,*] LT 1.)
print,ind[0]
print,all_seds[*,ind[0]]
;        5960
;    0.0603282     0.103469     0.260063     0.171956    0.0811385      53.2903      43.1425      29.1814     0.689211     0.402946     0.368217     0.300235      1.35131     0.478370      1.73637      1.34104
print,ind[20]
print,all_seds[*,ind[20]]
;       13439
;     0.625768     0.544996     0.378714     0.838214     0.591592      92.5699      65.3618      41.9704      1.00230     0.566440     0.614708     0.448454      2.74491      1.15966      3.66431      2.25899

;stop

IF keyword_set(brute_force) THEN goto,brute_force

;toto=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars,info_only=info_only,info_st=info_st)
;stop

;==== SFH ?? (Not sure what this is)
;file='/Volumes/PILOT_FLIGHT1/PHANGS_MUSE/DR2p2/coopt/MUSEDAP/fiducial/NGC0628-0.92asec_table_SFH.fits'
;toto=mrdfits(file,0,h) ;This is ?
;st_sfh=mrdfits(file,1,h) ;This is ?
;help,st_sfh,/str
;** Structure <478ef448>, 10 tags, length=72, data length=68, refs=1:
;   ID              LONG                 0
;   BIN_ID          LONG             47068
;   X               DOUBLE           121.80000
;   Y               DOUBLE          0.60000000
;   FLUX            DOUBLE           24.835750
;   SNR             DOUBLE           6.8757743
;   XBIN            DOUBLE           121.99930
;   YBIN            DOUBLE           2.9014085
;   SNRBIN          DOUBLE           103.86098
;   NSPAX           LONG               284
;print,minmax(st_sfh.id)
;           0        1425
;print,minmax(st_sfh.bin_id)
;       47028       47068


;vid=2000L ;selected voronoi ID
;vid=3000L ;selected voronoi ID
;vid=10000L
;vid=15000L ;selected voronoi ID
;vid=20000L ;selected voronoi ID
;vid=25000L
;vid=30000L ;selected voronoi ID
vid=35000L ;selected voronoi ID
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_plugin_phangs_stellar_continuum_1']         ;stellar continuum amplitude
iv =   [1.6, [2.2e-3, 2.2e-3],5.e-5]
llims=[1.e-2,[1.e-4,1.e-4],1.e-10]
;This allows fitting E(B-V)
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
     'dustem_plugin_phangs_stellar_continuum_2']         ;E(B-V)
iv =   [1.6, [2.2e-3, 2.2e-3],5.e-5,0.01]
llims=[1.e-2,[1.e-4,1.e-4],1.e-10,0.]
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)

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']

;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots
;stop
!dustem_nocatch=1
!dustem_verbose=use_verbose
IF keyword_set(noobj) THEN !dustem_noobj=1
!EXCEPT=2 ; for debugging
;=== 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
vids=[12640,5960,13439,12631,12306,4699,[2000,3000,10000,15000,20000,25000,30000L,35000,400000,450000]]
Nvids=n_elements(vids)
one_sed=(*jwst_seds[0])[0]
Nfilters=n_elements(astrosat_filters)+n_elements(muse_filters)+n_elements(filters)
seds=replicate(one_sed,Nfilters)
FOR ii=0L,Nvids-1 DO BEGIN
	vid=vids[ii]
	;===== get 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]
	;===== fixed parameters
	;fpd=phangs_stellar_continuum_plugin_weight2params(weights,parameter_values=fiv,redenning=reddening,/force_include_reddening)
	fpd=phangs_stellar_continuum_plugin_weight2params(weights,parameter_values=fiv)
	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
	;== SET THE OBSERVATIONAL STRUCTURE
	dustem_set_data,m_fit=sed,m_show=sed
    ;get a first guess at parameters (must be done after dustem_set_data)
	new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol)
	;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE
	;== ADJUSTED DURING THE FIT
	dustem_init_params,use_model,pd,new_iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims
    tit='Voronoi '+strtrim(vid,2)
	;===  RUN THE FIT
	;stop
	t1=systime(0,/sec)
	res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
                      ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $
                      ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
                      ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plots)
	t2=systime(0,/sec)
    IF not keyword_set(nostop) THEN stop
ENDFOR

message,'========== Finished individual fitting',/continue
stop

;==== Brute force
brute_force:

;stop

;==== This is a test to brute-force fit PAH abudance oand G0

Nvor=max(voronoi_id)

;ebv=st_muse_weights.reddening
;Rv=3.1
;Av_o_NH=1./2.  ;2 mag per 1e21 H/cm2
;;Do an Av map
;NH_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))  ;in 1e21 H/cm2
;FOR vid=0LL,Nvor-1 DO BEGIN
;	IF vid mod 100 EQ 0 THEN print,strtrim(vid/Nvor*100.,2)+' %'
;	ind=where(voronoi_id EQ vid)
;    Nh_map[ind]=Rv*ebv[vid]/Av_o_NH
;ENDFOR

;stop

dir=!dustem_wrap_soft_dir+'/Grids/'
table_name=dir+'DBP90_JWST_G0_YPAH_YVSG_4Phangs.fits'
;table_name=dir+'TEST_DBP90_JWST_G0_YPAH_YVSG_4Phangs.fits'
;Nvor=20000 ;for test

GOs=fltarr(Nvor)
Ypahs=fltarr(Nvor)
YVSGs=fltarr(Nvor) 
facts=fltarr(Nvor)
chi2s=fltarr(Nvor)
rchi2s=fltarr(Nvor)
dGOs=fltarr(Nvor)
dYpahs=fltarr(Nvor)
dYVSGs=fltarr(Nvor) 

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

;==== This is all the filters in a given SED
all_filters=(*seds_ptr[0]).filter
;==== filters used for brute force fit
use_filters_names=['F360M','F0770W','F1000W','F1130W','F2100W']
use_filters=dustem_filter_names2filters(use_filters_names)
use_Nfilters=n_elements(use_filters)
ind_filters=[-1]
FOR i=0L,use_Nfilters-1 DO BEGIN
	indd=where(all_filters EQ use_filters[i],count)
	IF count NE 0 THEN ind_filters=[ind_filters,indd]
ENDFOR
ind_filters=ind_filters[1:*]

;show_sed=1
show_sed=0
FOR vid=0LL,Nvor-1 DO BEGIN
   IF vid mod 10 EQ 0 THEN BEGIN
   	 message,'Doing sed '+strtrim(vid,2)+' '+strtrim(1.*vid/Nvor*100.,2)+' %',/continue
   ENDIF
   ;sed=all_seds[ind_filters,vid]
   sed=*seds_ptr[vid]
   ;=== restrict sed to requested filters
   sed=sed[ind_filters]
   ;=== normalize to NH=1.e20 H/cm2
   NH_value=la_mul(la_mean(use_NHmap[*all_seds_indices[vid]]),10.) ;NH value in 1.e20 H/cm2
   sed.stokesI=la_div(sed.stokesI,NH_value)
   sed.sigmaII=la_div(sed.sigmaII,NH_value^2)
   ;index=where(voronoi_id EQ vid,count)
   ;sed=dustem_sed_extractor(jwst_images,index,filters,/total_intensity_only)
   ;stop
   params=dustem_brute_force_fit(sed,table_name,use_filters,fact=fact,chi2=chi2,/normalize,rchi2=rchi2,show_sed=show_sed $
   								,params_hit=params_hit,params_uncertainties=params_uncertainties,params_min=params_min,params_max=params_max)
   GOs[vid]=params[0]
   Ypahs[vid]=params[1]/fact ;extenssive quantities must be devided by normalization factor
   Yvsgs[vid]=params[2]/fact
   facts[vid]=fact
   chi2s[vid]=chi2
   rchi2s[vid]=rchi2
   dGOs[vid]=params_uncertainties[0]
   dYpahs[vid]=params_uncertainties[1]/fact ;extenssive quantities must be devided by normalization factor
   dYvsgs[vid]=params_uncertainties[2]/fact
   print,params[0],params[1],params[2],fact,chi2
ENDFOR

stop

;save,GOs,Ypahs,Yvsgs,file=dir+source_name+'DBP90_JWST_G0_YPAH_YVSG.fits'
save,GOs,Ypahs,Yvsgs,file=dir+source_name+'TEST_DBP90_JWST_G0_YPAH_YVSG.fits'

;stop
;=== make maps of parameters
chi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
rchi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
GO_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
Ypah_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
YVSG_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
dGO_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
dYpah_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
dYVSG_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
fact_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))
FOR vid=0LL,Nvor-1 DO BEGIN
	IF vid mod 100 EQ 0 THEN print,(1.*vid)/Nvor*100.
	;ind=where(voronoi_id EQ vid)
    GO_map[*all_seds_indices[vid]]=GOs[vid]
    Ypah_map[*all_seds_indices[vid]]=Ypahs[vid]
    YVSG_map[*all_seds_indices[vid]]=Yvsgs[vid]
    dGO_map[*all_seds_indices[vid]]=dGOs[vid]
    dYpah_map[*all_seds_indices[vid]]=dYpahs[vid]
    dYVSG_map[*all_seds_indices[vid]]=dYvsgs[vid]
    chi2_map[*all_seds_indices[vid]]=chi2s[vid]
    rchi2_map[*all_seds_indices[vid]]=rchi2s[vid]
    fact_map[*all_seds_indices[vid]]=facts[vid]
ENDFOR

;stop
jwst_4_coutours=fltarr((size(jwst_images))[1],(size(jwst_images))[2],2)
jwst_4_coutours[*,*,0]=jwst_images[*,*,0,5]
jwst_4_coutours[*,*,1]=jwst_images[*,*,0,5]
levs=create_contour_st(jwst_4_coutours,lev1=[1,5])
levs.(0).color=100
levs.(1).color=200
levs.(1).thick=2

obp=[1.1,0.,1.15,1]
image_cont20,jwst_4_coutours[*,*,0],href,/square,imrange=[-1.,10],image_color_table='jpbloadct',/silent,title='',off_bar=obp,axis_color_table=1,levels=levs

window,0,xsize=800,ysize=1000
aa=alog10(GO_map)
wset,0 & image_cont20,GO_map,href,/square,imrange=[-1.,100],image_color_table='jpbloadct',/silent,title='G0',off_bar=obp,axis_color_table=1
wset,0 & image_cont20,aa,href,/square,imrange=[0.5,2.5],image_color_table='jpbloadct',/silent,title='G0',off_bar=obp,axis_color_table=1,levels=levs
window,1,xsize=800,ysize=1000
wset,1 & image_cont20,Ypah_map,href,/square,imrange=[-0.005,0.5],image_color_table='jpbloadct',/silent,title='Ypah',off_bar=obp,axis_color_table=1,levels=levs
window,2,xsize=800,ysize=1000
wset,2 & image_cont20,Yvsg_map,href,/square,imrange=[-0.005,0.5],image_color_table='jpbloadct',/silent,title='Yvsg',off_bar=obp,axis_color_table=1,levels=levs
window,3,xsize=800,ysize=1000
wset,3 & image_cont20,alog10(chi2_map),href,/square,imrange=[2,6],image_color_table='jpbloadct',/silent,title='chi2',off_bar=obp,axis_color_table=1,levels=levs
window,4,xsize=800,ysize=1000
image_cont20,alog10(fact_map),href,/square,imrange=[-1,5],image_color_table='jpbloadct',/silent,title='fact',off_bar=obp,axis_color_table=1,levels=levs

stop

wset,1 & image_cont20,Ypah_map,href,/square,imrange=[-0.0001,0.005],image_color_table='jpbloadct',/silent,title='Ypah',off_bar=obp,axis_color_table=1
wset,1 & image_cont20,Ypah_map,href,/square,imrange=[-0.0001,0.1],image_color_table='jpbloadct',/silent,title='Ypah',off_bar=obp,axis_color_table=1

wset,2 & image_cont20,Yvsg_map,href,/square,imrange=[-0.001,0.015],image_color_table='jpbloadct',/silent,title='Yvsg',off_bar=obp,axis_color_table=1

window,5
!p.multi=[0,3,2]
histGO=histogram(GOs,locations=GOv,Nbins=30)
cgplot,GOv,histGO,psym=10,/ylog,xtit='G0',yrange=[1,1.e5]
histYpah=histogram(Ypahs,locations=Ypahsv,Nbins=10)
cgplot,Ypahsv,histYpah,psym=10,/ylog,xtit='Ypah'
histYvsg=histogram(Yvsgs,locations=Yvsgsv,Nbins=10)
cgplot,Yvsgsv,histYvsg,psym=10,/ylog,xtit='Yvsg'
histchi2=histogram(chi2s,locations=chisv,Nbins=10)
cgplot,chisv,histchi2,psym=10,/ylog,xtit='Chi2'
histfacts=histogram(facts,locations=factsv,Nbins=10)
cgplot,factsv,histfacts,psym=10,/ylog,xtit='facts'


save,GO_map,Ypah_map,Yvsg_map,href,chi2_map,fact_map,GOs,Ypahs,Yvsgs,chi2s,facts,file=data_dir+'ngc0628_brute_force_fit_results.sav'

restore,data_dir+'ngc0628_brute_force_fit_results.sav',/verb


Nid=la_max(voronoi_id)
ids=fltarr(Nid)
snrs=fltarr(Nid)
counts=fltarr(Nid)
FOR i=0L,Nid-1 DO BEGIN
	ind=where(voronoi_id EQ i,count)
	ids[i]=i
	snrs[i]=snr[ind[0]]
	counts[i]=count
ENDFOR

cgplot,ids,snrs
cgplot,ids,counts
cgplot,counts,snrs

ind1=where(snr EQ max(snr))
ind=where(voronoi_id EQ voronoi_id[ind1],count)
print,count

ind=where(voronoi_id EQ 0,count)


file='/Volumes/PILOT_FLIGHT1/PHANGS_MUSE/DR2p2/coopt/MUSEDAP/fiducial/NGC0628-0.92asec_ppxf_SFH.fits'
st=mrdfits(file,1,h)
help,st,/str


END