OV-GSO-DC / dustem-wrapper_idl

Commit 4b745987f0e1bf4d8254497bfd69539ec33441bd

Authored by Ilyes Choubani
1 parent 96a72c8b
Exists in master

minor update

Showing 6 changed files with 81 additions and 39 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_compute_polext.pro
  Diff comments   View file @4b74598
... ... @@ -13,14 +13,12 @@ ENDIF
13 13  
14 14 POLEXT_spec = sti.polext.ext_tot * (*!dustem_HCD)/1.0e21 ;
15 15  
16   -;Hard-coded test to set to zero potential negative values in the polarized extinction arrays:
  16 +;Hard-coded test to set to potential zero negative values in the polarized extinction arrays:
17 17 ;so far this has been seen in absorption arrays
18 18 ;This test needs to include all grain species that polarize
19 19  
20 20 IF !run_pol THEN BEGIN
21 21  
22   - ;THIS FALSE CORRECT IT USING POLEXT_spec
23   - ;POLEXT_spec.abs_grain[2,*]
24 22 ind_ngtv_plxt = where (POLEXT_spec LT 0, ct_ngtv_plxt)
25 23  
26 24 IF ct_ngtv_plxt NE 0 THEN (POLEXT_spec)[ind_ngtv_plxt] = 0.
... ...
src/idl/dustem_compute_stokes.pro
  Diff comments   View file @4b74598
... ... @@ -104,6 +104,9 @@ ENDIF
104 104 ;Hard-coded block so that atan(U,Q) does not return wrong values when Q eq -0
105 105 testnzero = where(Q_spec EQ -0., ct_nzero)
106 106 if ct_nzero ne 0 then Q_spec[testnzero] = 0
  107 +testnzero = where(U_spec EQ -0., ct_nzero)
  108 +if ct_nzero ne 0 then U_spec[testnzero] = 0
  109 +
107 110  
108 111 PSI_spec = 0.5*atan(U_spec,Q_spec)/!dtor
109 112 ;polar_iqu2ippsi,stI,Q_spec,U_spec,p_bidon,PSI_spec ;same but this line runs several commands for each iteration of this current procedure
... ...
src/idl/dustem_init.pro
  Diff comments   View file @4b74598
... ... @@ -243,7 +243,9 @@ defsysv, '!dustem_mlog', 0 ;necessary for the plotting of negative values o
243 243 defsysv, '!dustem_psi', 0.
244 244 defsysv, '!dustem_psi_ext', 0.
245 245  
  246 +;useful pointers (current dustem output structure and current mpfit params values)
246 247 defsysv,'!dustem_current',ptr_new()
  248 +defsysv,'!dustem_current_params',ptr_new()
247 249  
248 250 ;IDL> help,*!dustem_data,/str
249 251 ;** Structure <2173b20>, 4 tags, length=448, data length=448, refs=1:
... ...
src/idl/dustem_init_parinfo.pro
  Diff comments   View file @4b74598
... ... @@ -6,6 +6,8 @@ PRO dustem_init_parinfo,pd,iv, $
6 6 (*!dustem_fit).param_descs=ptr_new(pd)
7 7  
8 8 (*!dustem_fit).param_init_values=ptr_new(iv)
  9 +
  10 +!dustem_current_params = ptr_new(iv/iv)
9 11 ;stop
10 12 dustem_set_func_ind,pd,iv
11 13  
... ...
src/idl/dustemcgwin_dataset.pro
  Diff comments   View file @4b74598
... ... @@ -68,7 +68,10 @@ if keyword_set(dataset) then begin
68 68  
69 69 'SED': begin
70 70  
71   - IF ct_m NE 0 then xr=xr_m & yr=yr_m
  71 + IF ct_m NE 0 then begin
  72 + xr=xr_m
  73 + yr=yr_m
  74 + endif
