OV-GSO-DC / dustem-wrapper_idl

Commit 3ff9cac318ad787ffbfc10ca4aba383c1c9808ac

Authored by Ilyes Choubani
1 parent 6cde48ea
Exists in master

general update

Showing 13 changed files with 652 additions and 443 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_check_data.pro
  Diff comments   View file @3ff9cac
1   -FUNCTION dustem_check_data, sed=sed,ext=ext,fpol=fpol,ist=ist,pst=pst,fpst=fpst,qst=qst,ust=ust,exst=exst,pexst=pexst,qexst=qexst,uexst=uexst
2   -
3   -
4   -;THIS PROCEDURE CHECKS THE DATA PROVIDED BY THE USER AND RETURNS DIFFERENT STRUCTURES TO DEFINE (SET) THE COMPOSITE DATA FITTING STRUCTURE (!dustem_data)
5   -
6   -
7   -;write help section here
8   -;I don't think this function should return anything
9   -;This replaces a big portion of dustem_set_data, JP was right.
  1 +FUNCTION dustem_check_data,data_sed=data_sed,data_ext=data_ext,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,psi_em=psi_em,qsed=qsed,used=used,qext=qext,uext=uext,psi_ext=psi_ext
10 2  
11 3  
12 4 ;CHECKING THE FORMAT OF THE SED STRUCTURE
13   -IF KEYWORD_SET(sed) THEN BEGIN
  5 +IF KEYWORD_SET(data_sed) THEN BEGIN
  6 +
14 7 ans=''
15   - wavs=sed.wave
  8 + wavs=data_sed.wave
16 9 ;=== Impose central wavelengths for photometric channels
17   - ind=where(sed.filter NE 'SPECTRUM',count)
  10 + ind=where(data_sed.filter NE 'SPECTRUM',count)
18 11 IF count NE 0 THEN BEGIN
19   - wavs(ind)=dustem_filter2wav(sed(ind).filter)
  12 + wavs(ind)=dustem_filter2wav(data_sed(ind).filter)
20 13 ENDIF
21   - ;=== define observations structure
22   - obs_st={instru_names:sed.instru,filt_names:sed.filter,wav:wavs, values:sed.StokesI,sigma:sqrt(abs(sed.sigmaII))}
  14 + ;=== Initializing observation structure
  15 + ;it doesn't matter that data.StokesI is used, further test will take care of what is necessary.
  16 + obs_EM={instru_names:data_sed.instru,filt_names:data_sed.filter,wav:wavs, values:data_sed.StokesI,sigma:sqrt(abs(data_sed.sigmaII))}
23 17 ;INITIALIZING OUTPUT STRUCTURES
24   - ist = obs_st
25   - pst = ist
26   - fpst = ist
27   - qst = ist
28   - ust = ist
29   - ;This allows the user to fit polsed, qsed and used individually.
  18 + sed = obs_EM
  19 + polsed = sed
  20 + polfrac = sed
  21 + qsed = sed
  22 + used = sed
  23 + psi_em = sed
  24 +
30 25 ;Test if there is StokesI data. If not move to polsed testing.
31   - ind = where(sed.StokesI NE la_undef() and finite(sed.StokesI) EQ 1 and sed.sigmaII NE la_undef() and finite(sed.sigmaII) EQ 1,N_sed) ;Testing if the data and errors are set
  26 + ind = where(data_sed.StokesI NE la_undef() and finite(data_sed.StokesI) EQ 1 and data_sed.sigmaII NE la_undef() and finite(data_sed.sigmaII) EQ 1,N_sed) ;Testing if the data and errors are set
32 27 IF N_sed EQ 0 THEN BEGIN
33   - ist = ptr_new()
  28 + sed = ptr_new()
34 29 goto, pzone
35 30 ENDIF ELSE BEGIN
36   - ind = where(finite(sed.StokesI) EQ 0 or finite(sed.sigmaII) EQ 0, countnoki) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  31 + ind = where(finite(data_sed.StokesI) EQ 0 or finite(data_sed.sigmaII) EQ 0, countnoki) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
37 32 IF countnoki NE 0 THEN BEGIN
38 33 message, "SED format isn't met. NaN value(s) detected in StokesI data.",/continue
39 34 message, 'Without automatic format fitering you cannot proceed.',/continue
... ... @@ -45,27 +40,29 @@ IF KEYWORD_SET(sed) THEN BEGIN
45 40 ENDIF
46 41 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
47 42 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
48   - testi=where(sed.sigmaII NE 0. and finite(sed.sigmaII) EQ 1 and sed.sigmaII NE la_undef() and sed.StokesI NE la_undef() and finite(sed.StokesI) EQ 1,ctesti)
  43 + tdaerri=where((data_sed.sigmaII EQ la_undef() and data_sed.StokesI NE la_undef()) or (data_sed.sigmaII NE la_undef() and data_sed.StokesI EQ la_undef()),ctdaerri)
  44 + IF ctdaerri NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  45 + testi=where(data_sed.sigmaII NE 0. and finite(data_sed.sigmaII) EQ 1 and data_sed.sigmaII NE la_undef() and data_sed.StokesI NE la_undef() and finite(data_sed.StokesI) EQ 1,ctesti)
