OV-GSO-DC / dustem-wrapper_idl

Commit 5d34557b6b55768a68d50ef7b9885ee7e64768eb

Authored by Jean-Philippe Bernard
1 parent 21b2be95
Exists in master

First commit

Showing 2 changed files with 200 additions and 0 deletions   Show diff stats
Inline Side-by-side
Data/EXAMPLE_OBSDATA/example_sed3.xcat 0 → 100644
  Diff comments   View file @5d34557
... ... @@ -0,0 +1,20 @@
  1 +\Written by write_xcat.pro on Thu Sep 15 14:18:31 2022
  2 +|INSTRU |FILTER |WAVE |STOKESI |STOKESQ |STOKESU |LARGEP |SMALLP |PSI |SIGMAII |SIGMAQQ |SIGMAUU |SIGMAIQ |SIGMAIU |SIGMAQU |SIGMA_LARGEP |SIGMA_SMALLP |SIGMA_PSI |
  3 +|char |char |real |real |real |real |real |real |real |real |real |real |real |real |real |real |real |real |
  4 +|NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |NULL |
  5 +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
  6 +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
  7 +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
  8 +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
  9 +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
  10 +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
  11 +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
  12 +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
  13 +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
  14 +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
  15 +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
  16 +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
  17 +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
  18 +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
  19 +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
  20 +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
... ...
src/idl/dustem_fit_polarisation_example.pro 0 → 100644
  Diff comments   View file @5d34557
