OV-GSO-DC / dustem-wrapper_idl

Commit 607060e51d9cbaf22569d6010e233afcc9a4879f

Authored by Ilyes Choubani
1 parent ba82c5cf
Exists in master

test version

Showing 18 changed files with 253 additions and 138 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_activate_plugins.pro
  Diff comments   View file @607060e
... ... @@ -41,6 +41,8 @@ ENDIF
41 41  
42 42 p_dim = p_min * (*(*!dustem_fit).param_init_values)
43 43  
  44 +st=dustem_run(p_dim)
  45 +
44 46  
45 47 f=1 ; Initializing the index that is associated to each type of plugin
46 48 FOR i=0L,n_elements(p_min)-1 DO BEGIN
... ... @@ -80,24 +82,33 @@ FOR i=0L,n_elements(p_min)-1 DO BEGIN
80 82  
81 83 ;==============Assigning data to the initialized plugin data arrays==============
82 84  
  85 +
  86 + print, 'ftn is:'
  87 + print, ftn
  88 + print, 'key is:'
  89 + print, index
  90 + print, 'val is:'
  91 + print, value
  92 +
  93 +
83 94 if strmid(ftn,0,13) eq 'dustem_plugin' then begin
84 95  
85 96  
86 97 k=where(strmid(tag_names(*!dustem_scope),0,8) eq strmid(strupcase(strmid(ftn,14)),0,8),counte) ; Selecting a plugin through matching the string name of the plugin form the scope system variable with the one read from the parameter description vector
87 98  
88   - ;Dry run of the plugins to obtain their scopes and run them accordingly (an advanced user might want to add their own lines here)
89   - str='toto='+ftn+'(scope=scope)'
  99 + ;Dry run of the plugins to obtain their scopes and run them accordingly (an advanced user might want to add their own lines here)
  100 + str='toto='+ftn+'(scope=scope)' & str=str(0)
90 101 toto=execute(str)
91 102  
92 103 (*!dustem_scope).(k)=ptr_new(scope)
93   -
94   -
  104 +
  105 + tgnms=tag_names((*!dustem_scope))
  106 +
95 107 str='((*!dustem_plugin).('+strtrim(k,2)+'))=ptr_new('+ftn+'(key=index,val=value)'+')' & str=str(0)
96 108 toto=execute(str) & IF !dustem_verbose NE 0 THEN message,strupcase(strmid(ftn,7)),/info
97   -
98   - ENDIF
99   -
100   - f=f+1 & i=i+count-1 ; Incrementing the parameter and same-type plugin indices
  109 +
  110 + ENDIF
  111 + f=f+1 & i=i+count-1 ; Incrementing the parameter and same-type plugin indices
101 112 ENDIF
102 113 ENDFOR
103 114 the_end:
... ...
src/idl/dustem_compute_stokes.pro
  Diff comments   View file @607060e
1   -FUNCTION dustem_compute_stokes,p_dim,stp,dustem_qsed=dustem_qsed,dustem_used=dustem_used,Q_spec=Q_spec,U_spec=U_spec,out_st=out_st,_extra=extra
  1 +FUNCTION dustem_compute_stokes,p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,out_st=out_st,_extra=extra
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -42,34 +42,26 @@ IF keyword_set(help) THEN BEGIN
42 42 goto,the_end
43 43 ENDIF
44 44  
45   -IF not keyword_set(stp) THEN BEGIN
46   - stp=dustem_run(p_dim)
47   - dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),stp
  45 +IF not keyword_set(st) THEN BEGIN
  46 + st=dustem_run(p_dim)
  47 + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st
48 48 ENDIF
49 49  
  50 +out=0.
  51 +
  52 +If keyword_set(result) THEN result=st
  53 +
  54 +fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(st.sed).wav)*1.e20/1.e7
  55 +stI = st.sed.em_tot * fact ;This is Total intensity I
  56 +P = st.polsed.em_tot * fact ;This is Polarized intensity P
50 57  
51   -If keyword_set(result) THEN result=stp
52 58  
53   -fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(stp.polsed).wav)*1.e20/1.e7
54   -stI = stp.sed.em_tot * fact ;This is Total intensity I
55   -P = stp.polsed.em_tot * fact ;This is Polarized intensity P
56   -out=0.
57 59 Nwaves=(size(stI))[1]
58 60  
59 61 frac=P/stI
60 62 tes=where(finite(frac) eq 0)
61 63 frac(tes)=0.
62 64  
63   -;if Q and U have already been computed elsewhere
64   -; FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
65   -;
66   -; IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'MODIFY_DUST_POLAR') THEN BEGIN
67   -;
68   -; Q_spec=(*(*!dustem_plugin).(i))[*,1]
69   -; U_spec=(*(*!dustem_plugin).(i))[*,2]
70   -;
71   -; ENDIF
72   -; ENDFOR
73 65  
74 66  
75 67 FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
... ... @@ -84,9 +76,11 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
84 76 ENDIF
85 77  
86 78 ENDFOR
87   -
88 79  
89   -IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0.,Nwaves)
  80 +IF ~isa(Q_spec) && ~isa(U_spec) THEN BEGIN
  81 +polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0,Nwaves)
  82 +;stop
  83 +ENDIF
