OV-GSO-DC / dustem-wrapper_idl

Commit b74c4452532fcca2650381651a2a6f4ce957edc4

Authored by Ilyes Choubani
1 parent 73acf100
Exists in master

Treating extinction fitting + added new procedures. Moving on to extinction plotting

Showing 7 changed files with 346 additions and 250 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_compute_ext.pro 0 → 100755
  Diff comments   View file @b74c445
... ... @@ -0,0 +1,103 @@
  1 +FUNCTION dustem_compute_ext,p_dim,st,EXT_spec,out_st=out_st,_extra=extra,wave0=wave0
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_compute_sed
  6 +; PURPOSE:
  7 +; Computes an EXT curve from a given Dustem spectrum
  8 +; CATEGORY:
  9 +; Dustem
  10 +; CALLING SEQUENCE:
  11 +; sed=dustem_compute_ext(p_dim[,st=][,_extra=][,/help])
  12 +; INPUTS:
  13 +; p_dim = parameter values (used only if st is not provided)
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; st = Dustem output structure, if provided, overrides p_dim
  16 +; OUTPUTS:
  17 +; ext = computed EXT for spectrum points in in !dustem_data (extinction)
  18 +; OPTIONAL OUTPUT PARAMETERS:
  19 +; out_st = same as ext
  20 +; ACCEPTED KEY-WORDS:
  21 +; help = If set, print this help
  22 +; COMMON BLOCKS:
  23 +; None
  24 +; SIDE EFFECTS:
  25 +; None
  26 +; RESTRICTIONS:
  27 +; The dustem idl wrapper must be installed
  28 +;-
  29 +
  30 +IF keyword_set(help) THEN BEGIN
  31 + doc_library,'dustem_compute_ext'
  32 + dustem_sed=0.
  33 + goto,the_end
  34 +ENDIF
  35 +
  36 +IF not keyword_set(st) THEN BEGIN
  37 +
  38 + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values) ; what about the 0/0 division case?
  39 + st=dustem_run(p_dim)
  40 +
  41 +ENDIF
  42 +
  43 +if not keyword_set(wave0) then wave0 = 0.
  44 +
  45 +
  46 +EXT_spec = st.ext.ext_tot * (*!dustem_HCD)/1.0e21 ;For any column density. Default is 10^20 H/cm^2
  47 +
  48 +
  49 +
  50 +;ADDING PLUGIN TO SPECTRUM----------------
  51 +scopes=tag_names((*!dustem_plugin))
  52 +IF scopes[0] NE 'NONE' THEN BEGIN
  53 +;IF ptr_valid(!dustem_plugin) THEN BEGIN
  54 + for i=0L,n_tags(*!dustem_plugin)-1 do begin
  55 + if total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) eq 'ADD_EXT') then EXT_spec+=(*(*!dustem_plugin).(i).spec)[*,0]
  56 + endfor
  57 +ENDIF
  58 +;------------------------------------------
  59 +
  60 +if not isarray(EXT_spec) THEN stop
  61 +
  62 +;COMPUTE THE MODEL SED
  63 +
  64 +dustem_ext = (*!dustem_data.ext).values * 0.
  65 +
  66 +;####For now no color corrections are performed (spectrum points)####
  67 +
  68 +; IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
  69 +; message,'DOING color correction calculations',/info
  70 +; ENDIF ELSE BEGIN
  71 +; message,'SKIPPING color correction calculations: we take previous values of cc',/info
  72 +; ENDELSE
  73 +
  74 +; ind_sed=where((*!dustem_data.sed).filt_names NE 'SPECTRUM' and (*!dustem_data.sed).wav gt wave0,count_sed)
  75 +; IF count_sed NE 0 THEN BEGIN
  76 +; filter_names=((*!dustem_data.sed).filt_names)[ind_sed]
  77 +; ssed=dustem_cc(st.sed.wav,SED_spec,filter_names,cc=cc)
  78 +; dustem_sed[ind_sed]=ssed
  79 +; ENDIF ELSE BEGIN
  80 +; print,"**** WARNING **** NO FILTER defined in DATA for COLOUR CORRECTION"
  81 +; ENDELSE
  82 +
  83 +;For spectrum data points, interpolate in log-log
  84 +;Linear interpolation leads to wrong values, in particular where few
  85 +;wavelengths points exist in the model (long wavelengths). (done automatically because the dustem grid is sampled in already log-log space)
  86 +ind_spec=where((*!dustem_data.ext).filt_names EQ 'SPECTRUM',count_spec)
  87 +IF count_spec NE 0 THEN BEGIN
  88 + dustem_ext[ind_spec]=interpol(EXT_spec,st.ext.wav,(((*!dustem_data.ext).wav)[ind_spec]))
  89 +ENDIF
  90 +out_st=st
  91 +
  92 +;clean pointers
  93 +heap_gc
  94 +
  95 +the_end:
  96 +
  97 +;VG : error in cc IRAS1 when no PAH
  98 +; j=where(dustem_sed ne dustem_sed,count)
  99 +; if count gt 0 then dustem_sed[j]=0.
  100 +;stop
  101 +RETURN,dustem_ext
  102 +
  103 +END
