PRO check_phangs_ssps_against_uv,show_map=show_map ;check_phangs_ssps_against_uv ;===== This computes a map of prediction by dustem_plugin_phangs_stellar_continuum ;===== constrained on JWST NIR data to astrosat data in the UV for NGC0628. ;===== uses the plugin in a mode where the fortran is not run, but yet, the output files of the fortran are read at each pixel, ;===== which is not efficient (supposed to run in 46hrs .... See check_phangs_ssps_against_uv_fast for an alternative use. window,0,xsize=900,ysize=1000 pdp_define_la_common obp=[1.1,0,1.15,1] source_name='ngc0628' data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/' ;dustem_init,show_plots=show_plots ;needed only for NHCO restore,data_dir+'ngc0628_jwst_images.sav',/verb ;needed for stellar parameters and voronoi bin info restore,data_dir+'ngc0628_muse_images.sav',/verb ;restore,data_dir+'ngc0628_muse_NH.sav',/verb ;restore,data_dir+'ngc0628_astrosat_data.sav',/verb ;stop ;=== restore SEDs in voronoi bins restore,data_dir+'ngc0628_jwst_seds_muse_pixels.sav',/verb jwst_seds=all_seds jwst_filters=(*jwst_seds[0]).filter restore,data_dir+'ngc0628_muse_seds_muse_pixels.sav',/verb muse_seds=all_seds muse_filters=(*muse_seds[0]).filter restore,data_dir+'ngc0628_astrosat_seds_muse_pixels.sav',/verb astrosat_seds=all_seds astrosat_filters=(*astrosat_seds[0]).filter all_filters=[astrosat_filters,muse_filters,jwst_filters] Nvor=max(voronoi_id) Nfilters=n_elements(jwst_filters)+n_elements(astrosat_filters)+n_elements(muse_filters) use_NHmap=NHCO ;used NH map in 1.e21 (from CO) ;galaxy_distance=10. ;MPc use_model='DBP90' ;Example with default keywords uses the DBP90 model ;use_model='MC10' ;Example with default keywords uses the DBP90 model ;use_model='DL07' ;Example with default keywords uses the DBP90 model use_polarization=0 ; initialize Dustemwrap in no polarization mode use_window=2 ; default graphics window number to use for plotting the results use_verbose=0 use_Nitermax=5 ; maximum number of iterations for the fit dustem_define_la_common ;parameters to fit pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction ;'(*!dustem_params).grains(2).mdust_o_mh', $ ;VSGs ;'(*!dustem_params).grains(3).mdust_o_mh', $ ;VSGs ;'(*!dustem_params).grains(4).mdust_o_mh', $ ;BGs 'dustem_plugin_phangs_stellar_continuum_1'] ;stellar continuum amplitude ;iv = [1.6, [2.2e-3, 2.2e-3, 5.7e-3, 5.7e-3],2.e-3] iv = [1.6, [2.2e-3, 2.2e-3],5.e-5] ;iv = [1.6, [2.2e-3, 2.2e-3,2e-3, 2e-3],5.e-5] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) ;llims=[1.e-2,[1.e-4,1.e-4,1.e-4,1.e-4],1.e-10] llims=[1.e-2,[1.e-4,1.e-4],1.e-10] ;llims=[1.e-2,[1.e-4,1.e-4,1.e-4,1.e-4],1.e-10] pd_to_update=['dustem_plugin_phangs_stellar_continuum_1','(*!dustem_params).grains(0).mdust_o_mh'] pd_filter_names=['NIRCAM11','MIRI2'] ;== INITIALISE DUSTEM dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots ;stop !quiet=1 !dustem_nocatch=1 !dustem_verbose=use_verbose IF keyword_set(noobj) THEN !dustem_noobj=1 ;!EXCEPT=2 ; for debugging (not sure hat this is) ;=== INFORMATION TO MAKE THE PLOT ;yr=[1.00e-2,10.0] ; y-axis limits yr=[1.00e-3,100.0] ; y-axis limits ;xr=[0.5,1.00e2] ; x-axis limits xr=[0.1,1.00e2] ; x-axis limits tit='FIT INTENSITY EXAMPLE' ; plot title ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title xtit=textoidl('\lambda (\mum)') ; x-axis title ;voronoi bins for test. Last series are randomly picked one_sed=(*jwst_seds[0])[0] Nfilters=n_elements(astrosat_filters)+n_elements(muse_filters)+n_elements(filters) seds=replicate(one_sed,Nfilters) diff_seds=fltarr(Nvor,Nfilters) diff_maps=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),Nfilters) ;difference maps diff_maps[*]=la_undef() ;imrange=[-1.,1.] ;MJy/sr imrange=[-0.3,0.2] ;MJy/sr t1=systime(0,/sec) first_vid=0L ;first_vid=5960L ;first_vid=12631L ;first_vid=12640L do_normalize=1 FOR vid=first_vid,Nvor-1 DO BEGIN IF vid mod 10 EQ 0 THEN BEGIN t2=systime(0,/sec) delta_t_hr=(t2-t1)/60.^2 ;elapsed time [hr] frac_perc=vid/Nvor*100 delta_t_remain=delta_t_hr/frac_perc*100.-delta_t_hr ;remaining time [hr] message,'done '+strtrim(frac_perc)+' % in '+strtrim(delta_t_hr,2)+' hr remaining '+strtrim(delta_t_remain,2)+' hr',/continue ind=where(diff_maps[*,*,0] NE la_undef(),count) print,la_min(diff_maps[*,*,0]),la_max(diff_maps[*,*,0]),count image_cont20,diff_maps[*,*,0],href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,imrange=imrange ;stop ENDIF ;index=vid ;===== get the SSP weights for the given voronoi bins weights=phangs_binid2weights(st_muse_weights,vid,st_templates,age_values,metalicity_values,reddening=reddening) ;===== extract SED Nf=0L & Nt=Nf+n_elements(astrosat_filters)-1 seds[Nf:Nt]=*astrosat_seds[vid] Nf=Nt+1 & Nt=Nf+n_elements(muse_filters)-1 seds[Nf:Nt]=*muse_seds[vid] Nf=Nt+1 & Nt=Nf+n_elements(jwst_filters)-1 seds[Nf:Nt]=*jwst_seds[vid] sed=seds ind=where(sed.stokesI EQ la_undef(),count) IF count NE 0 THEN goto,do_the_next ;needed to add this because astrosat had 0 variances sed.sigmaii=(5./100.*sed.stokesI)^2 ;===== set the fixed parameters and initial values fpd=['dustem_plugin_phangs_stellar_continuum_2'] ;This is E(B-V) fiv=[reddening] ;This is E(B-V) value Nage=n_elements(age_values) NZ=n_elements(metalicity_values) IF do_normalize EQ 1 THEN BEGIN NH_value=la_mul(la_mean(use_NHmap[*all_seds_indices[vid]]),10.) ;NH value in 1.e20 H/cm2 sed.stokesI=sed.stokesI/NH_value ;normalize SED to 1.e20 H/cm2 sed.sigmaII=sed.sigmaII/NH_value^2 ;normalize SED to 1.e20 H/cm2 ENDIF FOR i=0L,Nage-1 DO BEGIN FOR j=0L,NZ-1 DO BEGIN IF weights[i,j] NE 0 AND bins[i,j] NE la_undef() THEN BEGIN ij=lonarr(1,2) & ij[0,0]=i & ij[0,1]=j ind=ij2index(ij,[Nage,Nz])+1 fpd=[fpd,'dustem_plugin_phangs_stellar_continuum_'+strtrim(ind[0]+2,2)] ;+2 is because amplitude,E(B-V) is parameter 0 of pluggin and bins start at zero fiv=[fiv,weights[i,j]] ENDIF ENDFOR ENDFOR ;all_params=[pd,fpd] ;all_values=[iv,fiv] ;try=dustem_parameter_description2type(all_params[0],string_name=string_name) ;stop ;key=lindgen(n_elements(fpd)) ;val=fiv ;toto=dustem_plugin_phangs_stellar_continuum(key=key,val=val) ;== SET THE OBSERVATIONAL STRUCTURE dustem_set_data,m_fit=sed,m_show=sed ;stop IF vid EQ first_vid THEN BEGIN ;input_dustem_st=0 use_previous_fortran = 0 ;==== This is necessary to intialise the phangs pluggin (read the SSPs in particular) ;CAUTION: not sure about the no_reset_plugin_structure keyword ... dustem_init_plugins,pd,fpd=fpd new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $ use_previous_fortran=use_previous_fortran,output_dustem_st=output_dustem_st,input_dustem_st=input_dustem_st,dustem_predicted_sed=dustem_predicted_sed,/no_reset_plugin_structure) ENDIF ELSE BEGIN use_previous_fortran = 1 ;Now, This computes the SED from the plugin without runing the Fortran or reading its output (!!) new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $ use_previous_fortran=use_previous_fortran,input_dustem_st=output_dustem_st,dustem_predicted_sed=dustem_predicted_sed,/no_reset_plugin_structure) ENDELSE ;new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $ ; use_previous_fortran=use_previous_fortran,output_dustem_st=output_dustem_st,input_dustem_st=input_dustem_st,dustem_predicted_sed=dustem_predicted_sed) ;print,vid ;print,iv ;print,new_iv ;print,dustem_predicted_sed ; 0.632216 0.641847 0.652599 0.659949 0.602021 0.911603 0.857628 0.873598 0.867334 0.468007 0.432245 0.358588 0.244511 0.179884 0.653411 0.207659 ;=== get a first guess at parameters (must be done after dustem_set_data) ;== compute predicted sed ;toto=dustem_compute_sed(new_iv) ;print,toto ; 0.919570 0.933578 0.949217 0.959909 0.875650 1.32596 1.24750 1.27080 1.26204 0.681323 0.772087 0.554530 0.839953 0.624881 3.04538 0.661734 ;toto=dustem_predicted_sed*new_iv[3]/iv[3] ;print,toto ; 0.919570 0.933578 0.949217 0.959909 0.875650 1.32594 1.24743 1.27065 1.26149 0.680648 0.610574 0.517396 0.293335 0.212426 0.680191 0.217372 ;here !dustem_plugin is not set anymore ;as input_dustem_st is set, dustem_activate_plugins is not run ;new_iv=dustem_first_guess(pd,iv,sed,pd_to_update,pd_filter_names,fpd=fpd,fiv=fiv,pol=pol, $ ; use_previous_fortran=use_previous_fortran,input_dustem_st=output_dustem_st,dustem_predicted_sed=dustem_predicted_sed) ;print,dustem_predicted_sed dustem_predicted_sed=dustem_predicted_sed*new_iv[3]/iv[3] ;This takes care of rescaling ;stop ;cgplot,sed.wave,sed.stokesI,/xlog,/ylog,psym=-4 ;cgoplot,sed.wave,dustem_predicted_sed,color='red' ;stop ;== compute difference between predicted and observed SED diff_sed=la_sub(sed.stokesI,dustem_predicted_sed) diff_seds[vid,*]=diff_sed ;=== update difference maps FOR j=0L,Nfilters-1 DO BEGIN map=diff_maps[*,*,j] map[*all_seds_indices[vid]]=diff_seds[vid,j] diff_maps[*,*,j]=map ENDFOR do_the_next: ENDFOR save,diff_maps,href,file=data_dir+'ngc0628_astrosat_minus_prediction.sav',/verb ;NHthreshold=0. & tit='' ;NHthreshold=0.1 & tit='NH > 1e20 H/cm2' ;NHthreshold=0.5 & tit='NH > 5e20 H/cm2' NHthreshold=1. & tit='NH > 1e21 H/cm2' ;NHthreshold=5. & tit='NH > 5e21 H/cm2' ;NHthreshold=10. & tit='NH > 1e22 H/cm2' ;NHthreshold=50. & tit='NH > 5e22 H/cm2' ;vmin=-0.1 & vmax=0.1 & Nbins=50 vmin=-1. & vmax=1. & Nbins=50 ;vmin=-5. & vmax=5. & Nbins=100 ;vmin=-10. & vmax=10. & Nbins=100 im=diff_maps[*,*,0] ind=where(im NE la_undef() AND use_NHmap GT NHthreshold) res0=histogram(im[ind],locations=x0,min=vmin,max=vmax,Nbins=Nbins) cgplot,x0,res0,/ylog,xtit='(data-model) normalized to 1e21 H/cm2 [MJy/sr]',ytit='N',psym=10,yrange=[min(res0),1.e7],/ysty,color='violet',tit=tit im=diff_maps[*,*,1] ind=where(im NE la_undef() AND use_NHmap GT NHthreshold) res1=histogram(im[ind],locations=x1,min=vmin,max=vmax,Nbins=Nbins) cgoplot,x1,res1,color='blue',psym=10 im=diff_maps[*,*,2] ind=where(im NE la_undef() AND use_NHmap GT NHthreshold) res2=histogram(im[ind],locations=x2,min=vmin,max=vmax,Nbins=Nbins) cgoplot,x2,res2,color='yellow',psym=10 im=diff_maps[*,*,3] ind=where(im NE la_undef() AND use_NHmap GT NHthreshold) res3=histogram(im[ind],locations=x3,min=vmin,max=vmax,Nbins=Nbins) cgoplot,x3,res3,color='brown',psym=10 im=diff_maps[*,*,4] ind=where(im NE la_undef() AND use_NHmap GT NHthreshold) res4=histogram(im[ind],locations=x4,min=vmin,max=vmax,Nbins=Nbins) cgoplot,x4,res4,color='red',psym=10 cgplot,use_NHmap,diff_maps[*,*,0],/xlog,yr=[-1.,1.],/ysty,psym=3,xr=[1.e-3,1.e2],xtit='NH/1e21',ytit='(data-model) normalized to 1e21 H/cm2 [MJy/sr]' cgplot,use_NHmap,diff_maps[*,*,0],/xlog,yr=[-10.,10.],/ysty,psym=3,xr=[1.e-3,1.e2],xtit='NH/1e21',ytit='(data-model) normalized to 1e21 H/cm2 [MJy/sr]' message,'========== Finished computing difference SEDs',/continue stop END