90 84  
91 85 FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
92 86  
... ... @@ -103,7 +97,6 @@ ENDFOR
103 97  
104 98 ;-----------------------------------------
105 99  
106   -;largeP=sqrt(StokesQ^2+StokesU^2)
107 100  
108 101 ;INITIALIZING THE STOKES SEDS
109 102 dustem_qsed = (*!dustem_data.sed).values * 0.
... ... @@ -127,11 +120,11 @@ IF count_polsed NE 0 THEN BEGIN ; INCLUDE COLOR CORRECTIONS OF (Q,U) for dust &
127 120 filter_names=((*!dustem_data.sed).filt_names)(ind_polsed)
128 121  
129 122 if isa(dustem_qsed) then begin
130   - sqsed=dustem_cc(stp.sed.wav,Q_spec,filter_names,cc=cc)
  123 + sqsed=dustem_cc(st.polsed.wav,Q_spec,filter_names,cc=cc)
131 124 dustem_qsed(ind_polsed)=sqsed
132 125 endif
133 126 if isa(dustem_used) then begin
134   - sused=dustem_cc(stp.sed.wav,U_spec,filter_names,cc=cc)
  127 + sused=dustem_cc(st.polsed.wav,U_spec,filter_names,cc=cc)
135 128 dustem_used(ind_polsed)=sused
136 129 endif
137 130 ;--------------------------------------------------------------------------
... ... @@ -145,10 +138,10 @@ ENDIF
145 138  
146 139 ind_spec=where((*!dustem_data.sed).filt_names EQ 'SPECTRUM',count_spec)
147 140 IF count_spec NE 0 THEN BEGIN
148   - if isa(dustem_qsed) then dustem_qsed(ind_spec)=interpol(Q_spec,stp.wav,(((*!dustem_data.sed).wav)(ind_spec)))
149   - if isa(dustem_used) then dustem_used(ind_spec)=interpol(U_spec,stp.wav,(((*!dustem_data.sed).wav)(ind_spec)))
  141 + if isa(dustem_qsed) then dustem_qsed(ind_spec)=interpol(Q_spec,st.polsed.wav,(((*!dustem_data.sed).wav)(ind_spec)))
  142 + if isa(dustem_used) then dustem_used(ind_spec)=interpol(U_spec,st.polsed.wav,(((*!dustem_data.sed).wav)(ind_spec)))
150 143 ENDIF
151   -out_st=stp
  144 +out_st=st
152 145  
153 146  
154 147 ;clean pointers
... ...
src/idl/dustem_fit_sed_polsed_readme.pro renamed to src/idl/dustem_fit_sed_qsed_used_readme.pro
  Diff comments   View file @607060e
1   -PRO dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,help=help
  1 +PRO dustem_fit_sed_qsed_used_readme,postcript=postcript,mode=mode,help=help
2 2  
3 3 ;This Readme describes how to fit total intensity and polarization SEDs jointly
4 4 ;Obviously, the dustem package must have been installed succesfully
... ... @@ -6,7 +6,7 @@ PRO dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,help=help
6 6  
7 7 ;+
8 8 ; NAME:
9   -; dustem_fit_sed_polsed_readme
  9 +; dustem_fit_sed_qsed_used_readme
10 10 ; PURPOSE:
11 11 ; This is an example of how to fit total intensity and polarization SEDs jointly
12 12 ; It is meant to be an example to follow when writing your own
... ... @@ -43,15 +43,16 @@ PRO dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,help=help
43 43 ;-
44 44  
45 45 IF keyword_set(help) THEN BEGIN
46   - doc_library,'dustem_fit_sed_polsed_readme'
  46 + doc_library,'dustem_fit_sed_qsed_used_readme'
47 47 goto,the_end
48 48 ENDIF
49 49  
  50 +
50 51 dustem_define_la_common
51 52  
52 53 ;=== initialise dustem
53 54 ;dustem_init,mode='COMPIEGNE_ETAL2010',/pol
54   -use_mode='G17_MODELA'
  55 +use_mode='G17_MODELA' ; try with model D
55 56  
56 57 dustem_init,mode=use_mode,/pol
57 58  
... ... @@ -61,6 +62,9 @@ help,!run_pol
61 62  
62 63 ;stop
63 64  
  65 +!dustem_verbose=1
  66 +!dustem_show_plot=1
  67 +
64 68  
65 69 ;filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
66 70 filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
... ... @@ -70,81 +74,135 @@ sed.filter=filters
70 74 sed.wave=dustem_filter2wav(filters)
71 75 sed.instru=dustem_filter2instru(filters)
72 76  
73   -;INITIALIZING IQU & P SEDS and associated errors
  77 +
  78 +
  79 +;INITIALIZING IQU SEDS and associated errors