... ...
src/idl/dustem_compute_polext.pro
  Diff comments   View file @b74c445
1   -FUNCTION dustem_compute_polext ,p_dim,sti,_extra=extra,out_st=out_st,dustem_qext,dustem_uext,Q_ext,U_ext
  1 +FUNCTION dustem_compute_polext ,p_dim,sti,POLEXT_spec,_extra=extra,out_st=out_st,dustem_qext,dustem_uext,Q_ext,U_ext
2 2  
3 3 ;CREATED AS A CONSEQUENCE OF THE INCLUSION OF THE PLUGINS
4 4  
... ... @@ -10,46 +10,21 @@ IF not keyword_set(sti) THEN BEGIN
10 10 sti=dustem_run(p_dim)
11 11 ENDIF
12 12  
13   -If ptr_valid(*!dustem_plugin) THEN BEGIN
14   - ;ADDING PLUGIN TO SPECTRUM----------------
15   - ;if n_tags(!dustem_data.polext) gt 1 then begin
16   - for i=0L,n_tags(*!dustem_scope)-1 do begin
17   - if total(strsplit((*(*!dustem_scope).(i)),'+',/extract) eq 'ADD_POLEXT') then sti.polext.ext_tot+=(*(*!dustem_plugin).(i))
18   - endfor
19   - ;endif
20   - ;------------------------------------------
21   -
22   -; ;=============================THIS BLOCK IS WRONG YOU HAVE TO CHANGE IT=============================
23   -;
24   -;
25   -; ;MODIFYING SPECTRUM VIA PLUGIN-------------- THIS BLOCK NEEDS TO BE UNDER THE DUSTEM_POLSED LINE
26   -;
27   -; for i=0L,n_tags(*!dustem_scope)-1 do begin
28   -;
29   -; if total(strsplit((*(*!dustem_scope).(i)),'+',/extract) eq 'MODIFY_POLEXT') then begin
30   -;
31   -; ;making sure it's the stokes plugin that gets used here
32   -;
33   -; iidx=where(tag_names(*!dustem_scope) eq 'STEXT')
34   -;
35   -;
36   -; IF iidx eq i THEN BEGIN ;not sure about this condition but you get the spirit
37   -;
38   -; Q_ext=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/2-1)
39   -; U_ext=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:*)
40   -; dustem_qext=interpol(Q_ext,sti.polext.wav,(*!dustem_data.polext).wav)
41   -; dustem_uext=interpol(U_ext,sti.polext.wav,(*!dustem_data.polext).wav)
42   -;
43   -; endif
44   -; endif
45   -; endfor
46   -; ;endif
47   -; ;-------------------------------------------
48   -; ;=======================================================================================================
49 13  
  14 +POLEXT_spec = sti.polext.ext_tot * (*!dustem_HCD)/1.0e21 ;
  15 +
  16 +
  17 +
  18 +;ADDING PLUGIN TO SPECTRUM----------------
  19 +scopes=tag_names((*!dustem_plugin))
  20 +IF scopes[0] NE 'NONE' THEN BEGIN
  21 +;IF ptr_valid(!dustem_plugin) THEN BEGIN
  22 + for i=0L,n_tags(*!dustem_plugin)-1 do begin
  23 + if total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) eq 'ADD_EXT') then POLEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,0]
  24 + endfor
50 25 ENDIF
51 26  
52   -dustem_polext = interpol(sti.polext.ext_tot,sti.polext.wav,(*!dustem_data.polext).wav)
  27 +dustem_polext = interpol(POLEXT_spec,sti.polext.wav,(*!dustem_data.polext).wav)
53 28  
54 29  
55 30 return, dustem_polext
... ...
src/idl/dustem_compute_stokes.pro
  Diff comments   View file @b74c445
1 1 FUNCTION dustem_compute_stokes,p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em,out_st=out_st,_extra=extra
2 2  
3   -;,SP_spec
4   -
5 3 ;+
6 4 ; NAME:
7 5 ; dustem_compute_stokes
... ... @@ -78,7 +76,7 @@ ENDFOR
78 76  
79 77 IF ~isa(Q_spec) && ~isa(U_spec) THEN BEGIN
80 78 polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(!dustem_psi,Nwaves);PSI_spec
81   - PSI_spec = 0.5*atan(U,Q)/degtorad
  79 + PSI_spec = 0.5*atan(U_spec,Q_spec)/degtorad