... ... @@ -0,0 +1,180 @@
  1 +PRO dustem_fit_polarisation_example,model=model $
  2 + ,sed_file=sed_file $
  3 + ,Nitermax=Nitermax $
  4 + ,postscript=postscript $
  5 + ,fits_save_and_restore=fits_save_and_restore $
  6 + ,help=help
  7 +
  8 +;This routine is an example of how to fit observational SEDs with
  9 +;DustEM and DustEMWrap. The objective is to illustrate how to use DustEMWrap
  10 +; (and not to do science -- the fit obtained by running this example is
  11 +; likely to be poor!)
  12 +;
  13 +;For this example, the code uses the SED in the file example_SED_1.xcat,
  14 +;which is distributed in the Data/EXAMPLE_OBSDATA/ directory
  15 +;
  16 +;The example SED has Stokes I photometric data points from SPITZER
  17 +;IRAC and MIPS and IRAS. Examples illustrating running DustEMWrap to
  18 +;fit data with spectral data, polarisation data and extinction data
  19 +;are provided in other _readme routines in the src/idl/ directory.
  20 +
  21 +;; ***COMMENT AH***
  22 +;; Remove or update the next few lines -->>
  23 +;; No spectrum data points.
  24 +;;SPECTRUM data points can be included and the corresponding filter
  25 +;;filed must read SPECTRUM. Note that its is note necessary
  26 +;;to use the .xcat file format, and data SED can be provided
  27 +;;manually, but the observation structure must have the structure shown below.
  28 +;;For the example to work, the DustEM and DustEMWrap packages must have
  29 +;;been configured and installed succesfully.
  30 +;;See dustem_cvs_readme.txt for install instructions).
  31 +;; *** END COMMENT AH***
  32 +
  33 +;+
  34 +; NAME:
  35 +; dustem_fit_polarisation_example
  36 +; PURPOSE:
  37 +; This is an example of how to fit SEDs with DustEMWrap.
  38 +; It is intended as an example to follow when writing your own
  39 +; programs to analyse data with DustEMWrap.
  40 +; CATEGORY:
  41 +; DustEMWrap, Distributed, High-Level, User Example
  42 +; CALLING SEQUENCE:
  43 +; dustem_fit_polarisation_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save_and_restore=][,/help]
  44 +; INPUTS:
  45 +; None
  46 +; OPTIONAL INPUT PARAMETERS:
  47 +; None
  48 +; OUTPUTS:
  49 +; None
  50 +; OPTIONAL OUTPUT PARAMETERS:
  51 +; Plots, Results save structure
  52 +; ACCEPTED KEY-WORDS:
  53 +; model = specifies the interstellar dust mixture used by DustEM
  54 +; 'MC10' model from Compiegne et al 2010 (default)
  55 +; 'DBP90' model from Desert et al 1990
  56 +; 'DL01' model from Draine & Li 2001
  57 +; 'WD01_RV5p5B' model from Weingartner & Draine 2002 with Rv=5.5
  58 +; 'DL07' model from Draine & Li 2007
  59 +; 'J13' model from Jones et al 2013, as updated in
  60 +; Koehler et al 2014
  61 +; 'G17_ModelA' model A from Guillet et al (2018). Includes
  62 +; polarisation. See Tables 2 and 3 of that paper for details.
  63 +; 'G17_ModelB' model B from Guillet et al (2018)
  64 +; 'G17_ModelC' model C from Guillet et al (2018)
  65 +; 'G17_ModelD' model A from Guillet et al (2018)
  66 +; sed_file = string naming the path to text file in .xcat format that
  67 +; describes the observational SED. If not set, the file
  68 +; 'Data/SEDs/sample_SED.xcat' is used.
  69 +; postscript = if set, final plot is saved as postscript in the
  70 +; current working directory
  71 +; Nitermax = maximum number of fit iterations. Default is 5.
  72 +; fits_save_and_restore = if set, saves results in a fits file, restore the file and plot restored results
  73 +; help = if set, print this help
  74 +; COMMON BLOCKS:
  75 +; None
  76 +; SIDE EFFECTS:
  77 +; None
  78 +; RESTRICTIONS:
  79 +; The DustEM fortran code must be installed
  80 +; The DustEMWrap IDL code must be installed
  81 +; PROCEDURES AND SUBROUTINES USED:
  82 +; *** COMMENT AH --> is this really NONE? ****
  83 +; EXAMPLES
  84 +; dustem_fit_polarisation_example
  85 +; dustem_fit_polarisation_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits'
  86 +; dustem_fit_polarisation_example,model='DBP90'
  87 +; MODIFICATION HISTORY:
  88 +; Written by JPB Apr-2022
  89 +; Evolution details on the DustEMWrap gitlab.
  90 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
  91 +;-
  92 +
  93 +IF keyword_set(help) THEN BEGIN
  94 + doc_library,'dustem_fit_polarisation_example'
  95 + goto,the_end
  96 +END
  97 +
  98 +IF keyword_set(model) THEN BEGIN
  99 + use_model=strupcase(model)
  100 +ENDIF ELSE BEGIN
  101 + use_model='G17_MODELC' ;Default
  102 +ENDELSE
  103 +
  104 +dustem_define_la_common
  105 +
  106 +;=== initialise dustem
  107 +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;
  108 +
  109 +!dustem_verbose=1
  110 +!dustem_show_plot=1
  111 +!dustem_nocatch=1
  112 +
  113 +filename=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/example_sed3.xcat'
  114 +sed=read_xcat(filename)
  115 +
  116 +dustem_set_data,sed
  117 +
  118 +;=== Set which parameters you want to fit
  119 +pd = [ $
  120 + '(*!dustem_params).G0', $ ;G0
  121 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  122 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction
  123 + '(*!dustem_params).grains(2).mdust_o_mh' $ ;PAH1 mass fraction
  124 + ;'(*!dustem_params).grains(3).mdust_o_mh',$ ;amCBEx
  125 + ;'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude
  126 + ;'dustem_plugin_synchrotron_4', $ ;Synchrotron polarization angle
  127 + ;'dustem_plugin_modify_dust_polarization_2', $ ;Polarization angle applied to the Fortran dust model via a core plugin
  128 + ;'dustem_plugin_modify_spinning_polarization_1', $
  129 + ;'dustem_plugin_modify_spinning_polarization_2' $
  130 + ]
  131 +
  132 +;p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,1.6500E-04,68.,4e-3,20]
  133 +;iv=p_truth+[0.,8.e-4,8.e-4,8.e-4,8.e-4,-10.,2e-4,10] ;shifted from solution
  134 +
  135 +;p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,68.,1.2,20]
  136 +p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04]
  137 +iv=p_truth+[0.,1.e-4,1.e-4,1.e-4] ;shifted from solution
  138 +
  139 +Npar=n_elements(pd)
  140 +ulimed=replicate(0,Npar)
  141 +llimed=replicate(1,Npar)
  142 +llims=replicate(0,Npar)
  143 +
  144 +;stop
  145 +
  146 +dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims;,tied=tied
  147 +;stop
  148 +dustem_init_plugins,pd
  149 +
  150 +;=== SET THE OBSERVATION STRUCTURE
  151 +
  152 +;ind=where(sed.wave GE 0.1)
  153 +
  154 +;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
  155 +
  156 +;=== RUN fit
  157 +tol=1.e-16
  158 +xtol=1.e-16
  159 +Nitermax=30; just an initialization
  160 +;for i=0L,n_tags(!dustem_data)-1 do begin
  161 +;if ptr_valid(!dustem_data.(i)) then Nitermax+=1
  162 +;endfor
  163 +;Nitermax=4 ;maximum number of iteration. This is the criterion which will stop the fit procedure
  164 +
  165 +xrange=[1.,2.e4]
  166 +yrange=[1e-4,1.e3]
  167 +
  168 +loadct,13
  169 +;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine
  170 +t1=systime(0,/sec)
  171 +title='CACA' ;does not work if title is undefined ....
  172 +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)
  173 +t2=systime(0,/sec)
  174 +
  175 +stop
  176 +
  177 +the_end:
  178 +
  179 +END
  180 +
... ...