74 80 sed.StokesI=sed.wave*0
75   -sed.StokesQ=sed.wave*0
76   -sed.StokesU=sed.wave*0
77   -;sed.largeP=sed.wave*0
  81 +;sed.StokesQ=sed.wave*0
  82 +;sed.StokesU=sed.wave*0
78 83 sed.sigmaII=sed.wave*0
79   -sed.sigmaQQ=sed.wave*0
80   -sed.sigmaUU=sed.wave*0
81   -;sed.sigma_largep=sed.wave*0
  84 +;sed.sigmaQQ=sed.wave*0
  85 +;sed.sigmaUU=sed.wave*0
82 86  
83 87  
84   -;st=dustem_set_data(ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI,qsed=qsed,used=used,qext=qext,uext=uext)
85 88  
86 89 ;INITIALIZING THE !dustem_data tags that will be internally used
87 90 st=dustem_set_data(sed=sed,rchi2_weight=rchi2_weight,f_HI=f_HI)
88 91  
89   -;st=dustem_set_data() ;this should be one of the tests
90 92  
91 93  
92 94 ;=== Set which parameters you want to fit
93 95 pd = [ $
94 96 '(*!dustem_params).G0', $ ;G0
95   - '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
96   - '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction
97   - '(*!dustem_params).grains(2).mdust_o_mh',$ ;PAH1 mass fraction
98   - '(*!dustem_params).grains(3).mdust_o_mh',$ ;amCBEx
  97 + '(*!dustem_params).grains(0).mdust_o_mh',$
  98 + '(*!dustem_params).grains(1).mdust_o_mh',$
  99 + ;'(*!dustem_params).grains(1).mdust_o_mh',$
  100 + '(*!dustem_params).grains(2).mdust_o_mh'$
99 101 ;'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude
100 102 ;'dustem_plugin_synchrotron_4', $ ;Synchrotron polarization angle
101   - 'dustem_plugin_modify_dust_polarization_2', $ ;Polarization angle applied to the Fortran dust model via a core plugin
102   - 'dustem_plugin_modify_spinning_polarization_1', $
103   - 'dustem_plugin_modify_spinning_polarization_2' $
  103 + ;'dustem_plugin_modify_dust_polarization_2', $ ;Polarization angle applied to the Fortran dust model via a core plugin
  104 + ;'dustem_plugin_modify_spinning_polarization_1', $
  105 + ;'dustem_plugin_modify_spinning_polarization_2' $
104 106 ]
105 107  
106   -;p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,68.,1.2,20]
107   -p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,1.6500E-04,68.,4e-3,20]
108   -iv=p_truth+[0.,8.e-4,8.e-4,8.e-4,8.e-4,-10.,2e-4,10] ;shifted from solution
  108 +;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,68.,0.01,30]
  109 +
  110 +p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04];,0.01,30.]
  111 +
  112 +;p_truth=[68.]
109 113  
110 114  
  115 +; fpd = [ $
  116 +; '(*!dustem_params).G0', $ ;G0
  117 +; '(*!dustem_params).grains(0).mdust_o_mh',$
  118 +; '(*!dustem_params).grains(1).mdust_o_mh',$
  119 +; '(*!dustem_params).grains(2).mdust_o_mh'$
  120 +; ; '(*!dustem_params).grains(3).mdust_o_mh',$
  121 +; ; 'dustem_plugin_modify_dust_polarization_2', $ ;Polarization angle applied to the Fortran dust model via a core plugin
  122 +; ; 'dustem_plugin_modify_spinning_polarization_1' $
  123 +; ]
  124 +; ;
  125 +; fiv=[1.,7.8000E-04,7.8000E-04,1.6500E-04];,68.,0.01]
  126 +
  127 +
  128 +fixed=[1,1,1,1,0];0,1,1]
  129 +
  130 +;dustem_init_parinfo,pd,p_truth;,fixed=fixed
  131 +
  132 +;dustem_init_fixed_params,fpd,fiv
  133 +
  134 +
  135 +iv=p_truth+[0.,5.000E-04,5.000E-04,5.000E-04];,5.e-4,1.] ;shifted from solution ,10.
  136 +