49 46 IF ctesti NE 0 THEN BEGIN
50   - instru_names = ist.instru_names(testi)
51   - filt_names = ist.filt_names(testi)
52   - wav = ist.wav(testi)
53   - values = ist.values(testi)
54   - sigma = ist.sigma(testi)
55   - ist={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  47 + instru_names = sed.instru_names(testi)
  48 + filt_names = sed.filt_names(testi)
  49 + wav = sed.wav(testi)
  50 + values = sed.values(testi)
  51 + sigma = sed.sigma(testi)
  52 + sed={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
56 53 ENDIF
57 54 ENDELSE
58 55 pzone:
59 56 ans=''
60   - pst.values=sed.largeP
61   - pst.sigma=sqrt(abs(sed.sigma_largeP))
62   - ;Testig if there is LargeP data. If not move to polfrac data (fpol)
63   - ind = where(sed.largeP NE la_undef() and finite(sed.largeP) EQ 1 and sed.sigma_largeP NE la_undef() and finite(sed.sigma_largeP) EQ 1,NP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta)
  57 + polsed.values=data_sed.largeP
  58 + polsed.sigma=sqrt(abs(data_sed.sigma_largeP))
  59 + ;Testing if there is LargeP data. If not move to polfrac data.
  60 + ind = where(data_sed.largeP NE la_undef() and finite(data_sed.largeP) EQ 1 and data_sed.sigma_largeP NE la_undef() and finite(data_sed.sigma_largeP) EQ 1,NP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta)
64 61 IF NP EQ 0 THEN BEGIN
65   - pst = ptr_new()
66   - goto, fpzone
  62 + polsed = ptr_new()
  63 + goto, psiemzone
67 64 ENDIF ELSE BEGIN
68   - ind = where(finite(sed.largeP) EQ 0 or finite(sed.sigma_largeP) EQ 0, countnokp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  65 + ind = where(finite(data_sed.largeP) EQ 0 or finite(data_sed.sigma_largeP) EQ 0, countnokp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
69 66 IF countnokp NE 0 THEN BEGIN
70 67 message, "SED format isn't met. NaN value(s) detected in LargeP data.",/continue
71 68 message, 'Without automatic format fitering you cannot proceed.',/continue
... ... @@ -77,28 +74,64 @@ IF KEYWORD_SET(sed) THEN BEGIN
77 74 ENDIF
78 75 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
79 76 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
80   - testp=where(sed.sigma_largeP NE 0. and finite(sed.sigma_largeP) EQ 1 and sed.sigma_largeP NE la_undef() and sed.largeP NE la_undef() and finite(sed.largeP) EQ 1,ctestp)
  77 + tdaerrp=where((data_sed.sigma_largeP EQ la_undef() and data_sed.largeP NE la_undef()) or (data_sed.sigma_largeP NE la_undef() and data_sed.largeP EQ la_undef()),ctdaerrp)
  78 + IF ctdaerrp NE 0 THEN message, 'Data and errors do not match. Undefinded value(s) detected. Aborting...'
  79 + testp=where(data_sed.sigma_largeP NE 0. and finite(data_sed.sigma_largeP) EQ 1 and data_sed.sigma_largeP NE la_undef() and data_sed.largeP NE la_undef() and finite(data_sed.largeP) EQ 1,ctestp)
81 80 IF ctestp NE 0 THEN BEGIN
82   - instru_names = pst.instru_names(testp)
83   - filt_names = pst.filt_names(testp)
84   - wav = pst.wav(testp)
85   - values = pst.values(testp)
86   - sigma = pst.sigma(testp)
87   - pst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  81 + instru_names = polsed.instru_names(testp)
  82 + filt_names = polsed.filt_names(testp)
  83 + wav = polsed.wav(testp)
  84 + values = polsed.values(testp)
  85 + sigma = polsed.sigma(testp)
  86 + polsed={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
88 87 ENDIF
89   - ENDELSE
  88 + ENDELSE ; The deprecation of the fitting of POLFRAC does not happen here. But rather in dustem_mpfit_run
  89 + psiemzone:
  90 + ans=''
  91 + psi_em.values = data_sed.psi
  92 + psi_em.sigma = sqrt(abs(data_sed.sigma_psi))
  93 + ;Testing if there is PSI (EMISSION) data. If not move to polfrac data.
  94 + ind = where(data_sed.psi NE la_undef() and finite(data_sed.psi) EQ 1 and data_sed.sigma_psi NE la_undef() and finite(data_sed.sigma_psi) EQ 1,NPSIEM) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta)
  95 + IF NPSIEM EQ 0 THEN BEGIN
  96 + psi_em = ptr_new()
  97 + goto, fpzone
  98 + ENDIF ELSE BEGIN
  99 + ind = where(finite(data_sed.psi) EQ 0 or finite(data_sed.sigma_psi) EQ 0, countnokpsiem) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  100 + IF countnokpsiem NE 0 THEN BEGIN
  101 + message, "SED format isn't met. NaN value(s) detected in LargeP data.",/continue
  102 + message, 'Without automatic format fitering you cannot proceed.',/continue
  103 + read, ans, prompt='Would you like automatic format filtering? (Y/N):'
  104 + IF strupcase(ans) EQ 'N' THEN BEGIN
  105 + message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
  106 + message, 'Aborting ...'
  107 + ENDIF
  108 + ENDIF
  109 + IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
  110 + ELSE message, 'Data format is met. Removal of undefined data points...',/continue
  111 + tdaerrpsiem=where((data_sed.sigma_psi EQ la_undef() and data_sed.psi NE la_undef()) or (data_sed.sigma_psi NE la_undef() and data_sed.psi EQ la_undef()),ctdaerrpsiem)
  112 + IF ctdaerrpsiem NE 0 THEN message, 'Data and errors do not match. Undefinded value(s) detected. Aborting...'
  113 + testpsiem=where(data_sed.sigma_psi NE 0. and finite(data_sed.sigma_psi) EQ 1 and data_sed.sigma_psi NE la_undef() and data_sed.psi NE la_undef() and finite(data_sed.psi) EQ 1,ctestpsiem)
  114 + IF ctestpsiem NE 0 THEN BEGIN
  115 + instru_names = psi_em.instru_names(testpsiem)
  116 + filt_names = psi_em.filt_names(testpsiem)
  117 + wav = psi_em.wav(testpsiem)
  118 + values = psi_em.values(testpsiem)
  119 + sigma = psi_em.sigma(testpsiem)
  120 + psi_em={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  121 + ENDIF
  122 + ENDELSE ; The deprecation of POLFRAC does not happen here. More like in dustem_mpfit_run
90 123 fpzone:
91 124 ans=''
92   - fpst.values=sed.smallp
93   - fpst.sigma=sqrt(abs(sed.sigma_smallp))
  125 + polfrac.values=data_sed.smallp
  126 + polfrac.sigma=sqrt(abs(data_sed.sigma_smallp))
94 127 ;Testig if there is smallp data. If not move to qsed data
95   - ind = where(sed.smallP NE la_undef() and finite(sed.smallP) EQ 1 and sed.sigma_smallP NE la_undef() and finite(sed.sigma_smallP) EQ 1,NSP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta)
  128 + ind = where(data_sed.smallP NE la_undef() and finite(data_sed.smallP) EQ 1 and data_sed.sigma_smallP NE la_undef() and finite(data_sed.sigma_smallP) EQ 1,NSP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta)
96 129 IF NSP EQ 0 THEN BEGIN
97   - fpst = ptr_new()
  130 + polfrac = ptr_new()
98 131 ;goto, qzone
99 132 goto, qzone
100 133 ENDIF ELSE BEGIN
101   - ind = where(finite(sed.smallP) EQ 0 or finite(sed.sigma_smallP) EQ 0, countnoksp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  134 + ind = where(finite(data_sed.smallP) EQ 0 or finite(data_sed.sigma_smallP) EQ 0, countnoksp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
102 135 IF countnoksp NE 0 THEN BEGIN
103 136 message, "SED format isn't met. NaN value(s) detected in SmallP data.",/continue
104 137 message, 'Without automatic format fitering you cannot proceed.',/continue
... ... @@ -110,27 +143,29 @@ IF KEYWORD_SET(sed) THEN BEGIN
110 143 ENDIF
111 144 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
112 145 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
113   - testsp=where(sed.sigma_smallP NE 0. and finite(sed.sigma_smallP) EQ 1 and sed.sigma_smallP NE la_undef() and sed.smallP NE la_undef() and finite(sed.smallP) EQ 1,ctestsp)
  146 + tdaerrsp=where((data_sed.sigma_smallp EQ la_undef() and data_sed.smallP NE la_undef()) or (data_sed.sigma_smallp NE la_undef() and data_sed.smallP EQ la_undef()),ctdaerrsp)
  147 + IF ctdaerrsp NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  148 + testsp=where(data_sed.sigma_smallP NE 0. and finite(data_sed.sigma_smallP) EQ 1 and data_sed.sigma_smallP NE la_undef() and data_sed.smallP NE la_undef() and finite(data_sed.smallP) EQ 1,ctestsp)
114 149 IF ctestsp NE 0 THEN BEGIN
115   - instru_names = fpst.instru_names(testsp)
116   - filt_names = fpst.filt_names(testsp)
117   - wav = fpst.wav(testsp)
118   - values = fpst.values(testsp)
119   - sigma = fpst.sigma(testsp)
120   - fpst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  150 + instru_names = polfrac.instru_names(testsp)
  151 + filt_names = polfrac.filt_names(testsp)
  152 + wav = polfrac.wav(testsp)
  153 + values = polfrac.values(testsp)
  154 + sigma = polfrac.sigma(testsp)
  155 + polfrac={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
121 156 ENDIF
122 157 ENDELSE
123 158 qzone:
124 159 ans=''
125   - qst.values=sed.StokesQ
126   - qst.sigma=sqrt(abs(sed.sigmaQQ))
  160 + qsed.values=data_sed.StokesQ
  161 + qsed.sigma=sqrt(abs(data_sed.sigmaQQ))
127 162 ;Testig if there is StokesQ data. If not move to used data
128   - ind = where(sed.StokesQ NE la_undef() and finite(sed.StokesQ) EQ 1 and sed.sigmaQQ NE la_undef() and finite(sed.sigmaQQ) EQ 1,NQ)
  163 + ind = where(data_sed.StokesQ NE la_undef() and finite(data_sed.StokesQ) EQ 1 and data_sed.sigmaQQ NE la_undef() and finite(data_sed.sigmaQQ) EQ 1,NQ)
129 164 IF NQ EQ 0 THEN BEGIN
130   - qst = ptr_new()
  165 + qsed = ptr_new()
131 166 goto, uzone
132 167 ENDIF ELSE BEGIN
133   - ind = where(finite(sed.StokesQ) EQ 0 or finite(sed.sigmaQQ) EQ 0, countnokq) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  168 + ind = where(finite(data_sed.StokesQ) EQ 0 or finite(data_sed.sigmaQQ) EQ 0, countnokq) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
134 169 IF countnokq NE 0 THEN BEGIN
135 170 message, "SED format isn't met. NaN value(s) detected in StokesQ data.",/continue
136 171 message, 'Without automatic format fitering you cannot proceed.',/continue
... ... @@ -142,27 +177,29 @@ IF KEYWORD_SET(sed) THEN BEGIN
142 177 ENDIF
143 178 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
144 179 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
145   - testq=where(sed.sigmaQQ NE 0. and finite(sed.sigmaQQ) EQ 1 and sed.sigmaQQ NE la_undef() and sed.StokesQ NE la_undef() and finite(sed.StokesQ) EQ 1,ctestq)
  180 + tdaerrq=where((data_sed.sigmaQQ EQ la_undef() and data_sed.StokesQ NE la_undef()) or (data_sed.sigmaQQ NE la_undef() and data_sed.StokesQ EQ la_undef()),ctdaerrq)
  181 + IF ctdaerrq NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  182 + testq=where(data_sed.sigmaQQ NE 0. and finite(data_sed.sigmaQQ) EQ 1 and data_sed.sigmaQQ NE la_undef() and data_sed.StokesQ NE la_undef() and finite(data_sed.StokesQ) EQ 1,ctestq)
146 183 IF ctestq NE 0 THEN BEGIN
147   - instru_names = qst.instru_names(testq)
148   - filt_names = qst.filt_names(testq)
149   - wav = qst.wav(testq)
150   - values = qst.values(testq)
151   - sigma = qst.sigma(testq)
152   - qst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  184 + instru_names = qsed.instru_names(testq)
  185 + filt_names = qsed.filt_names(testq)
  186 + wav = qsed.wav(testq)
  187 + values = qsed.values(testq)
  188 + sigma = qsed.sigma(testq)
  189 + qsed={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
153 190 ENDIF
154 191 ENDELSE
155 192 uzone:
156 193 ans=''
157   - ust.values=sed.StokesU
158   - ust.sigma=sqrt(abs(sed.sigmaUU))
  194 + used.values=data_sed.StokesU
  195 + used.sigma=sqrt(abs(data_sed.sigmaUU))
159 196 ;Testig if there is StokesU data. If not return empty value and end loop.
160   - ind = where(sed.StokesU NE la_undef() and finite(sed.StokesU) EQ 1 and sed.sigmaUU NE la_undef() and finite(sed.sigmaUU) EQ 1,NU)
  197 + ind = where(data_sed.StokesU NE la_undef() and finite(data_sed.StokesU) EQ 1 and data_sed.sigmaUU NE la_undef() and finite(data_sed.sigmaUU) EQ 1,NU)
161 198 IF NU EQ 0 THEN BEGIN
162   - ust = ptr_new()
  199 + used = ptr_new()
163 200 goto, end_sed
164 201 ENDIF ELSE BEGIN
165   - ind = where(finite(sed.StokesU) EQ 0 or finite(sed.sigmaUU) EQ 0, countnoku) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  202 + ind = where(finite(data_sed.StokesU) EQ 0 or finite(data_sed.sigmaUU) EQ 0, countnoku) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
166 203 IF countnoku NE 0 THEN BEGIN
167 204 message, "SED format isn't met. NaN value(s) detected in StokesU data.",/continue
168 205 message, 'Without automatic format fitering you cannot proceed.',/continue
... ... @@ -174,46 +211,46 @@ IF KEYWORD_SET(sed) THEN BEGIN
174 211 ENDIF
175 212 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
176 213 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
177   - testu=where(sed.sigmaUU NE 0. and finite(sed.sigmaUU) EQ 1 and sed.sigmaUU NE la_undef() and sed.StokesU NE la_undef() and finite(sed.StokesU) EQ 1,ctestu)
  214 + tdaerru=where((data_sed.sigmaUU EQ la_undef() and data_sed.StokesU NE la_undef()) or (data_sed.sigmaUU NE la_undef() and data_sed.StokesU EQ la_undef()),ctdaerru)
  215 + IF ctdaerru NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  216 + testu=where(data_sed.sigmaUU NE 0. and finite(data_sed.sigmaUU) EQ 1 and data_sed.sigmaUU NE la_undef() and data_sed.StokesU NE la_undef() and finite(data_sed.StokesU) EQ 1,ctestu)
178 217 IF ctestu NE 0 THEN BEGIN
179   - instru_names = ust.instru_names(testu)
180   - filt_names = ust.filt_names(testu)
181   - wav = ust.wav(testu)
182   - values = ust.values(testu)
183   - sigma = ust.sigma(testu)
184   - ust={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  218 + instru_names = used.instru_names(testu)
  219 + filt_names = used.filt_names(testu)
  220 + wav = used.wav(testu)
  221 + values = used.values(testu)
  222 + sigma = used.sigma(testu)
  223 + used={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
185 224 ENDIF
186 225 ENDELSE
187 226 end_sed:
188 227 ENDIF
189 228  
190 229  
191   -
192   -;CHECKING THE FORMAT OF THE EXT STRUCTURE
193   -IF KEYWORD_SET(ext) THEN begin
  230 +IF KEYWORD_SET(data_ext) THEN begin
  231 +
194 232 ans=''
195   - wavs=ext.wave
  233 + wavs=data_ext.wave
196 234 ;=== Impose central wavelengths for photometric channels
197   - ind=where(ext.filter NE 'SPECTRUM',count)
  235 + ind=where(data_ext.filter NE 'SPECTRUM',count)
198 236 IF count NE 0 THEN BEGIN
199   - wavs(ind)=dustem_filter2wav(ext(ind).filter)
  237 + wavs(ind)=dustem_filter2wav(data_ext(ind).filter)
200 238 ENDIF
201 239 ;=== define observations structure
202   - obs_st={instru_names:ext.instru,filt_names:ext.filter,wav:wavs, values:ext.EXT_I,sigma:sqrt(abs(ext.sigextII))}
  240 + obs_EXT={instru_names:data_ext.instru,filt_names:data_ext.filter,wav:wavs, values:data_ext.EXT_I,sigma:sqrt(abs(data_ext.sigextII))}
203 241 ;INITIALIZING OUTPUT STRUCTURES
204   - exst = obs_st
205   - pexst = exst
206   - qexst = exst
207   - uexst = exst
208   - ;This allows the user to fit polext, qext and uext individually.
209   - ;Test if there is extinction data. If not move to polext testing.
210   - ;this condition was updated because you cannot have data without errors and errors without data
211   - ind = where(ext.EXT_I NE la_undef() and finite(ext.EXT_I) EQ 1 and ext.sigextII NE la_undef() and finite(ext.sigextII) EQ 1,N_ext) ;Testing if the data and errors are set
  242 + ext = obs_EXT
  243 + polext = ext
  244 + qext = ext
  245 + uext = ext
  246 + psi_ext = ext
  247 +
  248 + ind = where(data_ext.EXT_I NE la_undef() and finite(data_ext.EXT_I) EQ 1 and data_ext.sigextII NE la_undef() and finite(data_ext.sigextII) EQ 1,N_ext) ;Testing if the data and errors are set
212 249 IF N_ext EQ 0 THEN BEGIN
213   - exst = ptr_new()
  250 + ext = ptr_new()
214 251 goto, pexzone
215 252 ENDIF ELSE BEGIN
216   - ind = where(finite(ext.EXT_I) EQ 0 or finite(ext.sigextII) EQ 0, countnokex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  253 + ind = where(finite(data_ext.EXT_I) EQ 0 or finite(data_ext.sigextII) EQ 0, countnokex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
217 254 IF countnokex NE 0 THEN BEGIN
218 255 message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue
219 256 message, 'Without automatic format fitering you cannot proceed.',/continue
... ... @@ -226,29 +263,31 @@ IF KEYWORD_SET(ext) THEN begin
226 263 ENDIF
227 264 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
228 265 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
229   - testex=where(ext.sigextII NE 0. and finite(ext.sigextII) EQ 1 and ext.sigextII NE la_undef() and ext.EXT_I NE la_undef() and finite(ext.EXT_I) EQ 1,ctestex)
  266 + tdaerrex=where((data_ext.sigextII EQ la_undef() and data_ext.EXT_I NE la_undef()) or (data_ext.sigextII NE la_undef() and data_ext.EXT_I EQ la_undef()),ctdaerrex)
  267 + IF ctdaerrex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  268 + testex=where(data_ext.sigextII NE 0. and finite(data_ext.sigextII) EQ 1 and data_ext.sigextII NE la_undef() and data_ext.EXT_I NE la_undef() and finite(data_ext.EXT_I) EQ 1,ctestex)
230 269 IF ctestex NE 0 THEN BEGIN
231   - instru_names = exst.instru_names(testex)
232   - filt_names = exst.filt_names(testex)
233   - wav = exst.wav(testex)
234   - values = exst.values(testex)
235   - sigma = exst.sigma(testex)
236   - exst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  270 + instru_names = ext.instru_names(testex)
  271 + filt_names = ext.filt_names(testex)
  272 + wav = ext.wav(testex)
  273 + values = ext.values(testex)
  274 + sigma = ext.sigma(testex)
  275 + ext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
237 276 ENDIF
238 277 ENDELSE
239 278 pexzone:
240 279 ans=''
241   - pexst.values=ext.EXT_P
242   - pexst.sigma=sqrt(abs(ext.sigextP))
  280 + polext.values=data_ext.EXT_P
  281 + polext.sigma=sqrt(abs(data_ext.sigextP))
243 282 ;Testig if there is EXT_P data. If not move to EXT_Q testing
244   - ind = where(ext.EXT_P NE la_undef() and finite(ext.EXT_P) EQ 1 and ext.sigextP NE la_undef() and finite(ext.sigextII) EQ 1,NP_ext) ;Testing if the data and errors are set
  283 + ind = where(data_ext.EXT_P NE la_undef() and finite(data_ext.EXT_P) EQ 1 and data_ext.sigextP NE la_undef() and finite(data_ext.sigextP) EQ 1,NP_ext) ;Testing if the data and errors are set
245 284 IF NP_ext EQ 0 THEN BEGIN
246   - pexst = ptr_new()
  285 + polext = ptr_new()
247 286 goto, qexzone
248 287 ENDIF ELSE BEGIN
249   - ind = where(finite(ext.EXT_P) EQ 0 or finite(ext.sigextP) EQ 0, countnokpex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  288 + ind = where(finite(data_ext.EXT_P) EQ 0 or finite(data_ext.sigextP) EQ 0, countnokpex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
250 289 IF countnokpex NE 0 THEN BEGIN
251   - message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue
  290 + message, "SED format isn't met. NaN value(s) detected in EXT_P data.",/continue
252 291 message, 'Without automatic format fitering you cannot proceed.',/continue
253 292 ans=''
254 293 read, ans, prompt='Would you like automatic format filtering? (Y/N):'
... ... @@ -259,30 +298,67 @@ IF KEYWORD_SET(ext) THEN begin
259 298 ENDIF
260 299 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
261 300 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
262   - testpex=where(ext.sigextP NE 0. and finite(ext.sigextP) EQ 1 and ext.sigextP NE la_undef() and ext.EXT_P NE la_undef() and finite(ext.EXT_P) EQ 1,ctestpex)
  301 + tdaerrpex=where((data_ext.sigextP EQ la_undef() and data_ext.EXT_P NE la_undef()) or (data_ext.sigextP NE la_undef() and data_ext.EXT_P EQ la_undef()),ctdaerrpex)
  302 + IF ctdaerrpex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  303 + testpex=where(data_ext.sigextP NE 0. and finite(data_ext.sigextP) EQ 1 and data_ext.sigextP NE la_undef() and data_ext.EXT_P NE la_undef() and finite(data_ext.EXT_P) EQ 1,ctestpex)
263 304 IF ctestpex NE 0 THEN BEGIN
264   - instru_names = pexst.instru_names(testpex)
265   - filt_names = pexst.filt_names(testpex)
266   - wav = pexst.wav(testpex)
267   - values = pexst.values(testpex)
268   - sigma = pexst.sigma(testpex)
269   - pexst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  305 + instru_names = polext.instru_names(testpex)
  306 + filt_names = polext.filt_names(testpex)
  307 + wav = polext.wav(testpex)
  308 + values = polext.values(testpex)
  309 + sigma = polext.sigma(testpex)
  310 + polext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  311 + ENDIF
  312 + ENDELSE
  313 + psiexzone:
  314 + ans=''
  315 + psi_ext.values=data_ext.psi
  316 + psi_ext.sigma=sqrt(abs(data_ext.sigextpsi))
  317 + ;Testig if there is EXT_P data. If not move to EXT_Q testing
  318 + ind = where(data_ext.psi NE la_undef() and finite(data_ext.psi) EQ 1 and data_ext.sigextpsi NE la_undef() and finite(data_ext.sigextpsi) EQ 1,NPSIEXT) ;Testing if the data and errors are set
  319 + IF NPSIEXT EQ 0 THEN BEGIN
  320 + psi_ext = ptr_new()
  321 + goto, qexzone
  322 + ENDIF ELSE BEGIN
  323 + ind = where(finite(data_ext.psi) EQ 0 or finite(data_ext.sigextpsi) EQ 0, countnokpsiex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  324 + IF countnokpsiex NE 0 THEN BEGIN
  325 + message, "SED format isn't met. NaN value(s) detected in PSI_EXT data.",/continue
  326 + message, 'Without automatic format fitering you cannot proceed.',/continue
  327 + ans=''
  328 + read, ans, prompt='Would you like automatic format filtering? (Y/N):'
  329 + IF strupcase(ans) EQ 'N' THEN BEGIN
  330 + message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
  331 + message, 'Aborting ...'
  332 + ENDIF
  333 + ENDIF
  334 + IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
  335 + ELSE message, 'Data format is met. Removal of undefined data points...',/continue
  336 + tdaerrpsiex=where((data_ext.sigextpsi EQ la_undef() and data_ext.psi NE la_undef()) or (data_ext.sigextpsi NE la_undef() and data_ext.psi EQ la_undef()),ctdaerrpsiex)
  337 + IF ctdaerrpsiex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  338 + testpsiex=where(data_ext.sigextpsi NE 0. and finite(data_ext.sigextpsi) EQ 1 and data_ext.sigextpsi NE la_undef() and data_ext.psi NE la_undef() and finite(data_ext.psi) EQ 1,ctestpsiex)
  339 + IF ctestpsiex NE 0 THEN BEGIN
  340 + instru_names = psi_ext.instru_names(testpsiex)
  341 + filt_names = psi_ext.filt_names(testpsiex)
  342 + wav = psi_ext.wav(testpsiex)
  343 + values = psi_ext.values(testpsiex)
  344 + sigma = psi_ext.sigma(testpsiex)
  345 + psi_ext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
270 346 ENDIF
271 347 ENDELSE
272 348 qexzone:
273 349 ans=''
274   - qexst.values=ext.EXT_Q
275   - qexst.sigma=sqrt(abs(ext.sigextQQ))
  350 + qext.values=data_ext.EXT_Q
  351 + qext.sigma=sqrt(abs(data_ext.sigextQQ))
276 352 ;Testig if there is EXT_P data. If not move to EXT_U testing
277   - ind = where(ext.EXT_Q NE la_undef() and finite(ext.EXT_Q) EQ 1 and ext.sigextQQ NE la_undef() and finite(ext.sigextQQ) EQ 1,NQ_ext) ;Testing if the data and errors are set
  353 + ind = where(data_ext.EXT_Q NE la_undef() and finite(data_ext.EXT_Q) EQ 1 and data_ext.sigextQQ NE la_undef() and finite(data_ext.sigextQQ) EQ 1,NQ_ext) ;Testing if the data and errors are set
278 354 ;dustem_set_data.prostop
279 355 IF NQ_ext EQ 0 THEN BEGIN
280   - qexst = ptr_new()
  356 + qext = ptr_new()
281 357 goto, uexzone
282 358 ENDIF ELSE BEGIN
283   - ind = where(finite(ext.EXT_Q) EQ 0 or finite(ext.sigextQQ) EQ 0, countnokqex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  359 + ind = where(finite(data_ext.EXT_Q) EQ 0 or finite(data_ext.sigextQQ) EQ 0, countnokqex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
284 360 IF countnokqex NE 0 THEN BEGIN
285   - message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue
  361 + message, "SED format isn't met. NaN value(s) detected in EXT_Q data.",/continue
286 362 message, 'Without automatic format fitering you cannot proceed.',/continue
287 363 ans=''
288 364 read, ans, prompt='Would you like automatic format filtering? (Y/N):'
... ... @@ -293,29 +369,31 @@ IF KEYWORD_SET(ext) THEN begin
293 369 ENDIF
294 370 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
295 371 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
296   - testqex=where(ext.sigextQQ NE 0. and finite(ext.sigextQQ) EQ 1 and ext.sigextQQ NE la_undef() and ext.EXT_Q NE la_undef() and finite(ext.EXT_Q) EQ 1,ctestqex)
  372 + tdaerrqex=where((data_ext.sigextQQ EQ la_undef() and data_ext.EXT_Q NE la_undef()) or (data_ext.sigextQQ NE la_undef() and data_ext.EXT_Q EQ la_undef()),ctdaerrqex)
  373 + IF ctdaerrqex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  374 + testqex=where(data_ext.sigextQQ NE 0. and finite(data_ext.sigextQQ) EQ 1 and data_ext.sigextQQ NE la_undef() and data_ext.EXT_Q NE la_undef() and finite(data_ext.EXT_Q) EQ 1,ctestqex)
297 375 IF ctestqex NE 0 THEN BEGIN
298   - instru_names = qexst.instru_names(testqex)
299   - filt_names = qexst.filt_names(testqex)
300   - wav = qexst.wav(testqex)
301   - values = qexst.values(testqex)
302   - sigma = qexst.sigma(testqex)
303   - qexst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  376 + instru_names = qext.instru_names(testqex)
  377 + filt_names = qext.filt_names(testqex)
  378 + wav = qext.wav(testqex)
  379 + values = qext.values(testqex)
  380 + sigma = qext.sigma(testqex)
  381 + qext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
304 382 ENDIF
305 383 ENDELSE
306 384 uexzone:
307 385 ans=''
308   - uexst.values=ext.EXT_U
309   - uexst.sigma=sqrt(abs(ext.sigextUU))
  386 + uext.values=data_ext.EXT_U
  387 + uext.sigma=sqrt(abs(data_ext.sigextUU))
310 388 ;Testig if there is EXT_P data. If not return empty value and end loop.
311   - ind = where(ext.EXT_U NE la_undef() and finite(ext.EXT_U) EQ 1 and ext.sigextUU NE la_undef() and finite(ext.sigextUU) EQ 1,NU_ext) ;Testing if the data and errors are set
  389 + ind = where(data_ext.EXT_U NE la_undef() and finite(data_ext.EXT_U) EQ 1 and data_ext.sigextUU NE la_undef() and finite(data_ext.sigextUU) EQ 1,NU_ext) ;Testing if the data and errors are set
312 390 IF NU_ext EQ 0 THEN BEGIN
313   - uexst = ptr_new()
  391 + uext = ptr_new()
314 392 goto, end_ext
315 393 ENDIF ELSE BEGIN
316   - ind = where(finite(ext.EXT_U) EQ 0 or finite(ext.sigextUU) EQ 0, countnokuex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
  394 + ind = where(finite(data_ext.EXT_U) EQ 0 or finite(data_ext.sigextUU) EQ 0, countnokuex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
317 395 IF countnokuex NE 0 THEN BEGIN
318   - message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue
  396 + message, "SED format isn't met. NaN value(s) detected in EXT_U data.",/continue
319 397 message, 'Without automatic format fitering you cannot proceed.',/continue
320 398 ans=''
321 399 read, ans, prompt='Would you like automatic format filtering? (Y/N):'
... ... @@ -326,21 +404,22 @@ IF KEYWORD_SET(ext) THEN begin
326 404 ENDIF
327 405 IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
328 406 ELSE message, 'Data format is met. Removal of undefined data points...',/continue
329   - testuex=where(ext.sigextUU NE 0. and finite(ext.sigextUU) EQ 1 and ext.sigextUU NE la_undef() and ext.EXT_U NE la_undef() and finite(ext.EXT_U) EQ 1,ctestuex)
  407 + tdaerruex=where((data_ext.sigextUU EQ la_undef() and data_ext.EXT_U NE la_undef()) or (data_ext.sigextUU NE la_undef() and data_ext.EXT_U EQ la_undef()),ctdaerruex)
  408 + IF ctdaerruex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
  409 + testuex=where(data_ext.sigextUU NE 0. and finite(data_ext.sigextUU) EQ 1 and data_ext.sigextUU NE la_undef() and data_ext.EXT_U NE la_undef() and finite(data_ext.EXT_U) EQ 1,ctestuex)
330 410 IF ctestuex NE 0 THEN BEGIN
331   - instru_names = uexst.instru_names(testuex)
332   - filt_names = uexst.filt_names(testuex)
333   - wav = uexst.wav(testuex)
334   - values = uexst.values(testuex)
335   - sigma = uexst.sigma(testuex)
336   - uexst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
  411 + instru_names = uext.instru_names(testuex)
  412 + filt_names = uext.filt_names(testuex)
  413 + wav = uext.wav(testuex)
  414 + values = uext.values(testuex)
  415 + sigma = uext.sigma(testuex)
  416 + uext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
337 417 ENDIF
338 418 ENDELSE
339 419 end_ext:
340 420 ENDIF
341   -
342   -
343   -RETURN, 0.
344 421 ;stop
  422 +RETURN, 0.
  423 +
345 424  
346 425 END
... ...
src/idl/dustem_compute_polext.pro
  Diff comments   View file @3ff9cac
... ... @@ -10,7 +10,6 @@ IF not keyword_set(sti) THEN BEGIN
10 10 sti=dustem_run(p_dim)
11 11 ENDIF
12 12  
13   -
14 13 If ptr_valid(*!dustem_plugin) THEN BEGIN
15 14 ;ADDING PLUGIN TO SPECTRUM----------------
16 15 ;if n_tags(!dustem_data.polext) gt 1 then begin
... ... @@ -20,7 +19,6 @@ If ptr_valid(*!dustem_plugin) THEN BEGIN
20 19 ;endif
21 20 ;------------------------------------------
22 21  
23   -
24 22 ; ;=============================THIS BLOCK IS WRONG YOU HAVE TO CHANGE IT=============================
25 23 ;
26 24 ;
... ... @@ -51,8 +49,6 @@ If ptr_valid(*!dustem_plugin) THEN BEGIN
51 49  
52 50 ENDIF
53 51  
54   -
55   -
56 52 dustem_polext = interpol(sti.polext.ext_tot,sti.polext.wav,(*!dustem_data.polext).wav)
57 53  
58 54  
... ...
src/idl/dustem_compute_polsed.pro
  Diff comments   View file @3ff9cac
... ... @@ -59,6 +59,9 @@ frac=P/stI
59 59 tes=where(finite(frac) eq 0)
60 60 frac(tes)=0.
61 61  
  62 +;THIS LOOP IS WRONG AND DOES NOT MAKE SENSE
  63 +
  64 +
62 65 ;IF ptr_valid(!dustem_plugin) THEN BEGIN
63 66 scopes=tag_names((*!dustem_scope))
64 67 IF scopes[0] NE 'NONE' THEN BEGIN
... ... @@ -77,7 +80,7 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
77 80 ENDFOR
78 81  
79 82  
80   -IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0.,Nwaves)
  83 +IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(!dustem_psi,Nwaves)
81 84  
82 85  
83 86 FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
... ... @@ -119,7 +122,6 @@ ENDIF
119 122  
120 123 ;stop
121 124  
122   -
123 125 ;For spectrum data points, interpolate in log-log.
124 126 ;Linear interpolation leads to wrong values, in particular where few
125 127 ;wavelengths points exist in the model (long wavelengths).
... ...
src/idl/dustem_compute_sed.pro
  Diff comments   View file @3ff9cac
... ... @@ -62,7 +62,6 @@ spec = st.sed.em_tot * fact
62 62 ;if n_tags(!dustem_data.sed) gt 1 then begin ; executing only when !dustem_data.sed is set meaning that this allows for the creation of fake seds
63 63  
64 64  
65   -
66 65 scopes=tag_names((*!dustem_scope))
67 66 IF scopes[0] NE 'NONE' THEN BEGIN
68 67 ;IF ptr_valid(!dustem_plugin) THEN BEGIN
... ... @@ -70,6 +69,9 @@ IF scopes[0] NE 'NONE' THEN BEGIN
70 69 if total(strsplit((*(*!dustem_scope).(i)),'+',/extract) eq 'ADD_SED') then spec+=(*(*!dustem_plugin).(i))[*,0]
71 70 endfor
72 71 ENDIF
  72 +
  73 +
  74 +
73 75 ;endif
74 76 ;------------------------------------------
75 77  
... ... @@ -113,7 +115,7 @@ the_end:
113 115 ;VG : error in cc IRAS1 when no PAH
114 116 j=where(dustem_sed ne dustem_sed,count)
115 117 if count gt 0 then dustem_sed[j]=0.
116   -
  118 +;stop
117 119 RETURN,dustem_sed
118 120  
119 121 END
... ...
src/idl/dustem_compute_stokes.pro
  Diff comments   View file @3ff9cac
... ... @@ -51,7 +51,7 @@ out=0.
51 51  
52 52 If keyword_set(result) THEN result=st
53 53  
54   -fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(st.sed).wav)*1.e20/1.e7
  54 +fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/(st.sed).wav)*1.e20/1.e7
55 55 stI = st.sed.em_tot * fact ;This is Total intensity I
56 56 P = st.polsed.em_tot * fact ;This is Polarized intensity P
57 57  
... ... @@ -63,7 +63,6 @@ tes=where(finite(frac) eq 0)
63 63 frac(tes)=0.
64 64  
65 65  
66   -
67 66 FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
68 67  
69 68 IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN BEGIN
... ... @@ -78,7 +77,11 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
78 77 ENDFOR
79 78  
80 79 IF ~isa(Q_spec) && ~isa(U_spec) THEN BEGIN
81   - polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0,Nwaves)
  80 + ;maybe add a condition that takes into account the angles of the stokes parameters.
  81 + ;I believe it is best to hard code the angle information in here
  82 + ;this also includes its extraction?
  83 +
  84 + polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(!dustem_psi,Nwaves)
82 85 ;stop
83 86 ENDIF
84 87  
... ...
src/idl/dustem_init.pro
  Diff comments   View file @3ff9cac
... ... @@ -107,16 +107,16 @@ if keyword_set(pol) then begin
107 107 sed: ptr_new(), $
108 108 ext: ptr_new(), $
109 109 polext: ptr_new(), $
110   - polsed: ptr_new(), $
111   - polfrac: ptr_new(), $
  110 + ;polsed: ptr_new(), $ ;newly deprecated
  111 + ;polfrac: ptr_new(), $ ;newly deprecated
112 112 qsed: ptr_new(), $
113 113 used: ptr_new(), $
114 114 qext: ptr_new(), $
115 115 uext: ptr_new() $
116 116 }
117 117  
118   - rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.,qsed:0.,used:0.,qext:0.,uext:0.}
119   -
  118 + ;rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.,qsed:0.,used:0.,qext:0.,uext:0.} ;polsed and polfrac newly deprecated
  119 + rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,qsed:0.,used:0.,qext:0.,uext:0.}
120 120 endif else begin
121 121 defsysv, '!dustem_data', { $ ;Data to fit
122 122 sed: ptr_new(), $
... ... @@ -134,6 +134,24 @@ ENDFOR
134 134 instr+='}'
135 135 toto=execute(instr)
136 136  
  137 +
  138 +defsysv,'!dustem_show', { $ ;Data to fit
  139 + sed: ptr_new(), $
  140 + ext: ptr_new(), $
  141 + polext: ptr_new(), $
  142 + polsed: ptr_new(), $
  143 + polfrac: ptr_new(), $
  144 + psi_em: ptr_new(), $
  145 + psi_ext: ptr_new(), $
  146 + qsed: ptr_new(), $
  147 + used: ptr_new(), $
  148 + qext: ptr_new(), $
  149 + uext: ptr_new() $
  150 + }
  151 +
  152 +;defsysv, '!dustem_psi',0. ;TO BE DEPRECTED... CHECK IF U STARTED IMPLEMENTING THIS ;0 is default. Meanin No PSI value(s) detected in the data xcat file.
  153 +
  154 +
137 155 ;IDL> help,*!dustem_data,/str
138 156 ;** Structure <2173b20>, 4 tags, length=448, data length=448, refs=1:
139 157 ; INSTRU_NAMES STRING Array[14]
... ... @@ -284,6 +302,8 @@ ENDELSE
284 302 st_model=dustem_read_all(dir_in,/silent,kwords=kwords)
285 303 !dustem_params=ptr_new(st_model)
286 304  
  305 +;ADD A last wavelength here.
  306 +
287 307 ;=== create dynamical storage if needed
288 308 IF not file_test(!dustem_dat,/dir) THEN BEGIN
289 309 spawn,'mkdir '+!dustem_dat
... ...
src/idl/dustem_mpfit_run.pro
  Diff comments   View file @3ff9cac
... ... @@ -58,8 +58,8 @@ st_model=dustem_read_all(!dustem_soft_dir)
58 58 chi2_sed = 0
59 59 chi2_ext = 0
60 60 chi2_polext = 0
61   -chi2_polsed = 0
62   -chi2_polfrac = 0
  61 +;chi2_polsed = 0
  62 +;chi2_polfrac = 0
63 63 chi2_qsed =0
64 64 chi2_used =0
65 65 chi2_qext =0
... ... @@ -67,8 +67,8 @@ chi2_uext =0
67 67 rchi2_sed = 0
68 68 rchi2_ext = 0
69 69 rchi2_polext = 0
70   -rchi2_polsed = 0
71   -rchi2_polfrac = 0
  70 +;rchi2_polsed = 0
  71 +;rchi2_polfrac = 0
72 72 rchi2_qsed =0
73 73 rchi2_used =0
74 74 rchi2_qext =0
... ... @@ -76,8 +76,8 @@ rchi2_uext =0
76 76 wrchi2_sed = 0
77 77 wrchi2_ext = 0
78 78 wrchi2_polext = 0
79   -wrchi2_polsed = 0
80   -wrchi2_polfrac = 0
  79 +;wrchi2_polsed = 0
  80 +;wrchi2_polfrac = 0
81 81 wrchi2_qsed =0
82 82 wrich2_used =0
83 83 wrchi2_qext =0
... ... @@ -85,8 +85,8 @@ wrich2_uext =0
85 85 n_sed = 0
86 86 n_ext = 0
87 87 n_polext = 0
88   -n_polsed = 0
89   -n_polfrac = 0
  88 +;n_polsed = 0
  89 +;n_polfrac = 0
90 90 n_qsed =0
91 91 n_used =0
92 92 n_qext =0
... ... @@ -110,7 +110,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
110 110 CASE tagnames(ind_data_tag(i_tag)) OF
111 111  
112 112 'SED' : BEGIN
113   - dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  113 + dustem_sed = dustem_compute_sed(p_dim,st)
114 114 chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
115 115 n_sed = n_elements((*!dustem_data.sed).values)
116 116 rchi2_sed = chi2_sed / n_sed
... ... @@ -130,90 +130,13 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
130 130 (*!dustem_fit).rchi2=rchi2_sed
131 131 ; Plot if needed
132 132 IF !dustem_show_plot NE 0 THEN BEGIN
133   -
134   - ;title='DUSTEM WRAP('+strtrim(string(win),2)+')'
  133 +
135 134 window,win,title='DUSTEM WRAP (SED)'
136 135  
137 136 ;cgwindow,WMULTI=[2,2,1],/CURRENT
138   - ;I.C
139   - ;undeprecated the use of polfrac here
140   - ;**********
141   - ;I added this block so I can plot polfrac inside the sed window
142   - ;I will probably remove it
143   - ;**********
144   - ;test if polfrac is in tagnames:
145   - ;another problem is that sed and pol data don't have the same filters.
146   - varvar=where(tagnames eq 'POLFRAC', testfpol)
147   - if testfpol ne 0 then begin
148   - if !run_lin then begin
149   - dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here.
150   - ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat
151   - ind_fspec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_fspec) ;locating the spectrum points in polfrac data. Later testing will have to be on the wavelengths
152   -
153   - ;dustem_sed and dustem_polsed might not have the same length thus this block.
154   - ;I'm not so sure about this...
155   -
156   - If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin
157   - if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin
158   - dustem_polsed_x = dustem_polsed
159   - dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified
160   - ;matching the filter data points
161   - IF count_filt ne 0 then begin
162   - for i=0L,count_filt-1 do begin
163   - j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
164   - if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0))
165   - endfor
166   - endif
167   - ;matching the spectrum data points
168   - IF count_fspec ne 0 then begin
169   - for i=0L,count_fspec-1 do begin
170   - j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
171   - if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0))
172   - endfor
173   - endif
174   - endif ELSE begin
175   -
176   - ;I think I will have errors here!! The thing is I'm trying to make this work without (*!dustem_data.polsed) which is hard.
177   -
178   - dustem_polsed_x = dustem_sed
179   - dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified
180   - ;matching the filter data points
181   - IF count_filt ne 0 then begin
182   - for i=0L,count_filt-1 do begin
183   - j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
184   - if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0))
185   - endfor
186   - endif
187   - ;matching the spectrum data points
188   - IF count_fspec ne 0 then begin
189   - for i=0L,count_fspec-1 do begin
190   - j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
191   - if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0))
192   - endfor
193   - endif
194   -
195   - ENDELSE
196   - endif ELSE BEGIN
197   - dustem_polsed_x = dustem_polsed
198   - dustem_sed_x = dustem_sed
199   - ENDELSE
200   -
201   - dustem_polfrac = dustem_polsed_x/dustem_sed_x
202   - ;stop
203   - ;stop
204   - chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
205   - n_polfrac = n_elements((*!dustem_data.polfrac).values)
206   - rchi2_polfrac = chi2_polfrac / n_polfrac
207   - wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac
208   - message,"chi2 POLFRAC = "+strtrim(chi2_polfrac,2)+" rchi2 POLFRAC = "+strtrim(rchi2_polfrac,2),/continue ;," weighted rchi2 POLFRAC=",wrchi2_polfrac
209   - dustem_out = [ dustem_out, dustem_polfrac ]
210   - fpol=dustem_polfrac ; problem because extra structure does not get updated so you actually need the 'extra' procedure you haven't finihed coding.
211   - endif
212   - endif
213   - ;stop
  137 +
214 138 dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2,fpol=fpol
215 139  
216   - ;stop
217 140 ENDIF
218 141  
219 142 dustem_out = [ dustem_out, dustem_sed ]
... ... @@ -355,11 +278,67 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
355 278  
356 279 END
357 280  
358   -;USE OF POLFRAC DEPRECATED---------------------------------------
359   -
360   -; 'POLFRAC': BEGIN
361   -;
362   -; ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_POLSED AND DUSTEM_COMPUTE_SED
  281 +; 'POLFRAC': BEGIN ;Deprecate this. Fitting POLFRAC does not make any sense. It is a biased quantity.
  282 +; if !run_lin then begin
  283 +; dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here.
  284 +; ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat
  285 +; ind_fspec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_fspec) ;locating the spectrum points in polfrac data. Later testing will have to be on the wavelengths
  286 +;
  287 +; ;dustem_sed and dustem_polsed might not have the same length thus this block.
  288 +;
  289 +; If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin
  290 +; if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin
  291 +; dustem_polsed_x = dustem_polsed
  292 +; dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified
  293 +; ;matching the filter data points
  294 +; IF count_filt ne 0 then begin
  295 +; for i=0L,count_filt-1 do begin
  296 +; j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
  297 +; if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0))
  298 +; endfor
  299 +; endif
  300 +; ;matching the spectrum data points
  301 +; IF count_fspec ne 0 then begin
  302 +; for i=0L,count_fspec-1 do begin
  303 +; j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
  304 +; if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0))
  305 +; endfor
  306 +; endif
  307 +; endif ELSE begin
  308 +;
  309 +; dustem_polsed_x = dustem_sed
  310 +; dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified
  311 +; ;matching the filter data points
  312 +; IF count_filt ne 0 then begin
  313 +; for i=0L,count_filt-1 do begin
  314 +; j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
  315 +; if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0))
  316 +; endfor
  317 +; endif
  318 +; ;matching the spectrum data points
  319 +; IF count_fspec ne 0 then begin
  320 +; for i=0L,count_fspec-1 do begin
  321 +; j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
  322 +; if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0))
  323 +; endfor
  324 +; endif
  325 +;
  326 +; ENDELSE
  327 +; endif ELSE BEGIN
  328 +; dustem_polsed_x = dustem_polsed
  329 +; dustem_sed_x = dustem_sed
  330 +; ENDELSE
  331 +;
  332 +; dustem_polfrac = dustem_polsed_x/dustem_sed_x
  333 +; chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
  334 +; n_polfrac = n_elements((*!dustem_data.polfrac).values)
  335 +; rchi2_polfrac = chi2_polfrac / n_polfrac
  336 +; wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac
  337 +; message,"chi2 POLFRAC = "+strtrim(chi2_polfrac,2)+" rchi2 POLFRAC = "+strtrim(rchi2_polfrac,2),/continue ;," weighted rchi2 POLFRAC=",wrchi2_polfrac
  338 +; dustem_out = [ dustem_out, dustem_polfrac ]
  339 +; ;fpol=dustem_polfrac ; problem because extra structure does not get updated so you actually need the 'extra' procedure you haven't finihed coding.
  340 +; endif
  341 +;OLD DEPRECATED VERSION.