82 80 ENDIF
83 81  
84 82 FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
... ...
src/idl/dustem_compute_stokext.pro 0 → 100755
  Diff comments   View file @b74c445
... ... @@ -0,0 +1,161 @@
  1 +FUNCTION dustem_compute_stokext,p_dim,st,dustem_qext,dustem_uext,QEXT_spec,UEXT_spec,PSIEXT_spec,dustem_psi_ext,out_st=out_st,_extra=extra
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_compute_stokext
  6 +; PURPOSE:
  7 +; Provides Stokes Q and Stokes U "extinction" optical depths for a given dustem spectrum and wrapper SED
  8 +; CATEGORY:
  9 +; Dustem
  10 +; CALLING SEQUENCE:
  11 +; out=dustem_compute_stokext(p_dim[,st=][,_extra=][,/help])
  12 +; INPUTS:
  13 +; p_dim = parameter values
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; stp = Dustem structure
  16 +; OUTPUTS:
  17 +;
  18 +;
  19 +;
  20 +;
  21 +; OPTIONAL OUTPUT PARAMETERS:
  22 +; stp = Dustem output structure
  23 +; ACCEPTED KEY-WORDS:
  24 +; help = If set, print this help
  25 +; COMMON BLOCKS:
  26 +; None
  27 +; SIDE EFFECTS:
  28 +; None
  29 +; RESTRICTIONS:
  30 +; The dustem idl wrapper must be installed
  31 +; PROCEDURE:
  32 +; None
  33 +; EXAMPLES
  34 +;
  35 +; MODIFICATION HISTORY:
  36 +; Written by J.-Ph. Bernard
  37 +; see evolution details on the dustem cvs maintained at CESR
  38 +; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
  39 +
  40 +IF keyword_set(help) THEN BEGIN
  41 + doc_library,'dustem_compute_stokext'
  42 + out=0.
  43 + goto,the_end
  44 +ENDIF
  45 +
  46 +IF not keyword_set(st) THEN BEGIN
  47 + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values)
  48 + st=dustem_run(p_dim)
  49 +ENDIF
  50 +
  51 +out=0. ; no specific output is required.
  52 +
  53 +If keyword_set(result) THEN result=st
  54 +
  55 +degtorad = !pi/180
  56 +
  57 +EXT_spec = st.ext.ext_tot * (*!dustem_HCD)/1.d21 ;This is Total extinction (for a default value of 10^20 NH/cm^2)
  58 +POLEXT_spec = st.polext.ext_tot * (*!dustem_HCD)/1.d21 ;This is Polarized extinction
  59 +
  60 +Nwaves=(size(EXT_SPEC))[1]
  61 +
  62 +frac=POLEXT_spec/EXT_spec
  63 +tes=where(finite(frac) eq 0)
  64 +frac(tes)=0.
  65 +
  66 +FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
  67 +
  68 + IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'REPLACE_POLEXT') THEN BEGIN
  69 +
  70 + QEXT_spec = (*(*!dustem_plugin).(i).spec)[*,1]
  71 + UEXT_spec = (*(*!dustem_plugin).(i).spec)[*,2]
  72 + PSIEXT_spec = (*(*!dustem_plugin).(i).spec)[*,4]
  73 +
  74 + ENDIF
  75 +
  76 +ENDFOR
  77 +
  78 +IF ~isa(QEXT_spec) && ~isa(UEXT_spec) THEN BEGIN
  79 + polar_ippsi2iqu,EXT_spec,QEXT_spec,UEXT_spec,frac,replicate(!dustem_psi_ext,Nwaves);PSI_spec
  80 + PSIEXT_spec = 0.5*atan(UEXT_spec,QEXT_spec)/degtorad
  81 +ENDIF
  82 +
  83 +FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
  84 +
  85 + IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'ADD_POLEXT') THEN BEGIN
  86 +
  87 + QEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,1]
  88 + UEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,2]
  89 + PSIEXT_spec+= (*(*!dustem_plugin).(i).spec)[*,4]
  90 +
  91 + ENDIF
  92 +
  93 +ENDFOR
  94 +;-----------------------------------------
  95 +
  96 +;INITIALIZING THE STOKES SEDS
  97 +dustem_qext = (*!dustem_data.qext).values * 0.
  98 +dustem_uext = (*!dustem_data.uext).values * 0.
  99 +
  100 +if not isarray(EXT_spec) THEN stop
  101 +
  102 +
  103 +;COLOR CORRECTION NOT PERFORMED YET. ONLY SPECTRUM DATA POINTS ARE CONSIDERED FOR NOW
  104 +;
  105 +; ;Performing color corrections on the dustem spectra and generating the SEDs at the given filters
  106 +;
  107 +; IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
  108 +; message,'DOING color correction calculations',/info
  109 +; ENDIF ELSE BEGIN
  110 +; message,'SKIPPING color correction calculations',/info
  111 +; ENDELSE
  112 +;
  113 +; ind_qsed=where((*!dustem_data.qsed).filt_names NE 'SPECTRUM',count_qsed)
  114 +; ind_used=where((*!dustem_data.used).filt_names NE 'SPECTRUM',count_used)
  115 +;
  116 +;
  117 +; IF count_qsed NE 0 THEN BEGIN
  118 +; filter_names=((*!dustem_data.qsed).filt_names)(ind_qsed)
  119 +; if isa(dustem_qsed) then begin
  120 +; sqsed=dustem_cc(st.polsed.wav,Q_spec,filter_names,cc=cc)
  121 +; dustem_qsed(ind_qsed)=sqsed
  122 +; endif
  123 +; ENDIF
  124 +;
  125 +; IF count_used NE 0 THEN BEGIN
  126 +; filter_names=((*!dustem_data.used).filt_names)(ind_used)
  127 +; if isa(dustem_used) then begin
  128 +; sused=dustem_cc(st.polsed.wav,U_spec,filter_names,cc=cc)
  129 +; dustem_used(ind_used)=sused
  130 +; endif
  131 +; ENDIF
  132 +
  133 +
  134 +;For spectrum data points, interpolate in log-log
  135 +;Linear interpolation leads to wrong values, in particular where few
  136 +;wavelengths points exist in the model (long wavelengths).
  137 +
  138 +ind_qspec=where((*!dustem_data.qext).filt_names EQ 'SPECTRUM',count_qspec)
  139 +ind_uspec=where((*!dustem_data.uext).filt_names EQ 'SPECTRUM',count_uspec)
  140 +
  141 +IF count_qspec NE 0 THEN BEGIN
  142 + if isa(dustem_qext) then dustem_qext(ind_qspec)=interpol(QEXT_spec,st.polext.wav,(((*!dustem_data.qext).wav)(ind_qspec)))
  143 +ENDIF
  144 +
  145 +IF count_uspec NE 0 THEN BEGIN
  146 + if isa(dustem_uext) then dustem_uext(ind_uspec)=interpol(UEXT_spec,st.polext.wav,(((*!dustem_data.uext).wav)(ind_uspec)))
  147 +ENDIF
  148 +
  149 +out_st=st
  150 +
  151 +;generating interpolates for Psi_ext
  152 +dustem_psi_ext = 0.5*atan(dustem_uext,dustem_qext)/degtorad
  153 +
  154 +;clean pointers
  155 +heap_gc
  156 +
  157 +the_end:
  158 +
  159 +RETURN,out
  160 +
  161 +END