111 137  
112 138 Npar=n_elements(pd)
113   -ulimed=replicate(0,Npar)
114   -llimed=replicate(1,Npar)
115   -llims=replicate(0,Npar)
116 139  
117   -;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
118   -dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
119   -(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
120   -(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
  140 +ulimed=replicate(0,Npar);[0]
121 141  
  142 +;ulimed=[0,0,0,0,0,1,1,1]
  143 +;ulims=[0.,0.,0.,0.,0.,100.,1.,100.]
122 144  
123   -dustem_init_plugins,pd
  145 +llimed=replicate(1,Npar);[1]
  146 +llims=replicate(0,Npar);[10]
124 147  
125   -help,*!dustem_fit
  148 +dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims;,fixed=fixed
126 149  
127   -help,!run_pol
  150 +;dustem_init_fixed_params,fpd,fiv
128 151  
  152 +dustem_init_plugins,pd;, fpd=fpd
129 153  
130 154  
131   -sed.StokesI=dustem_compute_sed(p_truth,ssti) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
132   -sed.sigmaII=abs(sed.StokesI*0.000001)
133 155  
134   -;stop
  156 +sed.StokesI=dustem_compute_sed(p_truth,ste) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
  157 +sed.sigmaII=(sed.StokesI*0.3)^2
  158 +
  159 +; toto=dustem_compute_stokes(p_truth,ste,dustem_qsed,dustem_used) ;this procedure also allows for the extraction of the spectra
  160 +
  161 +; sed.StokesQ=dustem_qsed
  162 +; sed.sigmaQQ=(double((sed.StokesQ*0.03)^2)) > 1d-50
  163 +; ;tes=where(finite(sed.sigmaQQ) eq 0,coun)
  164 +; ;if coun ne 0 then sed.sigmaQQ(tes)=0.
  165 +
  166 +; sed.StokesU=dustem_used
  167 +; sed.sigmaUU=(double((sed.StokesU*0.03)^2)) > 1d-50
  168 +; ;tes=where(finite(sed.sigmaUU) eq 0,coun)
  169 +; ;if coun ne 0 then sed.sigmaUU(tes)=0.
  170 +
135 171  
  172 +;======================================================================================================================
136 173  
137   -;sed.largeP=dustem_compute_polsed(p_truth,sstp)
138   -;sed.sigma_largep=abs(sed.largeP*0.000001)
139 174  
140 175  
141   -toto=dustem_compute_stokes(p_truth,sstqu,dustem_qsed=dustem_qsed,dustem_used=dustem_used) ;this procedure also allows for the extraction of the spectra
  176 +; llimed=[1];replicate(1,Npar)
  177 +; llims=[1e-12];replicate(0,Npar)
  178 +
  179 +;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
  180 +
  181 +;dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims;,fixed=fixed
  182 +
142 183  
143   -sed.StokesQ=dustem_qsed
144   -sed.sigmaQQ=abs(sed.StokesQ*0.000001)
  184 +;dustem_init_plugins,pd,fpd=fpd
  185 +
  186 +
  187 +
  188 +;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
  189 +;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
  190 +
  191 +
  192 +;
  193 +; sed.StokesI=dustem_compute_sed(p_truth,ssti) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
  194 +; sed.sigmaII=(sed.StokesI*0.3)^2
  195 +;
  196 +;
  197 +; toto=dustem_compute_stokes(p_truth,sstqu,dustem_qsed,dustem_used) ;this procedure also allows for the extraction of the spectra
  198 +;
  199 +; sed.StokesQ=dustem_qsed
  200 +; sed.sigmaQQ=(double(sed.StokesQ)*0.4^2) > 1d-50
  201 +;
  202 +;
  203 +; sed.StokesU=dustem_used
  204 +; sed.sigmaUU=(double(sed.StokesU)*0.4^2) > 1d-50
145 205  
146   -sed.StokesU=dustem_used
147   -sed.sigmaUU=abs(sed.StokesU*0.000001)
148 206  
149 207 ;=== SET THE OBSERVATION STRUCTURE
150 208  
... ... @@ -153,16 +211,18 @@ ind=where(sed.wave GE 0.1)
153 211 st_bidon=dustem_set_data(sed=sed[ind],rchi2_weight=rchi2_weight,f_HI=f_HI) ; comment if second line is to be used
154 212  
155 213 ;=== RUN fit
156   -tol=1.e-16
157   -xtol=1.e-16
158   -Nitermax=30; just an initialization
  214 +tol=1.e-40
  215 +xtol=1.e-40
  216 +Nitermax=23; just an initialization
  217 +
  218 +
159 219 ;for i=0L,n_tags(!dustem_data)-1 do begin
160 220 ;if ptr_valid(!dustem_data.(i)) then Nitermax+=1
161 221 ;endfor
162 222 ;Nitermax=4 ;maximum number of iteration. This is the criterion which will stop the fit procedure
163 223  
164 224 xrange=[1.,2.e4]
165   -yrange=[1e-4,1.e3]
  225 +yrange=[1e-7,1.e3]
166 226  
167 227 loadct,13
168 228 ;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine
... ...
src/idl/dustem_fit_sed_readme.pro
  Diff comments   View file @607060e
... ... @@ -68,7 +68,7 @@ ENDIF
68 68 IF keyword_set(model) THEN BEGIN
69 69 use_model=strupcase(model)
70 70 ENDIF ELSE BEGIN
71   - use_model='COMPIEGNE_ETAL2010' ;Default is last dustem model
  71 + use_model='G17_MODELA';'COMPIEGNE_ETAL2010' ;Default is last dustem model
72 72 ENDELSE
73 73  
74 74 IF keyword_set(png) THEN BEGIN
... ...
src/idl/dustem_init.pro
  Diff comments   View file @607060e
... ... @@ -195,6 +195,7 @@ defsysv,'!dustem_filters',ptr_new() ;filter information (filled by dustem_fi
195 195 defsysv,'!dustem_verbose',0 ;1=verbose
196 196 defsysv,'!dustem_show_plot',1 ;0=show no plot
197 197  
  198 +
198 199 ;apparently not needed anymore. Was removed.
199 200 ;defsysv,'!ki',0. ;apparently needed by mpfit, but what is it ?
200 201  
... ... @@ -296,7 +297,7 @@ CASE !dustem_which OF
296 297 spawn,'rm -rf '+!dustem_dat+'/hcap'
297 298 spawn,'rm -rf '+!dustem_dat+'/oprop'
298 299 spawn,'rm -rf '+!dustem_dat+'/out'
299   - ;spawn,'mkdir '+!dustem_dat+'/plots' TO BE TESTED. NOT SURE ABOUT WHAT I AM DOING HERE. ASK JP ABOUT THIS
  300 + ;spawn,'rm -rf '+!dustem_dat+'/plots' TO BE TESTED. NOT SURE ABOUT WHAT I AM DOING HERE. ASK JP ABOUT THIS
300 301 ENDIF
301 302 END
302 303 ELSE: BEGIN
... ...
src/idl/dustem_init_fixed_params.pro
  Diff comments   View file @607060e
... ... @@ -22,23 +22,31 @@ ENDIF
22 22 (*!dustem_fit).param_init_values=ptr_new(fiv)
23 23 dustem_set_func_ind,fpd,fiv
24 24 ind=where(fpd NE '',count)
  25 +
  26 +
25 27 IF count NE 0 THEN BEGIN
26 28 dustem_set_params,[fiv]
27 29 ;dustem_create_functions,[fiv]/(*(*!dustem_fit).param_init_values)
  30 + ;dustem_init_plugins,fpd
28 31 dustem_activate_plugins,[fiv]/(*(*!dustem_fit).param_init_values)
29   - dustem_set_params,[fiv]
  32 + ;dustem_set_params,[fiv]
30 33 ENDIF
31 34  
  35 +
32 36 ;set the fixed system variables, used to keep memory of the fixed
33 37 ;parameters description and values
34 38 (*!dustem_fit).fixed_param_descs=ptr_new(fpd)
35 39 (*!dustem_fit).fixed_param_init_values=ptr_new(fiv)
36 40  
  41 +
37 42 ;put back in place fit parameter descriptions and values into system variables, if needed.
38 43 IF need_save EQ 1 THEN BEGIN
39 44 (*!dustem_fit).param_descs=ptr_new(pd)
40 45 (*!dustem_fit).param_init_values=ptr_new(iv)
41 46 (*!dustem_fit).param_func=ptr_new(find)
42   -ENDIF
  47 +ENDIF
  48 +
  49 +
  50 +
43 51  
44 52 END
... ...
src/idl/dustem_init_parinfo.pro
  Diff comments   View file @607060e
... ... @@ -32,7 +32,7 @@ IF keyword_set(up_limited) THEN BEGIN
32 32 ind=where(up_limited EQ 1,countup)
33 33 IF countup EQ 0 THEN goto,suite1
34 34 IF NOT keyword_set(up_limits) THEN BEGIN
35   - print,'up_limits must be set with up_limited'
  35 + message,'up_limits must be set with up_limited'
36 36 stop
37 37 ENDIF
38 38 IF countup NE 0 THEN BEGIN
... ... @@ -47,7 +47,7 @@ IF keyword_set(lo_limited) THEN BEGIN
47 47 ind=where(lo_limited EQ 1,countlo)
48 48 IF countlo EQ 0 THEN goto,suite2
49 49 IF NOT keyword_set(lo_limits) THEN BEGIN
50   - print,'lo_limits must be set with lo_limited'
  50 + message,'lo_limits must be set with lo_limited', /continue
51 51 stop
52 52 ENDIF
53 53 IF countlo NE 0 THEN BEGIN
... ...
src/idl/dustem_init_plugins.pro
  Diff comments   View file @607060e
1   -PRO dustem_init_plugins, pd
  1 +PRO dustem_init_plugins, pd, fpd=fpd
2 2  
3 3  
4 4 ;INITIALIZING THE SCOPES DATA ARRAYS--------------------
... ... @@ -9,6 +9,7 @@ plugin_names=['']
9 9 plugind_detect_string='dustem_plugin'
10 10  
11 11 Nplugins=0L
  12 +
12 13 for i=0L,n_elements(pd)-1 do begin
13 14 fi=strtrim(strmid(pd(i),0,strlen(plugind_detect_string)),2)
14 15 if fi eq plugind_detect_string then begin
... ... @@ -21,6 +22,32 @@ for i=0L,n_elements(pd)-1 do begin
21 22 ENDIF
22 23 ENDFOR
23 24  
  25 +
  26 +
  27 +IF KEYWORD_SET(fpd) THEN BEGIN
  28 +
  29 + for i=0L,n_elements(fpd)-1 do begin
  30 + fi=strtrim(strmid(fpd(i),0,strlen(plugind_detect_string)),2)
  31 + if fi eq plugind_detect_string then begin
  32 + if isa((*!dustem_fit).fixed_param_descs) then begin
  33 +
  34 + ftn = strmid((*(*!dustem_fit).fixed_param_descs)(i),0) ; String containing the name of the plugin and the keyword used (ie: dustem_create_continuum_2)
  35 +
  36 + endif else ftn = strmid((*(*!dustem_fit).param_descs)(i),0)
  37 +
  38 +
  39 + ii = strsplit(ftn,'_',count=countx)
  40 + ii = ii(countx-1)-1 ; Locating the last underscore to automate the extraction of the plugin's keyword
  41 + ftn = strmid(ftn,14,ii-14) ; String containing the name of the plugin without the associated keyword
  42 + plugin_names=[plugin_names,ftn]
  43 + Nplugins=Nplugins+1
  44 + ENDIF
  45 + ENDFOR
  46 +
  47 +
  48 +ENDIF
  49 +
  50 +
24 51 IF Nplugins EQ 0 THEN BEGIN
25 52 plugin_names=['NONE']
26 53 Nplugins=1
... ... @@ -63,10 +90,6 @@ IF Nplugins NE 0 THEN BEGIN
63 90  
64 91 ENDIF
65 92  
66   -
67   -;test=(n_tags(varvar) eq 1 && tag_exist(varvar,'varvar'))
68   -;if test then varvar=create_struct('varvar','ERROR IN DUSTEM_INIT_PLUGINS!')
69   -
70 93 defsysv, '!dustem_scope', ptr_new(varvar)
71 94 ;---------------------------------------------
72 95  
... ...
src/idl/dustem_mpfit_run.pro
  Diff comments   View file @607060e
... ... @@ -133,8 +133,8 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
133 133 / total(dustem_sed^2/(*!dustem_data.sed).sigma^2)
134 134  
135 135 print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
136   - (*!dustem_fit).chi2=chi2_sed
137   - (*!dustem_fit).rchi2=rchi2_sed
  136 + ;(*!dustem_fit).chi2=chi2_sed
  137 + ;(*!dustem_fit).rchi2=rchi2_sed
138 138  
139 139 ; Plot if needed
140 140 IF !dustem_show_plot NE 0 THEN BEGIN
... ... @@ -367,7 +367,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
367 367  
368 368 'QSED': BEGIN
369 369  
370   - toto = dustem_compute_stokes(p_dim,st,dustem_qsed=dustem_qsed,Q_spec=Q_spec)
  370 + toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
371 371 chi2_qsed =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2)
372 372 n_qsed = n_elements((*!dustem_data.qsed).values)
373 373 rchi2_qsed = chi2_qsed / n_qsed
... ... @@ -388,7 +388,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
388 388 titstq='Stokes Q Polarization Intensity'
389 389 ytitstq=textoidl('\nuQ_\nu/N_H (W/m^2/sr/H)')
390 390 xtit=textoidl('\lambda (\mum)')
391   - xr=[10,6e4]
  391 + xr=[10,2e4]
392 392 normpol=Q_spec*fact
393 393 normpol1=dustem_qsed*fact1
394 394 ;normpol1=(*!dustem_data.qsed).values*fact1
... ... @@ -419,7 +419,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
419 419  
420 420 'USED': BEGIN
421 421  
422   - toto = dustem_compute_stokes(p_dim,st,dustem_used=dustem_used,U_spec=U_spec)
  422 + toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
423 423 chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2)
424 424 n_used = n_elements((*!dustem_data.used).values)
425 425 rchi2_used = chi2_used / n_used
... ...
src/idl/dustem_old2new_sed_format.pro
  Diff comments   View file @607060e
... ... @@ -2,8 +2,12 @@ PRO dustem_old2new_sed_format
2 2  
3 3 dustem_init
4 4  
5   -dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
6   -file=dir+'Gal_composite_spectrum.xcat'
  5 +;dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
  6 +dir='/Users/ilyeschoubani/Downloads/test_wrap_NY/'
  7 +
  8 +file=dir+'SED_HD21.xcat'
  9 +
  10 +;file=dir+'Gal_composite_spectrum.xcat'
7 11  
8 12 str='cp '+file+' '+file+'.old'
9 13 spawn,str
... ...
src/idl/dustem_plot_fit_sed.pro
  Diff comments   View file @607060e
... ... @@ -115,7 +115,7 @@ use_lines=replicate(0,Ngrains)
115 115 norm = sed * 0. + 1
116 116  
117 117 xr=[1,6e4]
118   -yr=[1e-5,5e3]
  118 +yr=[1e-7,5e3]
119 119  
120 120 ;====== PLOT THE SED
121 121  
... ... @@ -198,7 +198,7 @@ cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot
198 198 ;----------------------------------------
199 199  
200 200 ;==== print the legend
201   -frmt0='(A32)'
  201 +frmt0='(A40)'
202 202 frmt1='(1E10.2)'
203 203 frmt2='(F7.2)'
204 204 use_legend_xpos=0.50
... ... @@ -234,7 +234,7 @@ IF keyword_set(res) THEN BEGIN
234 234 yypos=use_legend_ypos-(i-k)*use_legend_offset
235 235 xxpos=use_legend_xpos*0.26
236 236 ii = strsplit(string_name,'_',count=countx) & ii = ii(countx-1)-1 ; Locating the last underscore to automate the extraction of the plugin's keyword
237   - str=string(strmid(string_name,7,ii-7)+' ('+strmid(string_name,ii+1)+') '+ ' = ',format=frmt0)+string(res[i],format=frmt1)
  237 + str=string(strmid(string_name,7,ii-7)+' ('+strmid(string_name,ii+1)+') '+' = ',format=frmt0)+string(res[i],format=frmt1)
238 238  
239 239  
240 240 IF keyword_set(errors) THEN BEGIN
... ...
src/idl/dustem_plugin_modify_dust_polarization.pro
  Diff comments   View file @607060e
... ... @@ -38,19 +38,18 @@ IF keyword_set(help) THEN BEGIN
38 38 ENDIF
39 39  
40 40  
41   -scope='REPLACE_QSED+REPLACE_USED'
42   -;scope=''
43   -
44 41  
  42 +IF not KEYWORD_SET(st) THEN BEGIN
45 43 IF ptr_valid((*!dustem_fit).CURRENT_PARAM_VALUES) THEN BEGIN
46   - p_di=(*(*!dustem_fit).CURRENT_PARAM_VALUES)
47   - st=dustem_run(p_di)
  44 + p_di=(*(*!dustem_fit).CURRENT_PARAM_VALUES)
48 45  
49 46 ENDIF ELSE BEGIN
50   - p_di=(*(*!dustem_fit).PARAM_INIT_VALUES)
51   - st=dustem_run(p_di)
52   -
  47 + p_di=(*(*!dustem_fit).PARAM_INIT_VALUES)
53 48 ENDELSE
  49 +
  50 +st=dustem_run(p_di)
  51 +ENDIF
  52 +
54 53  
55 54 psi=0.
56 55 smallp=1.
... ... @@ -62,9 +61,11 @@ IF keyword_set(key) THEN BEGIN
62 61 ENDIF
63 62  
64 63  
65   -fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.e20/1.e7
66   -P=st.polsed.em_tot*fact
67   -I=st.sed.em_tot*fact
  64 +fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/((st.polsed).wav))*1.e20/1.e7
  65 +P=((st.polsed).em_tot)*fact
  66 +I=((st.sed).em_tot)*fact
  67 +
  68 +