363 342 ;
364 343 ; IF !run_lin then begin
365 344 ;
... ... @@ -401,45 +380,48 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
401 380  
402 381  
403 382 ;-----------------------------------------------------------------
404   - 'POLSED': BEGIN
405   -
406   - ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
407   -
408   - IF !run_lin THEN BEGIN
409   -
410   - dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
411   -
412   - chi2_polsed = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
413   - n_polsed = n_elements((*!dustem_data.polsed).values)
414   - rchi2_polsed = chi2_polsed / n_polsed
415   - wrchi2_polsed = rchi2_polsed * !fit_rchi2_weight.polsed
416   - message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+" rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
417   - print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
418   - ;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
419   -
420   - dustem_out = [ dustem_out, dustem_polsed ] ; classic command
421   - ;-----------------------------------------------
422   -
423   - ; Plot if needed
424   - IF !dustem_show_plot NE 0 THEN BEGIN
425   - win+=1
  383 +; 'POLSED': BEGIN
  384 +;
  385 +; ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
  386 +;
  387 +; IF !run_lin THEN BEGIN
  388 +; ;make polsed only used for when polsed is fitted. Otherwise it's the fitting of Q_SED AND USED that fits the polarization data
  389 +; dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  390 +;
  391 +; chi2_polsed = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
  392 +; n_polsed = n_elements((*!dustem_data.polsed).values)
  393 +; rchi2_polsed = chi2_polsed / n_polsed
  394 +; wrchi2_polsed = rchi2_polsed * !fit_rchi2_weight.polsed
  395 +; message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+" rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  396 +; print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
  397 +; ;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
  398 +;
  399 +; dustem_out = [ dustem_out, dustem_polsed ] ; classic command
  400 +; ;-----------------------------------------------
  401 +;
  402 +; ; Plot if needed
  403 +; IF !dustem_show_plot NE 0 THEN BEGIN
  404 +; win+=1
  405 +; ;
  406 +; xr=[10,8e3]
  407 +; yr=[1e-34,1e-28]