72 75  
73 76 vectw=DINDGEN(1001)*(alog10(1E7) - alog10(1E5))/(1000 - 1L) + alog10(1E5)
74 77 vectw=10^vectw[1:*]
... ... @@ -271,8 +274,12 @@ if keyword_set(dataset) then begin
271 274 ; xr=[0.01,30]
272 275 ; yr=[1.00E-06,10]
273 276  
274   - IF ct_x NE 0 then xr=xr_x & yr=yr_x
  277 + IF ct_x NE 0 then begin
275 278  
  279 + xr=xr_x
  280 + yr=yr_x
  281 +
  282 + ENDIF
276 283 vectw=DINDGEN(1001)*(alog10(1E7) - alog10(1E5))/(1000 - 1L) + alog10(1E5)
277 284 vectw=10^vectw[1:*]
278 285 vectw=[wavs,vectw]
... ... @@ -478,7 +485,10 @@ if keyword_set(dataset) then begin
478 485  
479 486 'POLEXT': begin
480 487  
481   - IF ct_x NE 0 then xr=xr_x & yr=yr_x
  488 + IF ct_x NE 0 then begin
  489 + xr=xr_x
  490 + yr=yr_x
  491 + ENDIF
482 492  
483 493 vectw=DINDGEN(1001)*(alog10(1E7) - alog10(1E5))/(1000 - 1L) + alog10(1E5)
484 494 vectw=10^vectw[1:*]
... ... @@ -675,8 +685,11 @@ if keyword_set(dataset) then begin
675 685  
676 686 'POLSED': begin
677 687  
678   - IF ct_m NE 0 then xr=xr_m & yr=yr_m
679   -
  688 + IF ct_m NE 0 then begin
  689 + xr=xr_m
  690 + yr=yr_m
  691 +
  692 + ENDIF
680 693 vectw=DINDGEN(1001)*(alog10(1E7) - alog10(1E5))/(1000 - 1L) + alog10(1E5)
681 694 vectw=10^vectw[1:*]
682 695 vectw=[wavs,vectw]
... ... @@ -864,8 +877,11 @@ if keyword_set(dataset) then begin
864 877  
865 878 'POLFRAC': begin
866 879  
867   - IF ct_m NE 0 then xr=xr_m & yr=yr_m
868   -
  880 + IF ct_m NE 0 then begin
  881 + xr=xr_m
  882 + yr=yr_m
  883 +
  884 + ENDIF
869 885 idx_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',ct_filt)
870 886 idx_spec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',ct_spec)
871 887  
... ... @@ -1009,8 +1025,12 @@ if keyword_set(dataset) then begin
1009 1025  
1010 1026 'FPOLEXT': BEGIN ; For completeness (not sure of the physical usefulness of this qauntity )
1011 1027  
1012   - IF ct_x NE 0 then xr=xr_x & yr=yr_x
  1028 + IF ct_x NE 0 then begin
1013 1029  
  1030 + xr=xr_x
  1031 + yr=yr_x
  1032 +
  1033 + ENDIF
1014 1034 idx_filt=where((*!dustem_data.fpolext).filt_names NE 'SPECTRUM',ct_filt)
1015 1035 idx_spec=where((*!dustem_data.fpolext).filt_names EQ 'SPECTRUM',ct_spec)
1016 1036  
... ... @@ -1159,8 +1179,12 @@ if keyword_set(dataset) then begin
1159 1179  
1160 1180 'PSI_EM': begin ;will be copying code from qsed and hopefully it will work. NB: THIS IS FALSE . YOU DO NOT NEED A LOGARITHMIC AXIS
1161 1181  
1162   - IF ct_m NE 0 then xr=xr_m & yr=yr_m
  1182 + IF ct_m NE 0 then begin
1163 1183  
  1184 + xr=xr_m
  1185 + yr=yr_m
  1186 +
  1187 + ENDIF
