PRO dustem_fit_sed_ext_pol_example,model=model $ ,sed=sed $ ,itermax=itermax $ ,postscript=postscript $ ,png=png $ ,save=save $ ,restore=restore $ ,help=help IF keyword_set(help) THEN BEGIN doc_library,'dustem_fit_sed_readme' goto,the_end END IF keyword_set(model) THEN BEGIN use_model=strupcase(model) ENDIF ELSE BEGIN use_model='G17_MODELA' ;Default is last dustem model ENDELSE use_polarization=0. ;default is no polarization in models use_window=2 ; default graphics window number to use for plotting the results ;=== Set the (model-dependent) parameters that you want to fit ;=== Refer to the DustEM and DustEMWrap userguides for an explanation ; of the different grain types CASE use_model OF 'DBP90':BEGIN 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', $ ;amCBEx 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0, 4.3e-4, 4.7e-4,6.4e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'DL01':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'WD01_RV5p5B':BEGIN ;; ***COMMENT AH*** ;; we need to implement this, or remove message, 'WD01 model not yet implemented in DustEMWrap',/info ;; ***END COMMENT AH*** END 'DL07':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001];,10,1.,10.,1.] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'MC10':BEGIN pd = [ $ ;'(*!dustem_params).gas.G0', $ ;G0 ; 'dustem_plugin_continuum_2', $ ;intensity of NIR continuum ; 'dustem_plugin_synchrotron_2',$ '(*!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', $ ;amCBEx '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(4).mdust_o_mh' $ ;aSilx ] ;initial parameter values for parameters to be fitted rv = [ $ ; 1.0, $ ; 0.002, $ ;intensity of NIR continuum ; 0.01,$ 7.8e-4, $ ;mass fraction of PAH0 7.8e-4, $ ;mass fraction of PAH1 1.65e-4, $ ;mass fraction of amCBEx 1.45e-3, $ ;mass fraction of amCBEx 7.8e-3 $ ;mass fraction of aSilx ] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) ;=== Fixed parameters ;parameter description of parameters to be set to a non-default value ; fpd=[ $ ; '(*!dustem_params).gas.G0' , $ ;multiplicative factor to total ISRF ; 'dustem_plugin_continuum_2' $ ;intensity of NIR continuum ; ] ; ;initial parameter values for fixed parameters ; fiv=[ $ ; 1. , $ ;multiplicative factor to total ISRF ; 3.e-3$ ;intensity of NIR continuum ; ] ;; ***END COMMENT AH*** END 'J13':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'G17_MODELA':BEGIN 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',$ ;aSil 'dustem_plugin_modify_dust_polx_2' ] ;'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [5.4e-4, 5.4e-4,1.8e-4,50] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) use_polarization=0 ;=== Fixed parameters ;parameter description of parameters to be set to a non-default value fpd=[ $ '(*!dustem_params).gas.G0' , $ ;multiplicative factor to total ISRF 'dustem_plugin_continuum_2' $ ;intensity of NIR continuum ] ;initial parameter values for fixed parameters fiv=[ $ 1. , $ ;multiplicative factor to total ISRF 3.e-3$ ;intensity of NIR continuum ] END 'G17_MODELB':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) use_polarization=1 END 'G17_MODELC':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) use_polarization=1 END 'G17_MODELD':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) use_polarization=1 END 'ELSE':BEGIN message,'model '+model+' unknown',/continue message,'Known models are MC10,DBP90,DL01,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue stop END ENDCASE ;== INITIALISE DUSTEM ;; ***COMMENT AH*** ;; do we need the use_polarization key word activated here? --> ;; len(kwords) must match Npar_dust (i.e. without plugins), otherwise ;; fails in dustem_read_grain -- should set this within the above case statement? ;; IC: yes good idea. Maybe by default set 'kwords' to the original dust keywords and let the user modify kwords(i) where i is the index of the dust species. dustem_init,model=use_model,/pol;,kwords=['logn','plaw','plaw'];,polarization=use_polarization ;; ***END COMMENT AH*** !dustem_nocatch=1 !dustem_verbose=1 !dustem_show_plot=1 !EXCEPT=2 ;so I can locate the plotting error... ;=== READ EXAMPLE DATA: EXTINCTION ;NB:HERE WE ARE NOT FITTING THE DATA. WE'RE JUST USING THE EXISTING FILTERS. dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' file=dir+'Mathis1990_DISM.xcat' ;We should add the serkowski xcat's content into the MATHIS one right? ;if keyword_set(ext) then file=ext ; not really needed for the tests ext=read_xcat(file,/silent) dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' file=dir+'example_SED_2.xcat' sed=read_xcat(file,/silent) ;=== GENERATE EXAMPLE DATA: EXTINCTION ;filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3'] ;Nfilt=n_elements(filters) ;ext=dustem_initialize_ext(Nfilt) ;ext.filter=replicate('SPECTRUM',Nfilt) ;Using these filters to set the wavelengths ;ext.wave=dustem_filter2wav(filters) ;ext.instru=replicate('EXTINCTION', Nfilt);don't really know what to put here ;dustem_filter2instru(filters) ;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro ext.EXT_I=1. for i=4l,n_tags(ext)-1 do begin ext.(i) = ext.EXT_I/1E13 endfor ;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro sed.StokesI=1. for i=4l,n_tags(sed)-1 do begin sed.(i) = sed.StokesI/1E10 endfor ; ind=where(spec.sigextII EQ 0.,count) ; IF count NE 0 THEN spec[ind].sigextII=(0.2*spec(ind).EXT_I)^2 ; ind=where(spec.instru EQ 'FIRAS',count) ; I do not know what this is for (WHY FIRAS)... maybe an old remnant of the older procedures ; IF count NE 0 THEN spec[ind].sigextII=(0.2*spec(ind).EXT_I)^2 ;###Before merging any xcat files we'll just displace the initial parameter array. ;###Creating fake extinction data: (already replaced iv by rv) ;###Creating observational strcture using xcat file (for simplicity) ;###Initializing the structure tags so that dustem_set_data sets said tags in !dustem_data ; for i=4l,n_tags(spec) -1 do begin ; ; ;initial loop ; ;actually only data pertaining to Q and U need to be set. ; ;but dustem_set_data will take care of the tags of the !dustem_data structure pertaining to other polarization quantities ; ;The content of the xcat file shouldn't not/isn't used after this in the main procedure ; if (tag_names(spec))[i] ne 'SIGEXTII' then begin ; spec.(i) = spec.(3)/1E12 ;just to fill the tags... ; ;This gives really high sigma values ; endif ; ; endfor ;###first call to dustem_set_data so that the !dustem_data tags that are used in the 'compute_' procedures are set ;dustem_set_data,m_fit,m_show,ext,ext dustem_set_data,sed,sed,ext,ext iv = rv+[10.00E-4,10.00E-4,10.00E-4,-30] ;###setting the initial parameter vector ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE ;== ADJUSTED DURING THE FIT dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims ;== INITIALIZE ANY PLUGINS dustem_init_plugins, pd,fpd dustem_init_fixed_params,fpd,fiv ;generation of the fake data (since dustem_set_data takes ONE structure containing all the info) ;!dustem_psi_ext=25 ;using the plugin instead ext.EXT_I = dustem_compute_ext(rv,st) toto = dustem_compute_stokext(rv,st,dustem_qext,dustem_uext) ext.EXT_Q = dustem_qext ext.EXT_U = dustem_uext sed.StokesI = dustem_compute_sed(rv,st) toto = dustem_compute_stokes(rv,st,dustem_qsed,dustem_used) sed.StokesQ = dustem_qsed sed.StokesU = dustem_used ;; ***END COMMENT AH*** ;== SET THE OBSERVATIONAL STRUCTURE dustem_set_data, sed,sed,ext,ext;sed=spec) ; dustem_init_fixed_params,fpd,fiv ;== RUN THE FIT tol=1.e-10 use_Nitermax=30 ;maximum number of iterations. IF keyword_set(itermax) THEN use_Nitermax=itermax xr_x = [0.01,30] yr_x = [1.00E-10,10] xr_m = [1.,5e5] yr_m = [5e-8,1.00e6] tit='Spectral Energy Distribution' ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') t1=systime(0,/sec) res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $ ,/xlog,/ylog,xr=xr,yr=yr,xr_m=xr_m,yr_m=yr_m,xr_x=xr_x,yr_x=yr_x,xtit=xtit,ytit=ytit,title=tit $ ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $ ,errors=errors,chi2=chi2,rchi2=rchi2) t2=systime(0,/sec) ;=== SAVE FIT RESULTS ;file_out='/tmp/DUSTEM_fit_example.sav' if keyword_set(save) then begin dir_sav='./' dustem_save_system_variables,dir_sav+save message,'Saved fit results in '+dir_sav+save,/continue if not keyword_set(restore) then $ restore=dir_sav+save end message,'The fit executed in '+strtrim(t2-t1,2)+' sec',/info ;====================================== ;====You could exit IDL here. The remaining lines of code (essentially ;====plotting the results) would work by returning to this point ;====restoring the IDL output file that was created above. ;====================================== ;file='/tmp/DUSTEM_fit_example.sav' if keyword_set(restore) then begin file=restore dustem_restore_system_variables,file end ;== PLOT THE FIT RESULTS RESTORED FROM .SAV FILE ; tit='Spectral Energy Distribution (Saved)' ; ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; xtit=textoidl('\lambda (\mum)') ; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values) ; chi2=(*!dustem_fit).chi2 ; rchi2=(*!dustem_fit).rchi2 ; ; res=*(*!dustem_fit).current_param_values ; chi2=(*!dustem_fit).chi2 ; rchi2=(*!dustem_fit).rchi2 ; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values) ;=== PLOT FIT RESULTS AND SAVE TO GRAPHICS FILE IF REQUESTED ; window,use_window ; loadct,13 IF keyword_set(postscript) THEN BEGIN ; dir_ps=!dustem_dat+'/Figures/' dir_ps='./' ; force_mkdir,dir_ps set_plot,'PS' ps_file=dir_ps+postscript device,filename=ps_file,/color ENDIF ;dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit ; dustem_sed_plot,*(*!dustem_fit).current_param_values $ ; ,ytit=ytit,xtit=xtit,title=tit $ ; ,yr=yr,xr=xr,/ysty,/xsty $ ; ,res=res,errors=errors,chi2=chi2,rchi2=rchi2 $ ; ,/xlog,/ylog,legend_xpos=legend_xpos,legend_ypos=legend_ypos IF keyword_set(postscript) THEN BEGIN device,/close set_plot,'X' message,'Wrote '+ps_file,/info stop ENDIF IF keyword_set(png) THEN BEGIN ; dir_png=!dustem_dat+'/Figures/' dir_png='./' ; force_mkdir,dir_png file_png=dir_png+png write_png,file_png,tvrd(/true) message,'Wrote '+file_png,/info ENDIF message,'Finished dustem_fit_sed_readme',/info the_end: END