68 69  
69 70 Nwaves=(size(I))[1]
70 71  
... ... @@ -81,6 +82,8 @@ out[*,0]=I
81 82 out[*,1]=Q
82 83 out[*,2]=U
83 84  
  85 +scope='REPLACE_QSED+REPLACE_USED'
  86 +
84 87 return, out
85 88  
86 89 the_end:
... ...
src/idl/dustem_plugin_modify_spinning_polarization.pro
  Diff comments   View file @607060e
... ... @@ -31,28 +31,32 @@ FUNCTION dustem_plugin_modify_spinning_polarization, key=key, val=val, scope=sco
31 31 ; This is a dustem plugin
32 32 ;-
33 33  
  34 +
34 35 IF keyword_set(help) THEN BEGIN
35 36 doc_library,'dustem_plugin_modify_spinning_polarization'
36 37 goto,the_end
37 38 ENDIF
38 39  
39 40  
40   -scope='ADD_QSED+ADD_USED'
  41 +
41 42  
42 43 IF ptr_valid((*!dustem_fit).CURRENT_PARAM_VALUES) THEN BEGIN
43 44 p_dii=(*(*!dustem_fit).CURRENT_PARAM_VALUES)
44   - st=dustem_run(p_dii)
45 45  
46 46 ENDIF ELSE BEGIN
47   - p_dii=(*(*!dustem_fit).PARAM_INIT_VALUES)
48   - st=dustem_run(p_dii)
49   -
50   - ENDELSE
  47 +
  48 +
  49 + p_dii=(*(*!dustem_fit).PARAM_INIT_VALUES)
  50 +
  51 +
  52 + ENDELSE
  53 +
  54 +st=dustem_run(p_dii)