1164 1188 idx_filt=where((*!dustem_data.psi_em).filt_names NE 'SPECTRUM',ct_filt)
1165 1189 idx_spec=where((*!dustem_data.psi_em).filt_names EQ 'SPECTRUM',ct_spec)
1166 1190  
... ... @@ -1295,25 +1319,10 @@ if keyword_set(dataset) then begin
1295 1319 endif
1296 1320  
1297 1321 endif
1298   -
1299   -
1300   -
1301   - ;stop
1302   -
1303   -
1304   -
1305   -
1306   -
1307   -
1308 1322  
1309 1323 endelse
1310   -
1311   -
1312   -
1313   - ;endelse
1314   -
  1324 +
1315 1325 endelse
1316   - ;stop
1317 1326  
1318 1327  
1319 1328 end
... ... @@ -1321,7 +1330,10 @@ if keyword_set(dataset) then begin
1321 1330 'PSI_EXT': begin ;MODIFYING THIS NOW - You should write a plugin for the modification of polarized extinction so that !dustem_psi_ext is different than !dustem_psi.
1322 1331  
1323 1332  
1324   - IF ct_x NE 0 then xr=xr_x & yr=yr_x
  1333 + IF ct_x NE 0 then begin
  1334 + xr=xr_x
  1335 + yr=yr_x
  1336 + ENDIF
1325 1337  
1326 1338 idx_filt=where((*!dustem_data.psi_ext).filt_names NE 'SPECTRUM',ct_filt)
1327 1339 idx_spec=where((*!dustem_data.psi_ext).filt_names EQ 'SPECTRUM',ct_spec)
... ... @@ -1341,7 +1353,9 @@ if keyword_set(dataset) then begin
1341 1353 ;##################TAKEN FROM SED
1342 1354  
1343 1355 ;Plotting of the spectra of the dust species
  1356 +
1344 1357 FOR i=0L,Ngrains-1 DO BEGIN
  1358 +
1345 1359 vecto = transpose((st.polext.(1))[i,*]+(st.polext.(2))[i,*])
1346 1360 vecte = transpose(((st.ext.(1))[i,*]+(st.ext.(2))[i,*]))
1347 1361 testneq = where(vecto ne 0,countneq)
... ... @@ -1349,7 +1363,7 @@ if keyword_set(dataset) then begin
1349 1363 ;testing whether the polext values contain negative ones
1350 1364 testpos = where(vecto LT 0, countpos)
1351 1365 IF countpos NE 0 THEN vecto[testpos] = 0
1352   -
  1366 +
1353 1367 if countneq ne 0 then begin
1354 1368  
1355 1369 vecfin = vecto*0.0
... ... @@ -1467,7 +1481,11 @@ if keyword_set(dataset) then begin
1467 1481  
1468 1482 'QSED': begin ; THE ONE I WILL BE WORKING ON
1469 1483  
1470   - IF ct_m NE 0 then xr=xr_m & yr=yr_m
  1484 + IF ct_m NE 0 then begin
  1485 +
  1486 + xr=xr_m
  1487 + yr=yr_m
  1488 + ENDIF
1471 1489  
1472 1490 vectw=DINDGEN(1001)*(alog10(1E7) - alog10(1E5))/(1000 - 1L) + alog10(1E5)
1473 1491 vectw=10^vectw[1:*]
... ... @@ -1821,7 +1839,12 @@ if keyword_set(dataset) then begin
1821 1839  
1822 1840 'USED': begin
1823 1841  
1824   - IF ct_m NE 0 then xr=xr_m & yr=yr_m
  1842 + IF ct_m NE 0 then begin
  1843 +
  1844 + xr=xr_m
  1845 + yr=yr_m
  1846 +
  1847 + ENDIF
1825 1848  
1826 1849 vectw=DINDGEN(1001)*(alog10(1E7) - alog10(1E5))/(1000 - 1L) + alog10(1E5)
1827 1850 vectw=10^vectw[1:*]
... ... @@ -2176,7 +2199,13 @@ if keyword_set(dataset) then begin
2176 2199  
2177 2200 'QEXT': begin ;Correct this like you corrected extinction...
2178 2201  
2179   - IF ct_x NE 0 then xr=xr_x & yr=yr_x
  2202 + IF ct_x NE 0 then begin
  2203 +
  2204 + xr=xr_x
  2205 + yr=yr_x
  2206 +
  2207 + ENDIF
  2208 +
