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,grid_brute_force=grid_brute_force ;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,/from_restore,/grid_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/' data_dir=!phangs_data_dir+'/ISRF/WORK/' grids_data_dir=!phangs_data_dir+'/ISRF/GRIDS/' ;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 IF keyword_set(grid_brute_force) THEN goto,grid_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 and G0 Using the Mathis ISRF 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/' ;==== This is to use the Mathis Field table_name=grids_data_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 stop ;==== Brute force fit using ISRF classes grid_brute_force: ;stop Nvor=max(voronoi_id) G0s=fltarr(Nvor)+la_undef() Ypahs=fltarr(Nvor)+la_undef() YVSGs=fltarr(Nvor)+la_undef() facts=fltarr(Nvor)+la_undef() chi2s=fltarr(Nvor)+la_undef() rchi2s=fltarr(Nvor)+la_undef() dG0s=fltarr(Nvor)+la_undef() dYpahs=fltarr(Nvor)+la_undef() dYVSGs=fltarr(Nvor)+la_undef() 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 model='DBP90' ;dir_grids=!dustem_wrap_soft_dir+'/Grids/' file=data_dir+'/ngc0682_isrf_classes_one_ratio.sav' restore,file,/verb ;% RESTORE: Restored variable: CLASSES. ;% RESTORE: Restored variable: VOR_CLASS. ;% RESTORE: Restored variable: MAP_CLASSES. ind=where(vor_class NE 0) class_min=min(vor_class[ind]) class_max=max(vor_class[ind]) restore,data_dir+'ngc0628_isrf_min_prediction.sav',/verb ;% RESTORE: Restored variable: ISRFS. ;% RESTORE: Restored variable: OBJECT_DISTANCE. ;% RESTORE: Restored variable: OBJECT_THICKNESS. ;% RESTORE: Restored variable: SOURCE_NAME. ;For test class_min=10 class_max=15 bidon=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths) fixed_parameters_description=['dustem_plugin_phangs_class_isrf_2'] ;This is the ISRF G0 factor file=!dustem_soft_dir+'/data/ISRF_MATHIS.DAT' ;file='/Users/jpb/Soft_Libraries/dustem_fortran/data/ISRF.DAT' readcol,file,Mathis_wavs,Mathis_ISRF Mathis_ISRF=interpol(Mathis_ISRF,Mathis_wavs,isrf_wavelengths) Mathis_1mic=interpol(Mathis_ISRF,isrf_wavelengths,1.0) FOR isrf_class=class_min,class_max DO BEGIN ;==== This is to use the grid for this ISRF class isrf_class_str='_isrfclass'+strtrim(isrf_class,2) table_name=grids_data_dir+model+'_MuseISRF_JWST_G0_YPAH_YVSG_4Phangs'+isrf_class_str+'.fits' ;select Voronoi bins with that ISRF class ind_class=where(vor_class EQ isrf_class,Nvor_class) IF Nvor_class NE 0 THEN BEGIN FOR vvid=0LL,Nvor_class-1 DO BEGIN IF vvid mod 100 EQ 0 THEN BEGIN message,'isrf#'+strtrim(isrf_class,2)+' Doing sed '+strtrim(vvid,2)+' '+strtrim(1.*vvid/Nvor_class*100.,2)+' %',/continue ENDIF ;=== Get the G0 value for that vid vid=ind_class[vvid] G0=interpol(ISRFs[*,vid],isrf_wavelengths,1.0)/Mathis_1mic fixed_parameters_values=[G0] 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) ind=where(sed.stokesI EQ la_undef(),count) IF count EQ 0 THEN BEGIN 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 $ ,fixed_parameters_description=fixed_parameters_description,fixed_parameters_values=fixed_parameters_values) ;This would be if there was no fixed_parameters ;G0s[vid]=params[0] ;Ypahs[vid]=params[1]/fact ;extenssive quantities must be divided by normalization factor ;Yvsgs[vid]=params[2]/fact G0s[vid]=G0 Ypahs[vid]=params[0]/fact ;extenssive quantities must be divided by normalization factor Yvsgs[vid]=params[1]/fact facts[vid]=fact chi2s[vid]=chi2 rchi2s[vid]=rchi2 ;This would be if there was no fixed_parameters ;dG0s[vid]=params_uncertainties[0] ;dYpahs[vid]=params_uncertainties[1]/fact ;extenssive quantities must be devided by normalization factor ;dYvsgs[vid]=params_uncertainties[2]/fact dG0s[vid]=0. ;should use difference between fixed param and nearest table neighbor dYpahs[vid]=params_uncertainties[0]/fact ;extenssive quantities must be devided by normalization factor dYvsgs[vid]=params_uncertainties[1]/fact ;print,params[0],params[1],params[2],fact,chi2 ENDIF ENDFOR ENDIF ELSE BEGIN message,'No Voronoi bins found with ISRF class '+strtrim(isrf_class,2),/continue ENDELSE ENDFOR ;stop ;save,GOs,Ypahs,Yvsgs,file=dir+source_name+'DBP90_JWST_G0_YPAH_YVSG.fits' save,G0s,Ypahs,Yvsgs,facts,chi2s,rchi2s,file=data_dir+source_name+'_DBP90_JWST_G0_YPAH_YVSG_ISRFclasses.fits' ;Nvor=max(voronoi_id) ;stop ;=== make maps of parameters chi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() rchi2_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() G0_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() Ypah_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() YVSG_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() dG0_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() dYpah_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() dYVSG_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() fact_map=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef() 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) G0_map[*all_seds_indices[vid]]=G0s[vid] Ypah_map[*all_seds_indices[vid]]=Ypahs[vid] YVSG_map[*all_seds_indices[vid]]=Yvsgs[vid] dG0_map[*all_seds_indices[vid]]=dG0s[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 ;FOR vid=0LL,Nvor-1 DO GO_map[*all_seds_indices[vid]]=GOs[vid] ;FOR vid=0LL,Nvor-1 DO Ypah_map[*all_seds_indices[vid]]=Ypahs[vid] ;FOR vid=0LL,Nvor-1 DO YVSG_map[*all_seds_indices[vid]]=Yvsgs[vid] ;FOR vid=0LL,Nvor-1 DO chi2_map[*all_seds_indices[vid]]=chi2s[vid] stop ;G0 histogram win=0L window,win,xsize=800,ysize=800 & win=win+1 ind=where(G0s NE la_undef()) res=histogram(G0s[ind],locations=xv) cgplot,xv,res,psym=10,title='G0s histogram',xtit='G0',ytit='Number',/ylog,yrange=[1,max(res)],/ysty ind=where(Ypahs NE la_undef()) res=histogram(Ypahs[ind],locations=xv) cgplot,xv,res,psym=10,title='Ypahs histogram',xtit='Ypah',ytit='Number',/ylog,yrange=[1,max(res)],/ysty window,win,xsize=800,ysize=900 & win=win+1 obp=[1.1,0.,1.15,1] image_cont20,G0_map,Href,/square,imrange=[-2,40],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='G0' image_cont20,la_log10(G0_map),Href,/square,imrange=[-0.2,2],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='G0' image_cont20,la_log10(Ypah_map),Href,/square,imrange=[-2.,6],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='log10(Ypah)' image_cont20,la_log10(YVSG_map),Href,/square,imrange=[-2.,6],axis_color_table=1,image_color_table='jpbloadct',/silent,off_bar=obp,title='log10(Yvsg)' image_cont20,chi2_map,Href,/square,imrange=[0.01,100]*1.e4,image_color_table='jpbloadct',/silent ;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=[-2.,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 END