51 55  
52   -;stop
53 56  
54 57 psi=0.
55 58 smallp=0. ;supposing the spinning dust is
  59 +
56 60 IF keyword_set(key) THEN BEGIN
57 61  
58 62 ind1=where(key EQ 1,count1)
... ... @@ -62,30 +66,42 @@ IF keyword_set(key) THEN BEGIN
62 66  
63 67 ENDIF
64 68  
  69 +print, 'smallp:'
  70 +print, smallp
65 71  
66   -fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(st.sed).wav)*1.e20/1.e7
  72 +print, 'psi'
  73 +print, psi
67 74  
68 75  
69   -SpinningI=(st.sed.em_grain_1)+(st.sed.em_grain_2)*fact
70 76  
71   -;760
72   -SpinningI[0:760]=0.
  77 +fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/((st.polsed).wav))*1.e20/1.e7
73 78  
  79 +filee=!dustem_dat+'out/SPIN.RES'
74 80  
75   -Nwaves=(size(SpinningI))[1]
  81 +readcol,filee,waves,em_g1,em_g2,em_g3,em_tot,format='D,D,D,D,D',delimiter=' ',silent=silent,skipline=8;comment='#'
  82 +spawn,'rm -rf '+!dustem_dat+'out/*SPIN*'
76 83  
  84 +Nwaves=n_elements(waves)
  85 +out=fltarr(Nwaves,3)