... ...
src/idl/dustem_init.pro
  Diff comments   View file @b74c445
... ... @@ -131,7 +131,7 @@ dustem_fit_st={data:ptr_new(), $ ;,wavelength:ptr_new(), ;because wavelength arr
131 131 defsysv, '!dustem_fit', ptr_new(dustem_fit_st) ;Data to fit
132 132  
133 133 if keyword_set(polarization) then begin
134   - defsysv, '!dustem_data', { $ ;Data to fit
  134 + defsysv, '!dustem_data', { $ ;Structure contaning the data to fit
135 135 sed: ptr_new(), $ ;
136 136 ext: ptr_new(), $ ;
137 137 qsed: ptr_new(), $ ;
... ... @@ -149,7 +149,7 @@ if keyword_set(polarization) then begin
149 149  
150 150 rchi2_weight = {sed: 0.,ext: 0.,qsed:0.,used:0.,polsed:0.,polfrac:0.,psi_em:0.,qext:0.,uext:0.,polext:0.,psi_ext:0.}
151 151  
152   - defsysv,'!dustem_show', { $ ;Data to fit
  152 + defsysv,'!dustem_show', { $ ;Structure containing all the data (to be fitted or not)
153 153 sed: ptr_new(), $ ;
154 154 ext: ptr_new(), $ ;
155 155 qsed: ptr_new(), $ ;
... ... @@ -166,13 +166,13 @@ if keyword_set(polarization) then begin
166 166  
167 167 endif else begin
168 168  
169   - defsysv, '!dustem_data', { $ ;Data to fit
  169 + defsysv, '!dustem_data', { $
170 170 sed: ptr_new(), $
171 171 ext: ptr_new() $
172 172 }
173 173 rchi2_weight = {sed: 0.,ext: 0.}
174 174  
175   - defsysv, '!dustem_show', { $ ;Data to fit
  175 + defsysv, '!dustem_show', { $
176 176 sed: ptr_new(), $
177 177 ext: ptr_new() $
178 178 }
... ... @@ -226,7 +226,11 @@ defsysv, '!dustem_mlog', 0 ;necessary for the plotting of negative values o
226 226 ;###################################################################
227 227  
228 228  
  229 +;These two sysvars allow for the polarization angle (whether in emission or extinction) to
  230 +;be updated in procedures (non-plugins) that create polarization data such as
  231 +;dustem_compute_stokes and dustem_compute_stokext
229 232 defsysv, '!dustem_psi', 0.
  233 +defsysv, '!dustem_psi_ext', 0.
230 234  
231 235  
232 236  
... ...
src/idl/dustem_mpfit_data.pro
  Diff comments   View file @b74c445
... ... @@ -56,8 +56,7 @@ use_Nitermax=3 ;Default maximum number of iterations
56 56 IF keyword_set(tol) THEN use_tol=tol
57 57 IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
58 58  
59   -; Build array concatenating observational datas
60   -; Order : that of !dustem_data (sed, ext, polext, polsed)
  59 +; Build array concatenating observational data
61 60  
62 61 if !run_pol then begin
63 62  
... ... @@ -68,8 +67,6 @@ if !run_pol then begin
68 67 if ctestpol ne 0 then begin
69 68 match,ind_data_tag,data_ind,ind1,ind2 ;not using match2 because when matching indices, if the first is 0 it doesn't count which is wrong.
70 69 ;I also only need the points that match.
71   - ;testind1 = where(ind1 ne -1,ctind1)
72   - ;if ctind1 ne 0 then ind_data_tag = ind_data_tag[testind1]
73 70 ind_data_tag = ind_data_tag[ind1]
74 71 count_data_tag = n_elements(ind_data_tag)
75 72 endif
... ...
src/idl/dustem_mpfit_run.pro
  Diff comments   View file @b74c445
... ... @@ -95,7 +95,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
95 95  
96 96 CASE tagnames(ind_data_tag(i_tag)) OF
97 97  
98   - 'SED' : BEGIN ;UPDATED. BUT ISSUE WITH THE CONCATENATION TO THE DUSTEM_OUT ARRAY AS OUTLINED BELOW
  98 + 'SED' : BEGIN
99 99  
100 100 dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
101 101 chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
... ... @@ -113,63 +113,33 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
113 113 / total(dustem_sed^2/(*!dustem_data.sed).sigma^2)
114 114  
115 115 print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
116   - (*!dustem_fit).chi2=chi2_sed
117   - (*!dustem_fit).rchi2=rchi2_sed
118   -
119   - ;I think this test isn't the correct one because when we want to show the SED data and not fit it we have an issue
120   - ;with the current hiding/showing formalism!!!!
  116 +
  117 + ;These two lines were put here so that the chi2 is updated everytime we plot
  118 + ;Moving it to the bottom so that all chi2 are taken into account
  119 +; (*!dustem_fit).chi2=chi2_sed
  120 +; (*!dustem_fit).rchi2=rchi2_sed
121 121  
122 122 dustem_out = [ dustem_out, dustem_sed ]
123 123  
124 124 END
125 125  
126   - 'EXT' : BEGIN ;NOT UPDATED
127   -
128   - ;ADDING PLUGIN TO SPECTRUM----------------
129   - IF ptr_valid(*!dustem_plugin) THEN BEGIN
130   - for i=0L,n_tags(*!dustem_scope)-1 do begin
131   - if total(strsplit(*(*!dustem_scope).(i),'+',/extract) eq 'ADD_EXT') then st.ext.ext_tot+=(*(*!dustem_plugin).(i))
132   - endfor
133   - endif
134   - ;------------------------------------------
135   -
136   - dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
137   -
  126 + 'EXT' : BEGIN
  127 + dustem_ext = dustem_compute_ext(p_dim,st,EXT_spec)
138 128 chi2_ext = total(((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2)
139 129 n_ext = n_elements((*!dustem_data.ext).values)
140 130 rchi2_ext = chi2_ext / n_ext
141   - wrchi2_ext = rchi2_ext * !fit_rchi2_weight.ext
142   - print,"chi2 EXT = ",chi2_ext," rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext
143   - ;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2
144   -
145   - ; Plot if needed
146   -; IF !dustem_show_plot NE 0 THEN BEGIN
147   -;
148   -; win+=1
149   -; ;dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[0.5,40],/xsty,win=win
150   -; xr=[0.3,40]
151   -; yr=[1.e-26,1.e-21]
152   -; dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
153   -; dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
154   -;
155   -; win+=1
156   -; ;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win
157   -; xr=[0.2,10]
158   -; yr=[0,3e-21]
159   -; dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
160   -; dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
161   -;
162   -; ENDIF
163   -
164   - ;dustem_plot_ext,st,/print_Rv ;commented this recently
  131 + wrchi2_ext = rchi2_ext * !fit_rchi2_weight.ext
  132 + print,"chi2 EXT = ",chi2_ext," rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext
  133 + ;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2
165 134  
166 135 dustem_out = [ dustem_out, dustem_ext ]
167 136  
168 137 END
169 138  
170   - 'POLEXT': BEGIN ;NOT UPDATED
  139 + 'POLEXT': BEGIN
171 140  
172 141  
  142 + ;Haven't investigated !run_circ=1 yet...
173 143 ;---------------the procedure within this block has not been updated for a vert long time
174 144 IF !dustem_show_plot NE 0 and !run_circ then begin
175 145 win+=1
... ... @@ -179,86 +149,16 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
179 149 ;--------------------------------------------------------------------------------------
180 150  
181 151  
182   -
183 152 IF !run_lin then begin
184 153 ;compute_polext needs updating. It's been months.
185   - dustem_polext = dustem_compute_polext(p_dim,st,dustem_qext,dustem_uext,Q_ext,U_ext)
  154 + dustem_polext = dustem_compute_polext(p_dim,st,POLEXT_spec)
186 155 chi2_polext = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2)
187   -; n_polext = n_elements((*!dustem_data.polext).values)
188   -; rchi2_polext = chi2_polext / n_polext
189   -; wrchi2_polext = rchi2_polext * !fit_rchi2_weight.polext
190   -; print,"chi2 POLEXT = ",chi2_polext," rchi2 POLEXT = ",rchi2_polext;," weighted rchi2 POLEXT=",wrchi2_polext
191   - ;print,"chi2 per wl: ",((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2
192   -
193   - ; Fit in (pmax,lmax,K)
194   - x=st.polext.wav
195   - logx=alog(x)
196   - logy=alog(st.polext.ext_tot)
197   -
198   - ; Whittet+1992 (Table 1)
199   - Uband = 0.36
200   - Bband = 0.43
201   - Vband = 0.55
202   - Rband = 0.63
203   - Iband = 0.78
204   - Jband = 1.21
205   - Hband = 1.64
206   - Kband = 2.04
207   - bands = [Uband,Bband,Vband,Rband,Iband,Jband,Hband]
208   - nbands = n_elements(bands)
209   - jband = indgen(nbands)*0
210   - for k=0,nbands-1 do begin
211   - xx=min(abs(st.polext.wav-bands(k)),j)
212   - jband(k) = j
213   - endfor
214   - ;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
215   - res = poly_fit(logx(jband),logy(jband),2)
216   -
217   - print
218   - print,"Serkowski Fit with free K over ",string(st.polext(jband).wav,format='(100F5.2)')
219   - print,"-------------------------"
220   - K = -res(2)
221   - lmax = exp(res(1)/(2*K))
222   - pmax = exp(res(0)+K*alog(lmax)^2)
223   - print,"* lmax = ",lmax," micron"
224   - print,"* K = ",K
225   - print,"* pmax = ",pmax," for NH=1e21"
226   - print
227   -
228   - ; Plot if needed
229   - IF !dustem_show_plot NE 0 THEN BEGIN
230   -
231   -; ; Polarization fraction in the UV-IR
232   -; win+=1
233   -;
234   -; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,50],/xsty,win=win;multi=0,ps=filename
235   -; ;
236   - ; Polarisation curve (Serkowski)
237   - win+=1
238   - ;xr=[0.3,1e4]
239   - xr=[1e-4,3.5]
240   - yr=[1e-30,5e-24]
241   -
242   - dustem_plot_polext,st,p_dim,dustem_polext,aligned=st_model.align.grains.aligned,win=win
243   -
244   -; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,/donotclose,multi=1
245   -; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,multi=2
246   -
247   -; win+=1
248   -; dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
249   -;
250   -
251   - ENDIF
252   -
253   - ;dustem_plot_polar,st,/print_ratio,cont=cont
254   -
255   - dustem_out = [ dustem_out, dustem_polext ]
256 156  
257 157 ENDIF
258 158  
259 159 END
260 160  
261   - 'QSED': BEGIN ;UPDATED
  161 + 'QSED': BEGIN
262 162  
263 163 toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em)
264 164  
... ... @@ -277,9 +177,8 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
277 177  
278 178 END
279 179  
280   - 'USED': BEGIN ;UPDATED
281   - ;THIS IMPLIES THAT QSED HAS BEEN RAN BEFORE
282   - ;toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
  180 + 'USED': BEGIN
  181 +
283 182 chi2_used = total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2)
284 183 n_used = n_elements((*!dustem_data.used).values)
285 184 rchi2_used = chi2_used / n_used
... ... @@ -290,18 +189,17 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
290 189  
291 190 END
292 191  
293   - 'POLSED': BEGIN ;UPDATED
294   -
295   - ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
  192 + 'POLSED': BEGIN
296 193  
297 194 IF !run_pol and !run_lin THEN BEGIN
298 195  
299   - dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec,dustem_polfrac);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
300   -
  196 + dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec,dustem_polfrac)
  197 +
  198 +
301 199 ENDIF
302 200 END
303 201  
304   - 'PSI_EM' : BEGIN ;This block is present if PSI_EM data is to be shown but not the rest of the stokes parameters ?
  202 + 'PSI_EM' : BEGIN
305 203  
306 204 IF !run_pol and !run_lin THEN BEGIN
307 205 degtorad = !pi/180
... ... @@ -310,58 +208,26 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
310 208 END
311 209  
312 210  
313   - 'POLFRAC': BEGIN ;We do not really need this block
314   -
315   - END
  211 +; 'POLFRAC': BEGIN ;We do not really need this block
  212 +;
  213 +; END
316 214  
317 215  
318   - 'QEXT': BEGIN ;NOT UPDATED ; Q - EXTINCTION CROSS SECTION
319   -
320   - chi2_qext =total(((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2)
  216 + 'QEXT': BEGIN ;UPDATED
  217 +
  218 + toto = dustem_compute_stokext(p_dim,st,dustem_qext,dustem_uext,QEXT_spec,UEXT_spec,PSIEXT_spec,dustem_psi_ext)
  219 + chi2_qext =total(((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2)
321 220 n_qext = n_elements((*!dustem_data.qext).values)
322 221 rchi2_qext = chi2_qext / n_qext
323 222 wrchi2_qext = rchi2_qext * !fit_rchi2_weight.qext
324 223 dustem_out = [ dustem_out, dustem_qext] ; classic command
325 224 message,"chi2 QEXT = "+strtrim(chi2_qext,2)+" rchi2 QEXT = "+strtrim(rchi2_qext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
326   - print,"chi2 per wl: ",((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2
327   - if !dustem_show_plot ne 0 then begin
328   -
329   - win+=1
330   - window ,win ,title='DUSTEM WRAP (QEXT)'
331   - titqex='Q-polarized Cross-section'
332   - ytitqex=textoidl('\sigma_{Q}/N_H (cm^2/H)')
333   - xtit=textoidl('\lambda^{-1} (\mum^{-1})')
334   - ;xr=[0.3,1e4]
335   - xr=[1e-4,3.5]
336   -
337   - yr=[1e-34,5e-24]
338   - indqex=where(Q_ext lt 0, countqex)
339   - ylog=1
340   - if countqex ne 0 then begin
341   - ylog=0
342   - yr=[-9e-27,5e-27]
343   - endif
344   - normpol = st.polext.ext_tot * 0. + 1.d21
345   - cgplot,1/st.polext.wav,Q_ext,xtit='',ytit=ytitqex,tit=titqex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
346   - cgoplot,1/st.polext.wav,Q_ext/normpol,color='black'
347   - cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,color='Indian Red',psym=16,symsize=1,thick=2
348   - err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,yrms=3.*((*!dustem_data.qext).sigma)/2/normpol,color=cgColor('Indian Red')
349   - cgoplot,1/(*!dustem_data.qext).wav,dustem_qext/normpol,psym=6,color='red',symsize=2
350   -
351   -
352   -
353   - cgplot,1/st.polsed.wav,Q_ext/Q_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
354   - cgoplot,1/st.polsed.wav,Q_ext/Q_ext,color='black'
355   - cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,color='Indian Red',psym=16,symsize=1,thick=2
356   - err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,yrms=3.*((*!dustem_data.qext).sigma)/2.,color=cgColor('Indian Red')
357   -
358   - endif
  225 + print,"chi2 per wl: ",((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2
359 226  
360 227 END
361 228 ;
362   - 'UEXT': BEGIN ;NOT UPDATED
363   -
364   -
  229 + 'UEXT': BEGIN ;UPDATED
  230 +
365 231 chi2_uext =total(((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2)
366 232 n_uext = n_elements((*!dustem_data.uext).values)
367 233 rchi2_uext = chi2_uext / n_uext
... ... @@ -369,48 +235,30 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
369 235 dustem_out = [ dustem_out, dustem_uext] ; classic command
370 236 message,"chi2 UEXT = "+strtrim(chi2_uext,2)+" rchi2 UEXT = "+strtrim(rchi2_uext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
371 237 print,"chi2 per wl: ",((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2
372   - if !dustem_show_plot ne 0 then begin
373   -
374   - win+=1
375   - window ,win ,title='DUSTEM WRAP (UEXT)'
376   - tituex='U-polarized Cross-section'
377   - ytituex=textoidl('\sigma_{U}/N_H (cm^2/H)')
378   - xtit=textoidl('\lambda^{-1} (\mum^{-1})')
379   - ;xr=[0.3,1e4]
380   - xr=[1e-4,3.5]
381   -
382   - yr=[1e-34,5e-24]
383   - induex=where(U_ext lt 0, countuex)
384   - ylog=1
385   - if countuex ne 0 then begin
386   - ylog=0
387   - yr=[-9e-27,5e-27]
388   - endif
389   - normpol = st.polext.ext_tot * 0. + 1.d21
390   - cgplot,1/st.polext.wav,U_ext,xtit='',ytit=ytituex,tit=tituex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
391   - cgoplot,1/st.polext.wav,U_ext/normpol,color='black'
392   - cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,color='Steel Blue',psym=16,symsize=1,thick=2
393   - err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
394   - cgoplot,1/(*!dustem_data.uext).wav,dustem_uext/normpol,psym=6,color='red',symsize=2
395   -
396   -
397   -
398   - cgplot,1/st.polsed.wav,U_ext/U_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
399   - cgoplot,1/st.polsed.wav,U_ext/U_ext,color='black'
400   - cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,color='Steel Blue',psym=16,symsize=1,thick=2
401   - err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
402   -
403   - endif
404   -
  238 +
405 239  
406 240 END
  241 +
  242 + 'PSI_EXT' : BEGIN ;This block is present if PSI_EM data is to be shown but not the rest of the stokes parameters ?
  243 +
  244 + IF !run_pol and !run_lin THEN BEGIN
  245 + degtorad = !pi/180
  246 + if ~isa(dustem_uext) or ~isa(dustem_qext) then toto = dustem_compute_stokext(p_dim,st,dustem_qext,dustem_uext,QEXT_spec,UEXT_spec,PSIEXT_spec,dustem_psi_ext)
  247 + ENDIF
  248 + END
407 249  
408   - ELSE :;IF tagnames(ind_data_tag(i_tag)) NE 'POLFRAC' THEN stop
  250 + ELSE : ;Do nothing
  251 +
409 252 ENDCASE
410 253  
411 254 ENDFOR
  255 +
  256 +
412 257 dustemwrap_plot,p_dim,st,dustem_sed,SED_spec,dustem_qsed,Q_spec,dustem_used,U_spec,dustem_polsed,P_spec,dustem_polfrac,SP_spec, dustem_psi_em, PSI_spec,_extra=_extra
413 258  
  259 +
  260 +
  261 +
414 262 DOF = n_elements(p_min)
415 263 total_chi2 = 0
416 264 total_rchi2 = 0
... ... @@ -463,6 +311,16 @@ if !run_pol then begin
463 311 + chi2_uext * !fit_rchi2_weight.uext
464 312 endif
465 313  
  314 +
  315 +
  316 +
  317 +stop
  318 +
  319 +;Using this expression for now for the weighting.
  320 +(*!dustem_fit).chi2 = total_wchi2
  321 +(*!dustem_fit).rchi2 = total_wchi2/(total_n-DOF)
  322 +
  323 +
466 324 message,'======================================================================',/continue
467 325 print,"Total weighted chi2 driving mpfitfun = ", total_wchi2
468 326 print,"Nb of data points = ", total_n
... ... @@ -473,7 +331,7 @@ if !run_pol then print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed $;
473 331 else print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed
474 332 print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF)
475 333 print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF)
476   -print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original
  334 +print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original ; I do not understand this line. Maybe redo some computation.
477 335 message,'======================================================================',/continue
478 336  
479 337 ;print,"Total reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac)/(n_sed+n_ext+n_polext+n_polsed+n_polfrac-DOF)
... ...