426 408 ;
427   - xr=[10,8e3]
428   - yr=[1e-34,1e-28]
429   -
430   - dustem_plot_polsed,st,p_dim,dustem_polsed,aligned=st_model.align.grains.aligned,win=win
431   -; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1
432   -; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2
433   -
434   -
435   -
436   - ENDIF
437   - ENDIF
438   - END
  409 +; dustem_plot_polsed,st,p_dim,dustem_polsed,aligned=st_model.align.grains.aligned,win=win
  410 +; ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1
  411 +; ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2
  412 +;
  413 +;
  414 +;
  415 +; ENDIF
  416 +; ENDIF
  417 +; END
439 418  
440   - 'QSED': BEGIN
  419 + 'QSED': BEGIN
441 420  
442 421 toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
  422 + ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
  423 + dustem_qsed(-1)=dustem_qsed(-2)
  424 + stop
443 425 chi2_qsed =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2)
444 426 n_qsed = n_elements((*!dustem_data.qsed).values)
445 427 rchi2_qsed = chi2_qsed / n_qsed
... ... @@ -455,16 +437,15 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
455 437 fact1 = 1.00E-20;(3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40 ; for data plotting
456 438  
457 439 ; fact=1/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.00e17
458   -; fact1=1.E-20
  440 +; fact1=1.E-20
459 441 ;
460 442 win+=1
461 443 window ,win ,title='DUSTEM WRAP (QSED)'
462 444 titstq='Stokes Q Polarization Intensity'
463 445 ytitstq=textoidl('Q_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')
464 446 xtit=textoidl('\lambda (\mum)')
465   - xr=[10,2e4]
466   - normpol=Q_spec*fact
467   - normpol1=dustem_qsed*fact1
  447 + xr=[10,1e6]
  448 +
468 449 ;normpol1=(*!dustem_data.qsed).values*fact1
469 450 yr=[1e-28,1e-17]
470 451 ;yr=[1e-34,1e-23]
... ... @@ -473,10 +454,14 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
473 454  
474 455 ylog=1
475 456 if countq ne 0 then begin
476   - ylog=0
477   - yr=[-2e-31,5e-31]
  457 +; ylog=0
  458 +; yr=[-1E-23,-1e-19];e-31
  459 + Q_spec=-Q_spec
  460 + dustem_qsed =-dustem_qsed
  461 + ((*!dustem_data.qsed).values)=-((*!dustem_data.qsed).values)
478 462 endif
479   -
  463 + normpol=Q_spec*fact
  464 + normpol1=dustem_qsed*fact1
480 465  
481 466  
482 467 cgplot,st.polsed.wav,Q_spec*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)'
... ... @@ -491,46 +476,46 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
491 476 err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1/normpol1,color=cgColor('Tomato')
492 477  
493 478 endif
494   -
495   - END
  479 + END