77 86  
78   -polar_ippsi2iqu,SpinningI,SpinningQ,SpinningU,replicate(smallp,Nwaves),replicate(psi,Nwaves)
79 87  
80   -out=fltarr(Nwaves,3)
  88 +out[*,0]=em_tot*fact
  89 +
  90 +polar_ippsi2iqu,out[*,0],SpinningQ,SpinningU,replicate(smallp,Nwaves),replicate(psi,Nwaves)
81 91  
82   -out[*,0]=SpinningI
83 92 out[*,1]=SpinningQ
84 93 out[*,2]=SpinningU
85 94  
86 95  
  96 +
  97 +scope='ADD_QSED+ADD_USED'
  98 +
  99 +
87 100 return, out
88 101  
  102 +
  103 +;clean pointers
  104 +heap_gc
89 105 the_end:
90 106  
91 107  
... ...
src/idl/dustem_read_qabs_lv.pro
  Diff comments   View file @607060e
... ... @@ -67,7 +67,7 @@ Nlines=n_elements(qsca_str)
67 67  
68 68 qabs_values=fltarr(Nlines,Ncol)
69 69 qsca_values=fltarr(Nlines,Ncol)
70   -FOR i=0L,Nlines-1 DO BEGIN
  70 +FOR i=0L,Nlines-2 DO BEGIN
