PRO dustem_run_prediction_example, model $ ,filters $ ,wave=wave $ ,next=next $ ,sed_outfile=sed_outfile $ ,ext_outfile=ext_outfile $ ,polarisation=polarisation $ ,show=show $ ,fits_save=fits_save $ ,noobj=noobj $ ,noexit=noexit $ ,help=help ;+ ; NAME: ; dustem_run_prediction_example ; ; PURPOSE: ; This is an example of how to run the DustEM fortran code using ; DustEMWrap and use it to make prediction for fluxes in a set of filters ; ; CATEGORY: ; DustEMWrap, Distributed, High-Level, User Example ; ; CALLING SEQUENCE: ; dustem_run_prediction_example,model,filters,show=show,postscript=postscript,help=help ; ; INPUTS: ; model = specifies the interstellar dust mixture used by ; DustEM. See userguide or dustem_test_model_exists.pro ; for more details about available models in current release. ; filters = specifies the filters for which flux predictions should ; be calculated. The filters will be shown on the plot. ; ; OPTIONAL INPUT PARAMETERS: ; wave = wavelength range specified as [minwave,maxwave] in microns over which ; to construct the extinction curve. The resulting curve ; will be defined on the vector : ; wave_vector=10^(alog10(wave[0])+alog10(wave[1])*findgen(Next)/float(Next)) ; default is [0.01,50] ; Next = number of measurements in the extinction curve. Default is 100. ; OUTPUTS: ; ; OPTIONAL OUTPUT PARAMETERS: ; sed_outfile = name of file to write predicted SED (in .xcat ; format). Default is './dustemwrap_predicted_sed.xcat' ; ext_outfile = name of file to write predicted extinction curve (in .xcat ; format). Default is './dustemwrap_predicted_ext.xcat' ; ; ACCEPTED KEY-WORDS: ; help = if set, print this help ; polarisation = if set, run DustemWrap in polarisation mode ; show = if set, show a plot of the predicted SED and the ; model spectra ; fits_save = if set, save the fit results in a binary ; FITS file. ; noobj = if set, runs with no object-oriented graphics ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; The DustEM fortran code must be installed ; The DustEMWrap IDL code must be installed ; ; PROCEDURE: ; None ; ; EXAMPLES ; dustem_run_prediction_example,'DBP90',['MIRI1','SPIRE3','NIKA21','IRAS1','IRAS2','IRAS3','SPIRE2','SPIRE1','PACS1','WISE3'] ; ; MODIFICATION HISTORY: ; Written JPB Apr-2011 ; Evolution details on the DustEMWrap gitlab. ; See http://dustemwrap.irap.omp.eu/ for FAQ and help. ;- IF keyword_set(help) THEN BEGIN doc_library,'dustem_run_prediction_example' goto,the_end ENDIF exists=dustem_test_model_exists(model) if exists ne 1 then $ message,'Unknown dust model' polct=dustem_test_model_exists(model,/pol,/silent) use_polarisation=0 if keyword_set(polarisation) then use_polarisation=1 ;== INITIALISE DUSTEM dustem_init,model=model,pol=use_polarisation !dustem_verbose=1 !dustem_show_plot=1 !dustem_which='RELEASE' dir_in=!dustem_soft_dir ;== READ DEFAULT DUST MODEL VALUES st_model=dustem_read_all(dir_in) IF keyword_set(sed_outfile) THEN BEGIN use_sed_outfile=sed_outfile ENDIF ELSE BEGIN use_sed_outfile='./dustmodel_prediction_sed.xcat' ENDELSE IF keyword_set(ext_outfile) THEN BEGIN use_ext_outfile=ext_outfile ENDIF ELSE BEGIN use_ext_outfile='./dustmodel_prediction_ext.xcat' ENDELSE ;=== initialize filters for the SED Nfilt=n_elements(filters) sed=dustem_initialize_sed(Nfilt) sed.filter=filters sed.wave=dustem_filter2wav(filters) sed.instru=dustem_filter2instru(filters) ;=== initialize IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro sed[*].StokesI=1.e-10 sed.StokesQ=sed.StokesI/100. sed.StokesU=sed.StokesI/100. sed.SigmaII=sed.StokesI/100. sed.SigmaQQ=sed.StokesI/100. sed.SigmaUU=sed.StokesI/100. sed.SigmaIQ=sed.StokesI/100. sed.SigmaIU=sed.StokesI/100. sed.SigmaQU=sed.StokesI/100. ;=== initialize wavelengths for the EXTINCTION use_next=100 ; number of extinction measurements use_wave=10^(alog10(0.01)+alog10(50)*findgen(use_next)/float(use_next)) ; wavelengths at which extinction is defined if keyword_set(next) then use_next=next if keyword_set(wave) then $ use_wave=range_gen(use_next,[wave[0],wave[1]],/log) ext=dustem_initialize_ext(use_next) ext.instru='EXTINCTION' ext.filter='SPECTRUM' ext.wave=use_wave ;=== initialize IQU and associated errors to avoid problems when checking EXT in dustem_set_data.pro ext[*].EXT_I=1.e-10 ext.EXT_Q=ext.EXT_I/100. ext.EXT_U=ext.EXT_I/100. ext.SIGEXTII=ext.EXT_I/1000. ext.SIGEXTQQ=ext.EXT_I/1000. ext.SIGEXTUU=ext.EXT_I/1000. ext.SIGEXTIQ=ext.EXT_I/1000. ext.SIGEXTIU=ext.EXT_I/1000. ext.SIGEXTQU=ext.EXT_I/1000. dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext ;=== initialize at least one parameter as 'free' ;=== this is only to keep dustemwrap happy in its management of variables ;=== no fitting is going to be done pd = ['(*!dustem_params).gas.G0'] ; G0 pval = [1.] ;=== declare some other parameters as fixed ;=== this shows how to add any non-dust plugins to the observed spectra/SED ;=== for the dust-model parameters, it is not necessary, but allows the ;=== user to (i) change default values and (2) see the adopted values in the graphical output windows fpd=['(*!dustem_params).G0', $ '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;VSG mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;,$ 'dustem_plugin_synchrotron_1', $ ; 'dustem_plugin_continuum_1', $ ; 'dustem_plugin_continuum_2'] ; fpval = [1., $ st_model.grains(0).mdust_o_mh,$ st_model.grains(1).mdust_o_mh,$ st_model.grains(2).mdust_o_mh,$ ;,$ 2.,$ 800,$ 0.05] ;=== initialise all this information into dustemwrap's brain dustem_init_params,model,pd,pval,fpd=fpd,fiv=fpval,pol=use_polarisation ;=== compute the predicted SED and extinction curve for the dust-model ;=== and any plugins dustem_Ised=dustem_compute_sed(pval,st=dummy) dustem_Qsed=dustem_Ised*0. dustem_Used=dustem_Ised*0. ; uncomment following lines if running a polarisation model and you ;want predictions of QU ;toto=dustem_compute_stokes(pval,st=dummy) ;this procedure also allows for the extraction of the spectra ;dustem_qsed = toto[0] ;dustem_used = toto[1] dustem_iext=dustem_compute_ext(pval,st=dummy) dustem_Qext=dustem_Iext*0. dustem_Uext=dustem_Iext*0. ; uncomment following lines if running a polarisation model and you ;want predictions of QU ;toto=dustem_compute_stokext(pval,st=dummy) ;this procedure also allows for the extraction of the spectra ;dustem_qext = toto[0] ;dustem_uext = toto[1] ;== Now we have a predicted SED (I), we fill up a structure that we ;== will write to the sed_outfile sed.stokesI=dustem_Ised sed.stokesQ=dustem_Qsed sed.stokesU=dustem_Used ;=== set uncertainties for I,Q,U to 0. sed.sigmaII=0. sed.sigmaQQ=0. sed.sigmaUU=0. ;=== set covariances to 0. sed.sigmaIQ=0. sed.sigmaIU=0. sed.sigmaQU=0. ;==== fill in dependent columns of the SED. sed=dustem_fill_sed_dependent_columns(sed) ;======== save the predicted SED to a file write_xcat,sed,use_sed_outfile message,'Wrote '+use_sed_outfile,/continue ;== Now we have a predicted extinction (IQU), we fill up a structure that we ;== will write to the ext_outfile ext.EXT_I=dustem_Iext ext.EXT_Q=dustem_Qext ext.EXT_U=dustem_Uext ;=== set uncertainties to 0. ext.SIGEXTII=abs(ext.EXT_I*0.000001) ext.SIGEXTQQ=abs(ext.EXT_Q*0.000001) ext.SIGEXTUU=abs(ext.EXT_U*0.000001) ;=== set covariances to 0. ext.SIGEXTIQ=0. ext.SIGEXTIU=0. ext.SIGEXTQU=0. ;==== fill in dependent columns of the EXT. ext=dustem_fill_ext_dependent_columns(ext) ;======== save the EXT to a file write_xcat,ext,use_ext_outfile message,'Wrote '+use_ext_outfile,/continue ;=== MAKE A PLOT TO VISUALISE THE OUTPUT IF !dustem_noobj THEN BEGIN dustemwrap_plot_noobj,pval,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit ENDIF ELSE BEGIN dustemwrap_plot,pval,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit ENDELSE IF keyword_set(fits_save) THEN BEGIN message,'Writing out results structure: '+fits_save,/info dustem_write_fits_table,filename=fits_save,help=help ENDIF the_end: message,'Finished dustem_run_prediction_example',/info if keyword_set(noexit) then stop END