From 5d34557b6b55768a68d50ef7b9885ee7e64768eb Mon Sep 17 00:00:00 2001 From: Jean-Philippe Bernard Date: Thu, 15 Sep 2022 14:40:18 +0200 Subject: [PATCH] First commit --- Data/EXAMPLE_OBSDATA/example_sed3.xcat | 20 ++++++++++++++++++++ src/idl/dustem_fit_polarisation_example.pro | 180 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 200 insertions(+), 0 deletions(-) create mode 100644 Data/EXAMPLE_OBSDATA/example_sed3.xcat create mode 100644 src/idl/dustem_fit_polarisation_example.pro diff --git a/Data/EXAMPLE_OBSDATA/example_sed3.xcat b/Data/EXAMPLE_OBSDATA/example_sed3.xcat new file mode 100644 index 0000000..d9bcede --- /dev/null +++ b/Data/EXAMPLE_OBSDATA/example_sed3.xcat @@ -0,0 +1,20 @@ +\Written by write_xcat.pro on Thu Sep 15 14:18:31 2022 +|INSTRU |FILTER |WAVE |STOKESI |STOKESQ |STOKESU |LARGEP |SMALLP |PSI |SIGMAII |SIGMAQQ |SIGMAUU |SIGMAIQ |SIGMAIU |SIGMAQU |SIGMA_LARGEP |SIGMA_SMALLP |SIGMA_PSI | +|char |char |real |real |real |real |real |real |real |real |real |real |real |real |real |real |real |real | +|NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL | +IRAS IRAS1 12.0000 0.0312461 0.000293617 0.000106868 -32768.0 -32768.0 -32768.0 3.12461e-08 2.93617e-10 1.06868e-10 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +IRAS IRAS2 25.0000 0.0332762 0.000312694 0.000113811 -32768.0 -32768.0 -32768.0 3.32762e-08 3.12694e-10 1.13811e-10 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +IRAS IRAS3 60.0000 0.121636 0.00114300 0.000416018 -32768.0 -32768.0 -32768.0 1.21636e-07 1.14300e-09 4.16018e-10 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +IRAS IRAS4 100.000 0.363825 0.00341883 0.00124435 -32768.0 -32768.0 -32768.0 3.63825e-07 3.41883e-09 1.24435e-09 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +PACS PACS3 160.000 0.546915 0.00513932 0.00187056 -32768.0 -32768.0 -32768.0 5.46915e-07 5.13932e-09 1.87056e-09 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +SPIRE SPIRE1 250.000 0.424589 0.00398983 0.00145218 -32768.0 -32768.0 -32768.0 4.24589e-07 3.98983e-09 1.45218e-09 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +SPIRE SPIRE2 350.000 0.244398 0.00229659 0.000835890 -32768.0 -32768.0 -32768.0 2.44398e-07 2.29659e-09 8.35890e-10 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +SPIRE SPIRE3 500.000 0.112675 0.00105880 0.000385371 -32768.0 -32768.0 -32768.0 1.12675e-07 1.05880e-09 3.85371e-10 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +HFI HFI2 550.078 0.0883461 0.000830182 0.000302162 -32768.0 -32768.0 -32768.0 8.83461e-08 8.30182e-10 3.02162e-10 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +HFI HFI3 849.270 0.0271650 0.000255268 9.29098e-05 -32768.0 -32768.0 -32768.0 2.71650e-08 2.55268e-10 9.29098e-11 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +HFI HFI4 1381.53 0.00613173 5.76194e-05 2.09717e-05 -32768.0 -32768.0 -32768.0 6.13173e-09 5.76194e-11 2.09717e-11 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +HFI HFI5 2096.45 0.00149945 1.40903e-05 5.12844e-06 -32768.0 -32768.0 -32768.0 1.49945e-09 1.40903e-11 5.12844e-12 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +HFI HFI6 2997.92 0.000522665 4.91144e-06 1.78762e-06 -32768.0 -32768.0 -32768.0 5.22665e-10 4.91144e-12 1.78762e-12 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +LFI LFI1 4285.71 0.000171805 1.61443e-06 5.87606e-07 -32768.0 -32768.0 -32768.0 1.71805e-10 1.61443e-12 5.87606e-13 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +LFI LFI2 6818.18 8.11966e-05 7.62998e-07 2.77709e-07 -32768.0 -32768.0 -32768.0 8.11966e-11 7.62998e-13 2.77709e-13 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 +LFI LFI3 10000.0 0.000125388 1.17826e-06 4.28851e-07 -32768.0 -32768.0 -32768.0 1.25388e-10 1.17826e-12 4.28851e-13 0.00000 0.00000 0.00000 -32768.0 -32768.0 -32768.0 diff --git a/src/idl/dustem_fit_polarisation_example.pro b/src/idl/dustem_fit_polarisation_example.pro new file mode 100644 index 0000000..e63c78f --- /dev/null +++ b/src/idl/dustem_fit_polarisation_example.pro @@ -0,0 +1,180 @@ +PRO dustem_fit_polarisation_example,model=model $ + ,sed_file=sed_file $ + ,Nitermax=Nitermax $ + ,postscript=postscript $ + ,fits_save_and_restore=fits_save_and_restore $ + ,help=help + +;This routine is an example of how to fit observational SEDs with +;DustEM and DustEMWrap. The objective is to illustrate how to use DustEMWrap +; (and not to do science -- the fit obtained by running this example is +; likely to be poor!) +; +;For this example, the code uses the SED in the file example_SED_1.xcat, +;which is distributed in the Data/EXAMPLE_OBSDATA/ directory +; +;The example SED has Stokes I photometric data points from SPITZER +;IRAC and MIPS and IRAS. Examples illustrating running DustEMWrap to +;fit data with spectral data, polarisation data and extinction data +;are provided in other _readme routines in the src/idl/ directory. + +;; ***COMMENT AH*** +;; Remove or update the next few lines -->> +;; No spectrum data points. +;;SPECTRUM data points can be included and the corresponding filter +;;filed must read SPECTRUM. Note that its is note necessary +;;to use the .xcat file format, and data SED can be provided +;;manually, but the observation structure must have the structure shown below. +;;For the example to work, the DustEM and DustEMWrap packages must have +;;been configured and installed succesfully. +;;See dustem_cvs_readme.txt for install instructions). +;; *** END COMMENT AH*** + +;+ +; NAME: +; dustem_fit_polarisation_example +; PURPOSE: +; This is an example of how to fit SEDs with DustEMWrap. +; It is intended as an example to follow when writing your own +; programs to analyse data with DustEMWrap. +; CATEGORY: +; DustEMWrap, Distributed, High-Level, User Example +; CALLING SEQUENCE: +; dustem_fit_polarisation_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save_and_restore=][,/help] +; INPUTS: +; None +; OPTIONAL INPUT PARAMETERS: +; None +; OUTPUTS: +; None +; OPTIONAL OUTPUT PARAMETERS: +; Plots, Results save structure +; ACCEPTED KEY-WORDS: +; model = specifies the interstellar dust mixture used by DustEM +; 'MC10' model from Compiegne et al 2010 (default) +; 'DBP90' model from Desert et al 1990 +; 'DL01' model from Draine & Li 2001 +; 'WD01_RV5p5B' model from Weingartner & Draine 2002 with Rv=5.5 +; 'DL07' model from Draine & Li 2007 +; 'J13' model from Jones et al 2013, as updated in +; Koehler et al 2014 +; 'G17_ModelA' model A from Guillet et al (2018). Includes +; polarisation. See Tables 2 and 3 of that paper for details. +; 'G17_ModelB' model B from Guillet et al (2018) +; 'G17_ModelC' model C from Guillet et al (2018) +; 'G17_ModelD' model A from Guillet et al (2018) +; sed_file = string naming the path to text file in .xcat format that +; describes the observational SED. If not set, the file +; 'Data/SEDs/sample_SED.xcat' is used. +; postscript = if set, final plot is saved as postscript in the +; current working directory +; Nitermax = maximum number of fit iterations. Default is 5. +; fits_save_and_restore = if set, saves results in a fits file, restore the file and plot restored results +; help = if set, print this help +; COMMON BLOCKS: +; None +; SIDE EFFECTS: +; None +; RESTRICTIONS: +; The DustEM fortran code must be installed +; The DustEMWrap IDL code must be installed +; PROCEDURES AND SUBROUTINES USED: +; *** COMMENT AH --> is this really NONE? **** +; EXAMPLES +; dustem_fit_polarisation_example +; dustem_fit_polarisation_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits' +; dustem_fit_polarisation_example,model='DBP90' +; MODIFICATION HISTORY: +; Written by JPB Apr-2022 +; 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_fit_polarisation_example' + goto,the_end +END + +IF keyword_set(model) THEN BEGIN + use_model=strupcase(model) +ENDIF ELSE BEGIN + use_model='G17_MODELC' ;Default +ENDELSE + +dustem_define_la_common + +;=== initialise dustem +dustem_init,mode=use_model,kwords=['logn-chrg-spin','plaw-pol','plaw-pol'],/pol;,kwords=['plaw-ed','logn','logn'];,/pol;kwords=['plaw-ed-chrg-spin'],/pol;,kwords=['plaw-ed'];kwords=['plaw-ed-chrg-spin'];,/pol; + +!dustem_verbose=1 +!dustem_show_plot=1 +!dustem_nocatch=1 + +filename=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/example_sed3.xcat' +sed=read_xcat(filename) + +dustem_set_data,sed + +;=== Set which parameters you want to fit +pd = [ $ + '(*!dustem_params).G0', $ ;G0 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction + '(*!dustem_params).grains(2).mdust_o_mh' $ ;PAH1 mass fraction + ;'(*!dustem_params).grains(3).mdust_o_mh',$ ;amCBEx + ;'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude + ;'dustem_plugin_synchrotron_4', $ ;Synchrotron polarization angle + ;'dustem_plugin_modify_dust_polarization_2', $ ;Polarization angle applied to the Fortran dust model via a core plugin + ;'dustem_plugin_modify_spinning_polarization_1', $ + ;'dustem_plugin_modify_spinning_polarization_2' $ + ] + +;p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,1.6500E-04,68.,4e-3,20] +;iv=p_truth+[0.,8.e-4,8.e-4,8.e-4,8.e-4,-10.,2e-4,10] ;shifted from solution + +;p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,68.,1.2,20] +p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04] +iv=p_truth+[0.,1.e-4,1.e-4,1.e-4] ;shifted from solution + +Npar=n_elements(pd) +ulimed=replicate(0,Npar) +llimed=replicate(1,Npar) +llims=replicate(0,Npar) + +;stop + +dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims;,tied=tied +;stop +dustem_init_plugins,pd + +;=== SET THE OBSERVATION STRUCTURE + +;ind=where(sed.wave GE 0.1) + +;st_bidon=dustem_set_data(m_fit=sed[ind],rchi2_weight=rchi2_weight,f_HI=f_HI) ; comment if second line is to be used + +;=== RUN fit +tol=1.e-16 +xtol=1.e-16 +Nitermax=30; just an initialization +;for i=0L,n_tags(!dustem_data)-1 do begin +;if ptr_valid(!dustem_data.(i)) then Nitermax+=1 +;endfor +;Nitermax=4 ;maximum number of iteration. This is the criterion which will stop the fit procedure + +xrange=[1.,2.e4] +yrange=[1e-4,1.e3] + +loadct,13 +;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine +t1=systime(0,/sec) +title='CACA' ;does not work if title is undefined .... +res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []',title=title) +t2=systime(0,/sec) + +stop + +the_end: + +END + -- libgit2 0.21.2