71 71 qabs_values(i,*)=str_sep(strcompress(strtrim(qabs_str(i),2)),' ')
72 72 qsca_values(i,*)=str_sep(strcompress(strtrim(qsca_str(i),2)),' ')
73 73 ENDFOR
... ...
src/idl/dustem_read_spin.pro
  Diff comments   View file @607060e
... ... @@ -14,7 +14,7 @@ first_char='#'
14 14 Ncomments=0L
15 15 WHILE first_char EQ '#' DO BEGIN
16 16 readf,unit,str
17   - first_char=strmid(str,0,1)
  17 + first_char=strmid(str,0,1)
18 18 comments=[comments,str]
19 19 IF first_char EQ '#' THEN Ncomments=Ncomments+1
20 20 ENDWHILE
... ...
src/idl/dustem_set_data.pro
  Diff comments   View file @607060e
1   -FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI
  1 + FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -210,9 +210,6 @@ ENDIF
210 210 ; ;VGe
211 211 ; ENDIF
212 212  
213   -
214   -
215   -
216 213 ;I'm not so sure about the extinction
217 214  
218 215  
... ... @@ -297,8 +294,6 @@ ENDIF
297 294 !fit_rchi2_weight.ext=1.
298 295  
299 296  
300   -
301   -
302 297 IF NP NE 0 THEN !fit_rchi2_weight.polsed=1.
303 298 IF NQ NE 0 THEN !fit_rchi2_weight.qsed=1.
304 299 IF NU NE 0 THEN !fit_rchi2_weight.used=1.
... ...
src/idl/dustem_set_func_ind.pro
  Diff comments   View file @607060e
... ... @@ -32,6 +32,7 @@ PRO dustem_set_func_ind,pd,iv
32 32 ; Written by NF,JPB,DP Jan-2007
33 33 ;-
34 34  
  35 +
35 36 nparam = n_elements(*(*!dustem_fit).param_descs)
36 37  
37 38 ;!func_ind=ptr_new(fltarr(nparam))
... ...
src/idl/dustem_set_params.pro
  Diff comments   View file @607060e
... ... @@ -30,7 +30,7 @@ FOR i=0L,Nparams-1 DO BEGIN
30 30 END
31 31 'PLUGIN': BEGIN
32 32  
33   - IF !dustem_verbose NE 0 THEN message,str,/info
  33 + ;IF !dustem_verbose NE 0 THEN message,str,/info
34 34 ;length=strlen((*(*!dustem_fit).param_descs)(i))
35 35  
36 36 strr = strmid((*(*!dustem_fit).param_descs)(i),14)
... ...