496 480  
497 481 'USED': BEGIN
498 482  
499 483 toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
500   - chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2)
  484 + stop
  485 + chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2)
501 486 n_used = n_elements((*!dustem_data.used).values)
502 487 rchi2_used = chi2_used / n_used
503   - wrchi2_used = rchi2_used * !fit_rchi2_weight.used
  488 + wrchi2_used = rchi2_used * !fit_rchi2_weight.used
504 489 dustem_out = [ dustem_out, dustem_used] ; classic command
505 490 message,"chi2 USED = "+strtrim(chi2_used,2)+" rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
506   - print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2
507   -
508   -
509   - if !dustem_show_plot ne 0 then begin
510   - fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40
511   - fact1 =1.00E-20;(3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting
  491 + print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2
  492 +
  493 +
  494 + if !dustem_show_plot ne 0 then begin
  495 + fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40
  496 + fact1 =1.00E-20;(3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting
512 497  
513   - win+=1
  498 + win+=1
514 499 window ,win ,title='DUSTEM WRAP (USED)'
515 500 titstu='Stokes U Polarization Intensity'
516 501 ytitstu=textoidl('U_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2')
517 502 xtit=textoidl('\lambda (\mum)')
518   - xr=[10,6e4]
519   - normpol2=U_spec*fact
  503 + xr=[10,1e6]
  504 + normpol2=U_spec*fact
520 505 normpol3=dustem_used*fact1
521 506 ;yr=[min(normpol2)-min(normpol2)/10,max(normpol2)+max(normpol2)*10]
522   -
523   - ;conditions here are on computed spectra and they should be on the data itself.
524   - ;This even impairs the visualization of the data. Correct this!!!!!
525 507  
526   - yr=[1e-28,1e-17];[1e-34,1e-23]
527   - indu=where(U_spec lt 0, countu)
528   - ;indu=where((*!dustem_data.used).values lt 0, countu)
529   - ylog=1
530   - if countu ne 0 then begin
531   - ylog=0
532   - yr=[-2e-31,5e-31]
533   - endif
  508 + ;conditions here are on computed spectra and they should be on the data itself.
  509 + ;This even impairs the visualization of the data. Correct this!!!!!
  510 +
  511 + yr=[1e-28,1e-17];[1e-34,1e-23]
  512 + indu=where(U_spec lt 0, countu)
  513 + ;indu=where((*!dustem_data.used).values lt 0, countu)
  514 + ylog=1
  515 + if countu ne 0 then begin
  516 + ylog=0
  517 + yr=[-2e-31,5e-31]
  518 + endif
534 519  
535 520 cgplot,st.polsed.wav,U_spec*fact,xtit='',ytit=ytitstu,tit=titstu,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
536 521 cgoplot,st.polsed.wav,U_spec*fact,color='black'
... ... @@ -547,10 +532,10 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
547 532  
548 533 ;cgwindow, 'plot',
549 534  
550   - endif
551   - END
552   -
553   -
  535 + endif
  536 + END
  537 +;
  538 +;
554 539  
555 540  
556 541 'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION
... ... @@ -670,16 +655,16 @@ if !fit_rchi2_weight.polext ne 0 then begin
670 655 total_chi2 += chi2_polext
671 656 total_rchi2 += rchi2_polext
672 657 endif
673   -if !fit_rchi2_weight.polsed ne 0 then begin
674   - total_n += n_polsed
675   - total_chi2 += chi2_polsed
676   - total_rchi2 += rchi2_polsed
677   -endif
678   -if !fit_rchi2_weight.polfrac ne 0 then begin
679   - total_n += n_polfrac
680   - total_chi2 += chi2_polfrac
681   - total_rchi2 += rchi2_polfrac
682   -endif
  658 +; if !fit_rchi2_weight.polsed ne 0 then begin
  659 +; total_n += n_polsed
  660 +; total_chi2 += chi2_polsed
  661 +; total_rchi2 += rchi2_polsed
  662 +; endif
  663 +; if !fit_rchi2_weight.polfrac ne 0 then begin
  664 +; total_n += n_polfrac
  665 +; total_chi2 += chi2_polfrac
  666 +; total_rchi2 += rchi2_polfrac
  667 +; endif
683 668 if !fit_rchi2_weight.qsed ne 0 then begin
684 669 total_n += n_qsed
685 670 total_chi2 += chi2_qsed
... ... @@ -708,8 +693,8 @@ total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $
708 693 + chi2_sed * !fit_rchi2_weight.sed
709 694 if !run_pol then begin
710 695 total_wchi2 += chi2_polext * !fit_rchi2_weight.polext $
711   - + chi2_polsed * !fit_rchi2_weight.polsed $
712   - + chi2_polfrac * !fit_rchi2_weight.polfrac $
  696 +; + chi2_polsed * !fit_rchi2_weight.polsed $
  697 +; + chi2_polfrac * !fit_rchi2_weight.polfrac $
713 698 + chi2_qsed * !fit_rchi2_weight.qsed $
714 699 + chi2_used * !fit_rchi2_weight.used $
715 700 + chi2_qext * !fit_rchi2_weight.qext $
... ... @@ -720,9 +705,9 @@ message,&#39;======================================================================&#39;
720 705 print,"Total weighted chi2 driving mpfitfun = ", total_wchi2
721 706 print,"Nb of data points = ", total_n
722 707 print,"DOF = ", DOF
723   -print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed,rchi2_polext,rchi2_polsed,rchi2_polfrac
724   -print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed,wrchi2_polext,wrchi2_polsed,wrchi2_polfrac
725   -if !run_pol then print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed,!fit_rchi2_weight.polext,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $
  708 +print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed,rchi2_polext;,rchi2_polsed,rchi2_polfrac
  709 +print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed,wrchi2_polext;,wrchi2_polsed,wrchi2_polfrac
  710 +if !run_pol then print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed,!fit_rchi2_weight.polext $;,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $
726 711 else print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed
727 712 print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF)
728 713 print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF)
... ...
src/idl/dustem_plot_fit_sed.pro
  Diff comments   View file @3ff9cac
