OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_compute_stokes.pro 4.23 KB
Edit Raw Normal View History
4fd64cbb   Ilyes Choubani   dustem_fit_sed_po... 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
 
FUNCTION dustem_compute_stokes,p_dim,stp,dustem_qsed=dustem_qsed,dustem_used=dustem_used,Q_spec=Q_spec,U_spec=U_spec,out_st=out_st,_extra=extra

;+
; NAME:
;    dustem_compute_stokes
; PURPOSE:
;    Provides Stokes Q and Stokes U dustem spectra and wrapper SEDs 
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    out=dustem_compute_polsed(p_dim[,st=][,_extra=][,/help])
; INPUTS:
;    p_dim      = parameter values
; OPTIONAL INPUT PARAMETERS:
;    stp        = Dustem structure
; OUTPUTS:
;    sed         : computed SED for filters in !dustem_data
;    dustem_qsed : SED in Q
;    dustem_used : SED in U
; OPTIONAL OUTPUT PARAMETERS:
;    stp        = Dustem output structure
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.

IF keyword_set(help) THEN BEGIN
    doc_library,'dustem_compute_stokes'
    out=0.
    goto,the_end
ENDIF

IF not keyword_set(stp) THEN BEGIN
    stp=dustem_run(p_dim)
    dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),stp
ENDIF


If keyword_set(result) THEN result=stp

fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(stp.polsed).wav)*1.e20/1.e7
stI = stp.sed.em_tot * fact   ;This is Total intensity I
P = stp.polsed.em_tot * fact   ;This is Polarized intensity P
out=0.
Nwaves=(size(stI))[1]

frac=P/stI
tes=where(finite(frac) eq 0)
frac(tes)=0.


FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN   
    
    IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'MODIFY_DUST_POLAR_ANGLE') THEN BEGIN
        
        Q_spec=(*(*!dustem_plugin).(i))[*,0]
        U_spec=(*(*!dustem_plugin).(i))[*,1]
        
    ENDIF 
ENDFOR    
    
IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0.,Nwaves) 
        
FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN

    IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_QSED') THEN BEGIN
        Q_spec+=(*(*!dustem_plugin).(i))[*,1]
    ENDIF
    
    IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_USED') THEN BEGIN
        U_spec+=(*(*!dustem_plugin).(i))[*,2]
    ENDIF              
    
ENDFOR


;-----------------------------------------

;largeP=sqrt(StokesQ^2+StokesU^2)

;INITIALIZING THE STOKES SEDS
dustem_qsed = (*!dustem_data.sed).values * 0.
dustem_used = dustem_qsed


if not isarray(stI) THEN stop 


;Performing color corrections on the dustem spectra and generating the SEDs at the given filters 

IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
    message,'DOING color correction calculations',/info
ENDIF ELSE BEGIN
          message,'SKIPPING color correction calculations',/info
      ENDELSE

ind_polsed=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_polsed)
  
IF count_polsed NE 0 THEN BEGIN ; INCLUDE COLOR CORRECTIONS OF (Q,U) for dust & synch
  filter_names=((*!dustem_data.polsed).filt_names)(ind_polsed)
      
  if isa(dustem_qsed) then begin 
    sqsed=dustem_cc(stp.sed.wav,Q_spec,filter_names,cc=cc)
    dustem_qsed(ind_polsed)=sqsed
  endif
  if isa(dustem_used) then begin
    sused=dustem_cc(stp.sed.wav,U_spec,filter_names,cc=cc)
    dustem_used(ind_polsed)=sused
  endif
  ;--------------------------------------------------------------------------
  
ENDIF

;For spectrum data points, interpolate in log-log
;Linear interpolation leads to wrong values, in particular where few
;wavelengths points exist in the model (long wavelengths).


ind_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',count_spec)
IF count_spec NE 0 THEN BEGIN
  if isa(dustem_qsed) then dustem_qsed(ind_spec)=interpol(Q_spec,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))  
  if isa(dustem_used) then dustem_used(ind_spec)=interpol(U_spec,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
ENDIF
out_st=stp

;out=fltarr(Nwaves)

; IF keyword_set(dustem_qsed) THEN out=dustem_qsed
; IF keyword_set(dustem_used) THEN out=dustem_used
; IF keyword_set(Q_spec) THEN out=Q_spec
; IF keyword_set(U_spec) THEN out=U_spec



;clean pointers
heap_gc

the_end:

RETURN,out

END