2180 2209 vectw=DINDGEN(1001)*(alog10(1E7) - alog10(1E5))/(1000 - 1L) + alog10(1E5)
2181 2210 vectw=10^vectw[1:*]
2182 2211 vectw=[wavs,vectw]
... ... @@ -2270,18 +2299,18 @@ if keyword_set(dataset) then begin
2270 2299 endif else begin ; normal plot
2271 2300  
2272 2301 if keyword_set(refresh) then begin ;The data points in the plot are being refreshed
2273   - ;stop
  2302 +
2274 2303 if keyword_set(positive_only) then begin
2275 2304  
2276 2305 if testpstv ne 0 then begin
2277   - ;stop
  2306 +
2278 2307 FOR i=0L,Ngrains-1 DO BEGIN
2279 2308  
2280 2309 vecto = transpose((st.polext.(1))[i,*]+(st.polext.(2))[i,*])
2281 2310 vecte = transpose((st.ext.(1))[i,*]+(st.ext.(2))[i,*])
2282 2311 testneq = where( vecto ne 0,countneq)
2283 2312  
2284   - ;stop
  2313 +
2285 2314 if countneq ne 0 then begin
2286 2315 vecfin = vecto*0.0 ;em_tot can be used here right?
2287 2316 ;stop
... ... @@ -2296,7 +2325,8 @@ if keyword_set(dataset) then begin
2296 2325 dustem_plot_mlog,1/st.polext.wav,specqgrain,ppositions=position,/positive_only,/overplot,noerase=1,color=use_cols[i],xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
2297 2326  
2298 2327 endif
2299   - ;stop
  2328 +
  2329 +
2300 2330 ENDFOR
2301 2331  
2302 2332 for i=0L,n_plgns-1 do begin
... ... @@ -2349,7 +2379,7 @@ if keyword_set(dataset) then begin
2349 2379 vecfin[testneq] = (vecto[testneq])/vecte[testneq];extra_spec[testneq];(st.sed.em_tot)[testneq];(st.sed.(i+1))[testneq];*fact;(st.sed.(i+1))[testneq]
2350 2380 ;stop
2351 2381 polar_ippsi2iqu,vecte*(*!dustem_HCD)/1.0e21,specqgrain,specugrain,vecfin,replicate(!dustem_psi_ext,n_elements(vecfin)) ;temporary solution
2352   - ;stop
  2382 + stop
2353 2383 ;dustem_plot_mlog,st.polsed.wav,specqgrain,ppositions=position,/overplot,noerase=1,color=use_cols[i],xr=xr,/xlog
2354 2384 dustem_plot_mlog,1/st.polext.wav,specqgrain,ppositions=position,/negative_only,/overplot,noerase=1,color=use_cols[i],xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
2355 2385 endif
... ... @@ -2552,7 +2582,13 @@ if keyword_set(dataset) then begin
2552 2582 end
2553 2583  
2554 2584 'UEXT': begin
2555   - IF ct_x NE 0 then xr=xr_x & yr=yr_x
  2585 + IF ct_x NE 0 then begin
  2586 +
  2587 + xr=xr_x
  2588 + yr=yr_x
  2589 +
  2590 + ENDIF
  2591 +
2556 2592 vectw=DINDGEN(1001)*(alog10(1E7) - alog10(1E5))/(1000 - 1L) + alog10(1E5)
2557 2593 vectw=10^vectw[1:*]
2558 2594 vectw=[wavs,vectw]
... ...
src/idl/from_cmtotal/dustem_mpfit.pro
  Diff comments   View file @4b74598
... ... @@ -3561,6 +3561,7 @@ function dustem_mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $
3561 3561 fnorm = fnorm1
3562 3562 iter = iter + 1
3563 3563 !dustem_iter.act=iter
  3564 + !dustem_current_params = ptr_new(x)
3564 3565 endif
3565 3566  
3566 3567 ;; Tests for convergence
... ...