... ... @@ -46,6 +46,9 @@
46 46 ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
47 47 ;-
48 48  
  49 +
  50 +;NOta bene: only delete the cgoplot commands
  51 +
49 52 ;stop
50 53  
51 54 IF keyword_set(help) THEN BEGIN
... ... @@ -136,7 +139,7 @@ IF keyword_set(ytit) THEN ytit = ytit ELSE ytit = textoidl(&#39;Brightness/N_H (MJy/
136 139  
137 140  
138 141 ;deffo define a title variable here. The procedures are not communicating with each other.
139   -
  142 +;pos=cgLayout()
140 143 ;###############################
141 144 if keyword_set(fpol) then begin
142 145 if !run_lin then begin
... ...
src/idl/dustem_plot_polsed.pro
  Diff comments   View file @3ff9cac
1   - PRO dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win
  1 +PRO dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win
2 2  
3 3  
4 4 IF keyword_set(ps) THEN BEGIN
... ...
src/idl/dustem_plugin_modify_dust_pol.pro
  Diff comments   View file @3ff9cac
... ... @@ -65,7 +65,6 @@ P=((st.polsed).em_tot)*fact
65 65 I=((st.sed).em_tot)*fact
66 66  
67 67  
68   -
69 68 Nwaves=(size(I))[1]
70 69  
71 70 frac=P/I*smallp
... ...
src/idl/dustem_plugin_stellar_population.pro
  Diff comments   View file @3ff9cac
... ... @@ -58,7 +58,7 @@ countf = 0.
58 58 countg = 0.
59 59 ;==============================
60 60  
61   -
  61 +;REPLACE THIS WHITH A PLUGIN PARAMETER THAT CAN BE FITTED
62 62 ;===========Specifying whether to add the Mathis field (and its scaling via G0) to the composite ISRF===================
63 63 ismathis=0 ;this value will have to be manually changed at the present moment
64 64 defsysv, '!ismathis', ismathis
... ... @@ -3177,14 +3177,12 @@ c(Ncomments-2)=&#39;# Nbr of points&#39;
3177 3177 c(Ncomments-1)='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'
3178 3178  
3179 3179 FOR i=0L,n_elements(comp_pop.popid)-1 DO BEGIN ; Looping over all the stellar populations
3180   -
3181   - ;The initial procedure had omega multiplied by a !pi factor. Its presence in the Planck (Astron) procedure makes for a good reason to discuss this with J.P.
3182   -
  3180 +
3183 3181 omega = ((comp_pop.radius)[i]/((comp_pop.distance)[i]*pc2cm))^2 ; Dilution factor of a stellar population
3184 3182  
3185   - Inu = planck(wave_angstrom,(comp_pop.temperature)[i])/(4.*!pi)*(wave_angstrom)^2/c2a ; ergs/cm2/s/Hz/sr
  3183 + Inu = planck(wave_angstrom,(comp_pop.temperature)[i])*(wave_angstrom)^2/c2a ; ergs/cm2/s/Hz/sr
3186 3184  
3187   - stellar_component=stellar_component+(comp_pop.nstars)[i]*omega*Inu
  3185 + stellar_component=stellar_component+(comp_pop.nstars)[i]*omega*Inu ;no *4 because of Lambert's cosine law.
3188 3186  
3189 3187 ;Rest of the lines of the new composite ISRF.DAT file
3190 3188  
... ... @@ -3199,10 +3197,15 @@ FOR i=0L,n_elements(comp_pop.popid)-1 DO BEGIN ; Looping over all the stellar po
3199 3197 ENDELSE
3200 3198 ENDFOR
3201 3199  
  3200 +;make sure which g0 is used!!!
  3201 +;alwayd devide by g0 even though it equals one.
  3202 +;
3202 3203  
3203   -IF !ismathis THEN st.isrf=ma_isrf.isrf+stellar_component/(st.gas.G0) ELSE st.isrf=stellar_component;/((*!dustem_params).gas.G0)
  3204 +IF !ismathis THEN st.isrf=ma_isrf.isrf+stellar_component/(st.gas.G0) ELSE st.isrf=stellar_component;/(st.gas.G0)
3204 3205  
3205   -;don't mind my tests I'm just making sure the ISRF file changes at all
  3206 +;stop
  3207 +
  3208 +;don't mind my tests I'm just making sure if the ISRF file changes at all
3206 3209 print,'stellar_component is:'
3207 3210  
3208 3211 print, stellar_component
... ... @@ -3229,6 +3232,10 @@ close,unit
3229 3232 free_lun,unit
3230 3233  
3231 3234 out=st.isrf
  3235 +
  3236 +;THIS BLOCK IS FALSE AND NEEDS TO BE CORRECTED. ASK JP WHAT HE WANTS FOR A DEFAULT RUN
  3237 +;(a run with no parameters for instance)
  3238 +
3232 3239 ;============This block creates a composite (9V) stellar population structure for when the function isn't used as a plugin===================
3233 3240 ENDIF ELSE BEGIN
3234 3241  
... ...
src/idl/dustem_sed_plot.pro
  Diff comments   View file @3ff9cac
1 1 PRO dustem_sed_plot,p_min,_extra=_extra,function_name=function_name,pol=pol,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png,fpol=fpol
  2 +;THIS procedure is incomplete. It needs to be updated.
2 3  
3 4 IF not keyword_set(function_name) THEN BEGIN
4 5 ;Run dustem with as an interface to mpfitfun; what does this mean? This whole procedure will have to be changed
... ...
src/idl/dustem_set_data.pro
  Diff comments   View file @3ff9cac
1   -FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,fpol=fpol,rchi2_weight=rchi2_weight,f_HI=f_HI
  1 +FUNCTION dustem_set_data,st_sed_fit=st_sed_fit,st_sed_show=st_sed_show,st_ext_fit=st_ext_fit,st_ext_show=st_ext_show,rchi2_weight=rchi2_weight,f_HI=f_HI,help=help
  2 +
  3 +
  4 +;update the help
  5 +
2 6  
3 7 ;+
4 8 ; NAME:
... ... @@ -8,7 +12,7 @@ FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,fpol=fpol,rchi2_weight=rchi2_
8 12 ; CATEGORY:
9 13 ; Dustem
10 14 ; CALLING SEQUENCE:
11   -; dustem_set_data,spec[,help=]
  15 +; dustem_set_data(data_sed=data_sed,data_ext=data_ext,[/help])
12 16 ; INPUTS:
13 17 ; spec: array of structures containing a SED
14 18 ; OPTIONAL INPUT PARAMETERS:
... ... @@ -42,8 +46,8 @@ FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,fpol=fpol,rchi2_weight=rchi2_
42 46 ;
43 47 ; V. Guillet (2012) : dustem_set_data is turned into a function which can be used
44 48 ; indifferently for sed, ext, polext
45   -; I. Choubani (2022): Adjusting the function to meet the format criteria set by the current release.
46   -; The function uses two main input structures: SED and EXT. Polarization data is included within each.
  49 +; I. Choubani (2022): Introducing the dustem_check_data procedure where the new format outlined by the latest release is taken into account.
  50 +; This function uses two main input structures: DATA_SED and DATA_EXT. Polarization data is included within each.
47 51 ;-
48 52  
49 53  
... ... @@ -52,132 +56,240 @@ IF keyword_set(help) THEN BEGIN
52 56 goto,the_end
53 57 ENDIF
54 58  
55   -;Cleaninig data structure
56   -!dustem_data.sed = ptr_new()
57   -!dustem_data.ext = ptr_new()
  59 +;Cleaninig data/showing structure
  60 +for i=0L,n_tags(!dustem_data)-1 do begin
  61 + !dustem_data.(i) = ptr_new()
  62 + !dustem_show.(i) = ptr_new()
  63 +endfor
  64 +
  65 +xx=dustem_check_data(data_sed=st_sed_fit,data_ext=st_ext_fit,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,psi_em=psi_em,qsed=qsed,used=used,qext=qext,uext=uext,psi_ext=psi_ext)
58 66  
  67 +;outputs whatever the input structure has, into organized seperate structures formatted for the fitting process.
59 68  
60   -;SETTING THE DATA STRUCTURES FOR THE FITTING PROCESS - SED
  69 +;If SED/EXT data is present in the structure it is fitted automatically.
  70 +;If the user wished to not fit it they will have to null the corresponding fields in the input fitting structure.
  71 +IF isa(sed) THEN !dustem_data.sed = ptr_new(sed)
  72 +IF isa(ext) THEN !dustem_data.ext = ptr_new(ext)
61 73  
62   -IF keyword_set(sed) THEN BEGIN
  74 +
  75 +;####################
  76 +;OLD VARIABLE. KEPT IT.
  77 +; If f_HI is specified, multiply the data by f_HI
  78 +If keyword_set(f_HI) and ptr_valid(!dustem_data.sed) then begin
  79 + if f_HI gt 0 then (*!dustem_data.sed).values *= f_HI else f_HI = 1.
  80 +endif else f_HI = 1.
  81 +if keyword_set(f_HI) and ptr_valid(!dustem_data.ext) then begin
  82 + if f_HI gt 0 then (*!dustem_data.ext).values *= f_HI else f_HI = 1.
  83 +endif else f_HI = 1.
  84 +defsysv,'!dustem_f_HI',f_HI
  85 +;####################
  86 +
  87 +;NOTA BENE: The user is bound to fit Q and U together.
  88 +if !run_pol then begin
63 89  
  90 + if keyword_set(st_sed_fit) then begin
64 91  
65   - ;maybe to activate fpol here you need to set it to one? CHECK THIS OUT.
66   - if keyword_set(fpol) then xx=dustem_check_data(sed=sed,ist=ist,pst=pst,fpol=fpol,fpst=fpst,qst=qst,ust=ust) else xx=dustem_check_data(sed=sed,ist=ist,pst=pst,fpst=fpst,qst=qst,ust=ust)
67   - ;stop
68   - ;=== fill !dustem_data
69   - IF isa(ist) THEN !dustem_data.sed = ptr_new(ist) ;ELSE !dustem_data.sed = ptr_new()
70   - ; If f_HI is specified, multiply the data by f_HI
71   - If keyword_set(f_HI) and ptr_valid(!dustem_data.sed) then begin
72   - if f_HI gt 0 then (*!dustem_data.sed).values *= f_HI else f_HI = 1.
73   - endif else f_HI = 1.
74   - defsysv,'!dustem_f_HI',f_HI
75   - if !run_pol then begin
76   - If isa(pst) then BEGIN
77   - ;=== fill !dustem_data
78   - !dustem_data.polsed = ptr_new(pst)
79   - ; If f_HI is specified, multiply the data by f_HI
80   - if keyword_set(f_HI) then if f_HI gt 0 then $
81   - (*!dustem_data.polsed).values *= f_HI
82   - ENDIF
83   - If isa(fpst) then BEGIN
84   - ;=== fill !dustem_data
85   - !dustem_data.polfrac = ptr_new(fpst)
86   - ; If f_HI is specified, multiply the data by f_HI
87   - if keyword_set(f_HI) then if f_HI gt 0 then $
88   - (*!dustem_data.polfrac).values *= f_HI
89   - ENDIF
90   - IF isa(qst) THEN BEGIN
91   - ;=== fill !dustem_data
92   - !dustem_data.qsed = ptr_new(qst)
93   - ; If f_HI is specified, multiply the data by f_HI
  92 + teststks = isa(qsed) and isa(used)
  93 +
  94 + IF (not teststks) then message, 'Stokes parameters are not set correctly. Please check the input structure. Aborting...'
  95 +
  96 + ;EMISSION
  97 + if teststks then nowfitting = 'QSED and USED'
  98 +
  99 + message, 'Polarization component(s) to fit: '+nowfitting ,/continue
  100 +
  101 + IF isa(qsed) THEN BEGIN
  102 + !dustem_data.qsed = ptr_new(qsed)
94 103 if keyword_set(f_HI) then if f_HI gt 0 then $
95 104 (*!dustem_data.qsed).values *= f_HI
96 105 ENDIF
97   - if isa(ust) then BEGIN
98   - ;=== fill !dustem_data
99   - !dustem_data.used = ptr_new(ust)
100   - ; If f_HI is specified, multiply the data by f_HI
  106 + if isa(used) then BEGIN
  107 + !dustem_data.used = ptr_new(used)
101 108 if keyword_set(f_HI) then if f_HI gt 0 then $
102 109 (*!dustem_data.used).values *= f_HI
103 110 ENDIF
104 111 endif
105   -ENDIF
106   -
107   -;SETTING THE DATA STRUCTURES FOR THE FITTING PROCESS - EXT
108   -
109   -IF keyword_set(ext) THEN BEGIN
110   - xx=dustem_check_data(ext=ext,exst=exst,pexst=pexst,qexst=qexst,uexst=uexst)
111   - ;=== fill !dustem_data
112   - IF isa(exst) THEN !dustem_data.ext = ptr_new(exst) ;ELSE !dustem_data.ext = ptr_new()
113   -
114   - ; If f_HI is specified, multiply the data by f_HI
115   - if keyword_set(f_HI) and ptr_valid(!dustem_data.ext) then begin
116   - if f_HI gt 0 then (*!dustem_data.ext).values *= f_HI else f_HI = 1.
117   - endif else f_HI = 1.
118   - defsysv,'!dustem_f_HI',f_HI
119 112  
120   - if !run_pol then begin
121   - if isa(pexst) then BEGIN
122   - !dustem_data.polext = ptr_new(pexst)
123   - ; If f_HI is specified, multiply the data by f_HI
  113 + if keyword_set(st_ext_fit) then begin
  114 + ;Two tests that determine what the user will eventually fit
  115 +
  116 + testexstks = isa(qext) and isa(uext)
  117 + testpexstks = (isa(qext) and isa(uext) and isa(polext))
  118 +
  119 + IF (not testexstks and not isa(polext)) then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
  120 +
  121 + IF testpexstks THEN $
  122 + message, 'Presence of both extinction polarization and stokes parameter values. Fitting structure is not correctly filtered. Aborting...'
  123 +
  124 + ;EXTINCTION - allowed for the default fitting of both because I do not know if the same formalism should be applied to the extinction data because it also implies deviding extinction data into two components.
  125 + if testexstks then nowfitting = 'QEXT and UEXT' else begin ;Serkowski is default
  126 + if isa(polext) then nowfitting = 'POLEXT'
  127 + endelse
  128 +
  129 + message, 'Extinction polarization component(s) to fit: '+nowfitting ,/continue
  130 + if isa(polext) then BEGIN
  131 + !dustem_data.polext = ptr_new(polext)
124 132 if keyword_set(f_HI) then if f_HI gt 0 then $
125 133 (*!dustem_data.polext).values *= f_HI
126 134 ENDIF
127   - IF isa(qexst) THEN BEGIN
128   - ;=== fill !dustem_data
129   - !dustem_data.qext = ptr_new(qexst)
130   - ; If f_HI is specified, multiply the data by f_HI
  135 + IF isa(qext) THEN BEGIN
  136 + !dustem_data.qext = ptr_new(qext)
131 137 if keyword_set(f_HI) then if f_HI gt 0 then $
132 138 (*!dustem_data.qext).values *= f_HI
133 139 ENDIF
134   - if isa(uexst) then BEGIN
135   - ;=== fill !dustem_data
136   - !dustem_data.uext = ptr_new(uexst)
137   - ; If f_HI is specified, multiply the data by f_HI
  140 + if isa(uext) then BEGIN
  141 + !dustem_data.uext = ptr_new(uext)
138 142 if keyword_set(f_HI) then if f_HI gt 0 then $
139 143 (*!dustem_data.uext).values *= f_HI
140 144 ENDIF
141 145 endif
142   -ENDIF
  146 +endif
  147 +
  148 +
  149 +if not keyword_set(st_sed_show) then st_sed_show=st_sed_fit
143 150  
144   -!fit_rchi2_weight.sed=1.
145   -!fit_rchi2_weight.ext=1.
  151 +;Setting of the structure that will be displayed
  152 +yy=dustem_check_data(data_sed=st_sed_show,data_ext=st_ext_show,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,qsed=qsed,used=used,qext=qext,uext=uext)
146 153  
147   -If !run_pol THEN BEGIN
148   - ;SED
149   - IF isa(pst) THEN !fit_rchi2_weight.polsed=1.
150   - IF isa(fpst) THEN !fit_rchi2_weight.polfrac=1.
151   - IF isa(qst) THEN !fit_rchi2_weight.qsed=1.
152   - IF isa(ust) THEN !fit_rchi2_weight.used=1.
  154 +IF isa(sed) THEN !dustem_show.sed = ptr_new(sed)
  155 +IF isa(ext) THEN !dustem_show.ext = ptr_new(ext)
  156 +
  157 +
  158 +;####################
  159 +;OLD VARIABLE. KEPT IT.
  160 +; If f_HI is specified, multiply the data by f_HI
  161 +If keyword_set(f_HI) and ptr_valid(!dustem_show.sed) then begin
  162 + if f_HI gt 0 then (*!dustem_show.sed).values *= f_HI else f_HI = 1.
  163 +endif else f_HI = 1.
  164 +if keyword_set(f_HI) and ptr_valid(!dustem_show.ext) then begin
  165 + if f_HI gt 0 then (*!dustem_show.ext).values *= f_HI else f_HI = 1.
  166 +endif else f_HI = 1.
  167 +defsysv,'!dustem_f_HI',f_HI
  168 +;####################
  169 +
  170 +if !run_pol then begin
153 171  
154   - ;EXT
155   - IF isa(pexst) THEN !fit_rchi2_weight.polext=1.
156   - IF isa(qexst) THEN !fit_rchi2_weight.qext=1.
157   - IF isa(uexst) THEN !fit_rchi2_weight.uext=1.
158   -ENDIF
  172 + if isa(st_sed_show) then begin ; because sometimes this keyword isn't set.
  173 +
  174 + teststks = isa(qsed) and isa(used)
  175 +
  176 + IF (not teststks and not isa(polsed)) or (not teststks and not isa(polfrac)) then message, 'Stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
  177 +
  178 + if teststks and isa(polsed) and not (isa(polfrac)) and not(isa(psi_em)) then nowshowing = 'QSED, USED and POLSED'
  179 + if teststks and isa(polsed) and not(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED, POLSED and PSI_EM'
  180 + if teststks and not(isa(polsed)) and isa(polfrac) and not(isa(psi_em)) then nowshowing = 'QSED, USED and POLFRAC'
  181 + if teststks and not(isa(polsed)) and isa(polfrac) and isa(psi_em) then nowshowing = 'QSED, USED, POLFRAC and PSI_EM'
  182 + if teststks and not(isa(polsed)) and not(isa(polfrac)) and not(isa(psi_em)) then nowshowing = 'QSED and USED'
  183 + if teststks and isa(polsed) and isa(polfrac) and not(isa(psi_em)) then nowshowing = 'QSED, USED, POLSED and POLFRAC'
  184 + if teststks and isa(polsed) and isa(polfrac) and isa(psi_em) then nowshowing = 'QSED, USED, POLFRAC, POLSED and PSI_EM'
  185 + if teststks and not(isa(polsed)) and not(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED and PSI_EM'
  186 +
159 187  
  188 + message, 'Polarization component(s) to show: '+nowshowing ,/continue
  189 +
  190 +
  191 + If isa(polsed) then BEGIN ;Apply changes to !dustem_data here.
  192 + !dustem_show.polsed = ptr_new(polsed)
  193 + if keyword_set(f_HI) then if f_HI gt 0 then $
  194 + (*!dustem_show.polsed).values *= f_HI
  195 + ENDIF
  196 + If isa(polfrac) then BEGIN
  197 + !dustem_show.polfrac = ptr_new(polfrac)
  198 + if keyword_set(f_HI) then if f_HI gt 0 then $
  199 + (*!dustem_show.polfrac).values *= f_HI
  200 + ENDIF
  201 + IF isa(qsed) THEN BEGIN
  202 + !dustem_show.qsed = ptr_new(qsed)
  203 + if keyword_set(f_HI) then if f_HI gt 0 then $
  204 + (*!dustem_show.qsed).values *= f_HI
  205 + ENDIF
  206 + if isa(used) then BEGIN
  207 + !dustem_show.used = ptr_new(used)
  208 + if keyword_set(f_HI) then if f_HI gt 0 then $
  209 + (*!dustem_show.used).values *= f_HI
  210 + ENDIF
  211 + if isa(psi_em) then BEGIN
  212 + !dustem_show.psi_em = ptr_new(psi_em)
  213 + if keyword_set(f_HI) then if f_HI gt 0 then $
  214 + (*!dustem_show.psi_em).values *= f_HI
  215 + ENDIF
  216 + endif
  217 + if isa(st_ext_show) then begin
  218 + teststks = isa(qext) and isa(uext)
  219 +
  220 + IF (not teststks and not isa(polext)) then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...'
  221 +
  222 + if teststks and isa(polext) and not(isa(psi_ext)) then nowshowing = 'QEXT, UEXT and POLEXT'
  223 + if teststks and isa(polext) and isa(psi_ext) then nowshowing = 'QEXT, UEXT, POLEXT and PSI_EXT'
  224 + if teststks and not(isa(polext)) and isa(psi_ext) then nowshowing = 'QEXT,UEXT and PSI_EXT'
  225 + if teststks and not(isa(polext)) and not(isa(psi_ext)) then nowshowing = 'QEXT and UEXT'
  226 +
  227 + message, 'Polarization component(s) to show: '+nowshowing ,/continue
  228 +
  229 + ;EXTINCTION
  230 + if isa(polext) then BEGIN
  231 + !dustem_show.polext = ptr_new(polext)
  232 + if keyword_set(f_HI) then if f_HI gt 0 then $
  233 + (*!dustem_show.polext).values *= f_HI
  234 + ENDIF
  235 + IF isa(qext) THEN BEGIN
  236 + !dustem_show.qext = ptr_new(qext)
  237 + if keyword_set(f_HI) then if f_HI gt 0 then $
  238 + (*!dustem_show.qext).values *= f_HI
  239 + ENDIF
  240 + if isa(uext) then BEGIN
  241 + !dustem_show.uext = ptr_new(uext)
  242 + if keyword_set(f_HI) then if f_HI gt 0 then $
  243 + (*!dustem_show.uext).values *= f_HI
  244 + ENDIF
  245 + if isa(psi_ext) then BEGIN
  246 + !dustem_show.psi_ext = ptr_new(psi_ext)
  247 + if keyword_set(f_HI) then if f_HI gt 0 then $
  248 + (*!dustem_show.psi_ext).values *= f_HI
  249 + ENDIF
  250 + endif
  251 +endif
160 252  
161 253 IF keyword_set(rchi2_weight) THEN BEGIN
162 254 !fit_rchi2_weight.sed = rchi2_weight.sed
163 255 !fit_rchi2_weight.ext = rchi2_weight.ext
164 256  
165 257 IF !run_pol THEN BEGIN
166   - IF isa(pst) THEN !fit_rchi2_weight.polsed = rchi2_weight.polsed
167   - IF isa(fpst) THEN !fit_rchi2_weight.polfrac = rchi2_weight.polfrac
168   - IF isa(qst) THEN !fit_rchi2_weight.qsed = rchi2_weight.qsed
169   - IF isa(ust) THEN !fit_rchi2_weight.used = rchi2_weight.used
  258 + ;IF isa(polsed) THEN !fit_rchi2_weight.polsed = rchi2_weight.polsed
  259 + ;IF isa(polfrac) THEN !fit_rchi2_weight.polfrac = rchi2_weight.polfrac
  260 + IF isa(qsed) THEN !fit_rchi2_weight.qsed = rchi2_weight.qsed
  261 + IF isa(used) THEN !fit_rchi2_weight.used = rchi2_weight.used
170 262  
171   - IF isa(pexst) THEN !fit_rchi2_weight.polext = rchi2_weight.polext
172   - IF isa(qexst) THEN !fit_rchi2_weight.qext = rchi2_weight.qext
173   - IF isa(uexst) THEN !fit_rchi2_weight.uext = rchi2_weight.uext
  263 + IF isa(polext) THEN !fit_rchi2_weight.polext = rchi2_weight.polext
  264 + IF isa(qext) THEN !fit_rchi2_weight.qext = rchi2_weight.qext
  265 + IF isa(uext) THEN !fit_rchi2_weight.uext = rchi2_weight.uext
174 266 ENDIF
175 267  
176   -ENDIF
  268 +ENDIF ELSE BEGIN
  269 +
  270 + !fit_rchi2_weight.sed=1.
  271 + !fit_rchi2_weight.ext=1.
  272 +
  273 + If !run_pol THEN BEGIN
  274 + ;IF isa(polsed) THEN !fit_rchi2_weight.polsed=1.
  275 + ;IF isa(polfrac) THEN !fit_rchi2_weight.polfrac=1.
  276 + IF isa(qsed) THEN !fit_rchi2_weight.qsed=1.
  277 + IF isa(used) THEN !fit_rchi2_weight.used=1.
  278 + IF isa(polext) THEN !fit_rchi2_weight.polext=1.
  279 + IF isa(qext) THEN !fit_rchi2_weight.qext=1.
  280 + IF isa(uext) THEN !fit_rchi2_weight.uext=1.
  281 + ENDIF
177 282  
178   -return,ptr_new(obs_st)
  283 +ENDELSE
179 284  
180 285  
181 286 the_end:
182 287  
  288 +;stop
  289 +return,0.
  290 +
  291 +
  292 +
  293 +
  294 +
183 295 END
... ...