OV-GSO-DC / dustem-wrapper_idl

Commit d0bfcf1e8d46b2f9cb4be5d06c94ab538c424fe8

Authored by Annie Hughes
1 parent 909b77f6
Exists in master

moved dustem_fit_J13.pro to ONHOLD dir outside git repo

Showing 1 changed file with 0 additions and 266 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_fit_J13.pro deleted
View file @909b77f
... ... @@ -1,266 +0,0 @@
1   -PRO dustem_fit_j13,help=help,noplot=noplot,mode=mode
2   -
3   -;+
4   -; NAME:
5   -; dustem_fit
6   -; PURPOSE:
7   -; Runs conjointly a fit on extinction, emission and polarization by extinction
8   -; CATEGORY:
9   -; Dustem
10   -; CALLING SEQUENCE:
11   -; sed=dustem_fit([/noplot,/help])
12   -; INPUTS:
13   -; None
14   -; OPTIONAL INPUT PARAMETERS:
15   -; noplot : do not display figures
16   -; OUTPUTS:
17   -; None
18   -; OPTIONAL OUTPUT PARAMETERS:
19   -; ACCEPTED KEY-WORDS:
20   -; help = If set, print this help
21   -; COMMON BLOCKS:
22   -; None
23   -; SIDE EFFECTS:
24   -; None
25   -; RESTRICTIONS:
26   -; The dustem idl wrapper must be installed
27   -; PROCEDURE:
28   -; None
29   -; EXAMPLES
30   -;
31   -; MODIFICATION HISTORY:
32   -; TLS added by D. Paradis.
33   -; Written by V. Guillet (2012)
34   -; see evolution details on the dustem cvs maintained at CESR
35   -; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
36   -;-
37   -
38   -IF keyword_set(mode) THEN BEGIN
39   - use_mode=strupcase(mode)
40   -ENDIF ELSE BEGIN
41   - use_mode='COMPIEGNE_ETAL2010' ;Default is last dustem model
42   -ENDELSE
43   -
44   -; Verbose
45   -defsysv, '!dustem_verbose', 1
46   -
47   -data_files = { $
48   - ; sed: !dustem_wrap_soft_dir + 'Data/SEDs/Gal_composite_spectrum.xcat', $
49   - sed: '/Users/paradis/Soft_Librairies/Planck/MW/HA_diffuse_COlevel.xcat' $
50   - ; sed: !dustem_wrap_soft_dir + 'Data/SEDs/spec_firas_chi2.xcat', $
51   - ; ext: !dustem_wrap_soft_dir + 'Data/EXTs/spec_sk-69228-2.xcat' $
52   - ; polext: !dustem_wrap_soft_dir + 'Data/POLEXTs/Serkowski_DISM.xcat', $
53   - ; polsed: '' $
54   - }
55   -
56   -; Local variable, to define global !fit_rchi2_weight
57   -rchi2_weight = { $
58   - sed: 1. $
59   - ; ext: 1. $
60   - ; polext: 1. $
61   - ; polsed: 1. $
62   - }
63   -
64   -; Flags to determine what to fit
65   -defsysv, '!dustem_data', { $ ;Data to fit
66   - sed: ptr_new() $
67   - ; ext: ptr_new() $
68   - ; polext: ptr_new() $
69   -; polsed: ptr_new() $
70   - }
71   -
72   -; Define global !fit_rchi2_weight with the same structure and field order like !dustem_data
73   -tagnames = tag_names(!dustem_data)
74   -instr="defsysv, '!fit_rchi2_weight', {"
75   -for i = 0, n_elements(tagnames)-1 do instr+=tagnames(i)+': rchi2_weight.'+tagnames(i)+', '
76   -instr+='empty: 0. }
77   -toto=execute(instr)
78   -
79   -;=== Set weight for reduced chi2 ponderation between emission, extinction, polarized extinction and emission ===
80   -; Calculate total according to fileds in !dustem_data structure
81   -total = 0
82   -for i = 0, n_tags(!dustem_data)-1 do total = total + !fit_rchi2_weight.(i)
83   -; values are relative to total
84   -for i = 0, n_tags(!dustem_data)-1 do !fit_rchi2_weight.(i) /= total
85   -
86   -dustem_init,mode=use_mode
87   -
88   -if keyword_set(noplot) then defsysv,'!dustem_show_plot', 0
89   -
90   -; Check that polarisation has been run (-pol in GRAIN.DAT) if fit on polarization is asked for
91   -if (tag_exist(!dustem_data,'polext') or tag_exist(!dustem_data,'polsed')) and !run_pol eq 0. then begin
92   - message,"pol option must be set in GRAIN.DAT to fit polarization",/info
93   - stop
94   -endif
95   -
96   -; Check that DTLS has been run (-dtls in GRAIN.DAT)
97   -if (tag_exist(!dustem_data,'dtls') ) and !run_tls eq 0. then begin
98   - message,"dtls option must be set in GRAIN.DAT to use the TLS model",/info
99   - stop
100   -endif
101   -
102   -;=== Read observational data ===
103   -for i = 0, n_elements(tagnames)-1 do begin
104   - instr='!dustem_data.('+strtrim(i,2)+') = dustem_set_data(read_xcat(data_files.'+tagnames(i)+',/silent))'
105   - toto=execute(instr)
106   -endfor
107   -
108   -;=== Parameters for fit ===
109   -pd_mass = [ $
110   - '(*!dustem_params).grains(0).mdust_o_mh', $ ;PAH0 mass fraction
111   - '(*!dustem_params).grains(1).mdust_o_mh', $ ;PAH1 mass fraction
112   - '(*!dustem_params).grains(2).mdust_o_mh' $ ;Small amCBEx
113   - ; '(*!dustem_params).grains(3).mdust_o_mh' $ ;Large amCBEx
114   - ]
115   -
116   -;tls=['(*!dustem_params).tls.a_dtls', $
117   -; '(*!dustem_params).tls.lc' $
118   -; ; '(*!dustem_params).tls.vm' $
119   -;]
120   -
121   -iv_tls=[10., 9.]
122   -
123   -iv_mass = [ $
124   - 2.23e-4, $
125   - 2.64e-6, $
126   - 3.5e-4 $
127   -; 5.7e-4 $
128   - ] ; DUST MASSES
129   -;iv_mass = [ $
130   -; 9.3E-4, $
131   -; 0.000449, $
132   -; 2.2e-3 $
133   -;; 0.0008 $
134   -;; 7.8e-4 $ ;7.8e-3
135   -; ] ; DUST MASSES
136   -;iv_mass = [ $
137   -; 1.1e-3, $
138   -; 1.8e-9, $
139   -; 6.e-3, $
140   -; 0.01 $
141   -; ] ; DUST MASSES
142   -
143   -pd_sdist = [ $
144   -;; '(*!dustem_params).grains(0).amin', $ ;PAH0 minmal size
145   -;; '(*!dustem_params).grains(0).amax' $ ;PAH0 maximal size
146   -; '(*!dustem_params).grains(0).alpha_o_a0', $ ;PAH0 mean size
147   - '(*!dustem_params).grains(0).alpha_o_a0' $ ;PAH0 mean size
148   -
149   -;; '(*!dustem_params).grains(0).at', $ ;PAH0 sigma size
150   -; '(*!dustem_params).grains(0).amin', $ ;PAH1 minmal size
151   -; '(*!dustem_params).grains(0).amax' $ ;PAH1 maximal size
152   -; '(*!dustem_params).grains(2).alpha_o_a0' $ ;PAH1 mean size
153   -;; '(*!dustem_params).grains(1).at', $ ;PAH1 sigma size
154   -;; '(*!dustem_params).grains(2).amin', $ ;VSG minimal size
155   -;; '(*!dustem_params).grains(2).amax', $ ;VSG maximal size
156   -; '(*!dustem_params).grains(3).amin', $ ;Carbon BG minimal size
157   -;; '(*!dustem_params).grains(3).amax', $ ;Carbon BG maximal size
158   -; '(*!dustem_params).grains(2).alpha_o_a0', $ ;Carbon BG power law slope
159   -; '(*!dustem_params).grains(3).alpha_o_a0' $ ;Carbon BG power law slope
160   -
161   -; '(*!dustem_params).grains(3).at', $ ;Carbon BG threshold size (exp. decay)
162   -; '(*!dustem_params).grains(3).ac', $ ;Carbon BG cut-off size (exp. decay)
163   -; '(*!dustem_params).grains(3).gamma', $ ;Carbon BG exp. decay slope
164   -; '(*!dustem_params).grains(4).amin', $ ;Silicate BG minimal size
165   -;; '(*!dustem_params).grains(4).amax', $ ;Silicate BG maximal size
166   -; '(*!dustem_params).grains(4).alpha_o_a0', $ ;Silicate BG power law slope
167   -; '(*!dustem_params).grains(4).at', $ ;Silicate BG threshold size (exp. decay)
168   -; '(*!dustem_params).grains(4).ac', $ ;Silicate BG cut-off size (exp. decay)
169   -; '(*!dustem_params).grains(4).gamma' $ ;Silicate BG exp. decay slope
170   - ]
171   -iv_sdist = [ $
172   -;3.5000E-08, 1.2000E-07, 6.4000E-08, 1.0000E-01 $ ; PAH0
173   -;;3.5000E-08, 1.2000E-07, 6.4000E-08, 1.0000E-01, $ ; PAH1
174   -;-4.98, 3.12e-7 $
175   --4. $
176   -;1.7e-8,-3.21 $ ; VSG
177   -;;6.0000E-07, 5.0000E-04,-2.6000E+00, 0.7000E-05, 0.7000E-05, 0.9000E+00, $ Carbon BG
178   -;;6.0000E-07, 5.0000E-04,-3.4000E+00, 3.0000E-05, 1.0000E-05, 0.4200E+00 $ Silicate BG
179   -;6.0000E-07, -2.6000E+00, 0.7000E-05, 0.7000E-05, 0.9000E+00, $ Carbon BG
180   -;6.0000E-07, -3.4000E+00, 3.0000E-05, 1.0000E-05, 0.4200E+0 $ Silicate BG
181   - ]
182   -
183   -pd_pol_alignment = [ $
184   - '(*!dustem_params).pol(4).atresh', $ ; Threshold size for grain alignment
185   - '(*!dustem_params).pol(4).pstiff', $ ; Stiffness of transition between non-aligned and aligned grains sizes
186   - '(*!dustem_params).pol(4).plev' $ ; Maximal alignment efficiency of grains
187   - ]
188   -iv_pol_alignment = [ $
189   - 0.1, $
190   - 0.5, $
191   - 1.0 $
192   - ]
193   -
194   -pd_G0=['(*!dustem_params).G0'] ;
195   -iv_G0=[4.]
196   -
197   -pd_cont = ['dustem_create_continuum_2']
198   -iv_cont = [0.0039]
199   -;stop
200   -;pd = [ pd_pol_alignment, pd_mass ]
201   -;iv = [ iv_pol_alignment, iv_mass ]
202   -
203   -; Parameters
204   -;pd = [ pd_G0, pd_mass, pd_pol_alignment, pd_sdist, pd_cont ]
205   -;iv = [ iv_G0, iv_mass, iv_pol_alignment, iv_sdist, iv_cont ]
206   -;pd = [ pd_mass, pd_sdist, pd_cont ]
207   -;iv = [ iv_mass, iv_sdist, iv_cont ]
208   -;pd = [ pd_mass, tls,pd_cont ]
209   -;iv = [ iv_mass, iv_tls,iv_cont ]
210   -pd = [pd_G0, pd_mass, pd_cont,pd_sdist];,tls]
211   -iv = [iv_G0, iv_mass, iv_cont,iv_sdist];,iv_tls]
212   -;pd = [ pd_cont,tls,pd_G0 ]
213   -;iv = [ iv_cont,iv_tls,iv_G0 ]
214   -
215   -Npar=n_elements(pd)
216   -ulimed=replicate(0,Npar)
217   -llimed=replicate(1,Npar)
218   -llims=replicate(0.,Npar)
219   -ulims=replicate(0.,Npar)
220   -
221   -; Upper limit on alignment efficiency
222   -;ulimed(1+n_elements(pd_mass)+2) = 1
223   -;ulims(1+n_elements(pd_mass)+2) = 1.
224   -; Upper limit on grain size
225   -; amCBEx
226   -;ulimed(1+n_elements(pd_mass)+n_elements(pd_pol_alignment)+9) = 1
227   -;ulims(1+n_elements(pd_mass)+n_elements(pd_pol_alignment)+9) = 9.9d-4
228   -
229   -;Npar=n_elements(pd)
230   -;ulimed=replicate(0,Npar)
231   -;llimed=replicate(1,Npar)
232   -;llims=replicate(0.,Npar)
233   -;ulims=replicate(0.,Npar)
234   -;stop
235   -; ulimed=[0 , 0. ,0 ,0, 0, 0,0,1] ; aucun des parametres n'a une valeur superieure limite
236   -;ulims=[0 , 0. ,0 ,0, 0, 0, 0,5.e-6]
237   - llimed=[1, 1, 1., 1.,1,1.] ; tous les parametres ont une valeur inferieure limite
238   - llims =[0.1,0, 0,0.,0.,-2];, 1.e-8];,0,0,3.]
239   -
240   -
241   -;=== Fixed parameters
242   -;fpd=['(*!dustem_params).G0'] ;Fixed parameter description
243   -;fiv=[1000]
244   -;dustem_init_fixed_params,fpd,fiv
245   -fpd=['(*!dustem_params).grains(3).mdust_o_mh'];,'(*!dustem_params).grains(3).mdust_o_mh'] ;Fixed parameter description
246   -fiv=[1.e-16];,1.e-16]
247   -dustem_init_fixed_params,fpd,fiv
248   -
249   -;== SET THE FITTED PARAMETERS
250   -dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
251   -
252   -;=== RUN fit
253   -tol=1.e-13
254   -Nitermax=20 ;maximum number of iteration. This is the criterium which will stop the fit procedure
255   -
256   -loadct,13
257   -!y.range=[1.e-3,100.] ;This is to ajust plot range from outside the routine
258   -t1=systime(0,/sec)
259   -res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax) ;,func_name=func_name,cf_min=cf_min)
260   -t2=systime(0,/sec)
261   -print,'executed in ',t2-t1,' sec'
262   -
263   -;st=dustem_read_all_res(!dustem_res,/silent)
264   -
265   -END
266   -