dustem_fillup_systvar_from_fits.pro
8.03 KB
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
PRO dustem_fillup_systvar_from_fits,syst_var,str_input_SED,used_pol
;IC: some way for the procedure to be aware of emission and extinction
test_m = fix(total(tag_names (str_input_sed) EQ 'STOKESI')) ;a proper boolean
Nsed=n_elements(str_input_SED)
sed=dustem_initialize_internal_sed(Nsed)
IF test_m THEN BEGIN
(syst_var).sed=ptr_new(sed)
(*syst_var.sed).instru_names=dustem_filter2instru(str_input_SED.filter)
(*syst_var.sed).filt_names=strtrim(str_input_SED.filter,2)
(*syst_var.sed).wav=str_input_SED.wavelength
(*syst_var.sed).values=str_input_SED.STOKESI
(*syst_var.sed).sigma=la_power(str_input_SED.varianceII,0.5)
ENDIF ELSE BEGIN
(syst_var).ext=ptr_new(sed)
(*syst_var.ext).instru_names=dustem_filter2instru(str_input_SED.filter)
(*syst_var.ext).filt_names=strtrim(str_input_SED.filter,2)
(*syst_var.ext).wav=str_input_SED.wavelength
(*syst_var.ext).values=str_input_SED.STOKESI_EXT
(*syst_var.ext).sigma=la_power(str_input_SED.varianceII_EXT,0.5)
ENDELSE
IF used_pol THEN BEGIN
;stop
IF test_m THEN BEGIN;IC: quickest way I found to make this work
;===== note: those should include only points with polarization data, not all sed points
indq=where(str_input_SED.STOKESQ NE la_undef(),Nqsed)
indu=where(str_input_SED.STOKESU NE la_undef(),Nused)
qsed=dustem_initialize_internal_sed(Nqsed)
used=dustem_initialize_internal_sed(Nused)
indpolar=where(str_input_SED.STOKESU NE la_undef() AND str_input_SED.STOKESQ NE la_undef(),Npolar)
psi_sed=dustem_initialize_internal_sed(Npolar)
smallp_sed=dustem_initialize_internal_sed(Npolar)
largep_sed=dustem_initialize_internal_sed(Npolar)
syst_var.qsed=ptr_new(qsed)
syst_var.used=ptr_new(used)
(*syst_var.qsed).instru_names=dustem_filter2instru(str_input_SED[indq].filter)
(*syst_var.qsed).filt_names=strtrim(str_input_SED[indq].filter,2)
(*syst_var.qsed).wav=str_input_SED[indq].wavelength
(*syst_var.qsed).values=str_input_SED[indq].STOKESQ
(*syst_var.qsed).sigma=la_power(str_input_SED[indq].varianceQQ,0.5)
;now U
(*syst_var.used).instru_names=dustem_filter2instru(str_input_SED[indu].filter)
(*syst_var.used).filt_names=strtrim(str_input_SED[indu].filter,2)
(*syst_var.used).wav=str_input_SED[indu].wavelength
(*syst_var.used).values=str_input_SED[indu].STOKESU
(*syst_var.used).sigma=la_power(str_input_SED[indu].varianceUU,0.5)
;=== add up dependent fields
syst_var.psi_em=ptr_new(psi_sed)
syst_var.polfrac=ptr_new(smallp_sed)
syst_var.polsed=ptr_new(largep_sed)
polar_iqu2ippsi,str_input_SED[indpolar].STOKESI,str_input_SED[indpolar].STOKESQ,str_input_SED[indpolar].STOKESU,p,psi,largep=largep
(*syst_var.psi_em).values=psi
(*syst_var.polfrac).values=p
(*syst_var.polsed).values=largep
CIQ=str_input_SED[indpolar].VARIANCEII & CIQ[*]=0. ;because input SEDs don't include covariances, yet
CIU=str_input_SED[indpolar].VARIANCEII & CIU[*]=0. ;because input SEDs don't include covariances, yet
CQU=str_input_SED[indpolar].VARIANCEII & CQU[*]=0. ;because input SEDs don't include covariances, yet
polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI $
,str_input_SED[indpolar].STOKESQ $
,str_input_SED[indpolar].STOKESU $
,str_input_SED[indpolar].VARIANCEII $
,str_input_SED[indpolar].VARIANCEQQ $
,str_input_SED[indpolar].VARIANCEUU $
,CIQ,CIU,CQU,variance_smallp,variance_psi,variance_largep,/lac,snr_largeP=snr_largeP
(*syst_var.psi_em).instru_names=dustem_filter2instru(str_input_SED[indpolar].filter)
(*syst_var.psi_em).filt_names=strtrim(str_input_SED[indpolar].filter,2)
(*syst_var.psi_em).wav=str_input_SED[indpolar].wavelength
(*syst_var.polfrac).instru_names = (*syst_var.psi_em).instru_names
(*syst_var.polfrac).filt_names = (*syst_var.psi_em).filt_names
(*syst_var.polfrac).wav = (*syst_var.psi_em).wav
(*syst_var.polsed).instru_names = (*syst_var.psi_em).instru_names
(*syst_var.polsed).filt_names = (*syst_var.psi_em).filt_names
(*syst_var.polsed).wav = (*syst_var.psi_em).wav
(*syst_var.psi_em).sigma=la_power(variance_psi,0.5)
(*syst_var.polfrac).sigma=la_power(variance_smallp,0.5)
(*syst_var.polsed).sigma=la_power(variance_largep,0.5)
ENDIF ELSE BEGIN
;===== note: those should include only points with polarization data, not all sed points
indq=where(str_input_SED.STOKESQ_EXT NE la_undef(),Nqsed)
indu=where(str_input_SED.STOKESU_EXT NE la_undef(),Nused)
qsed=dustem_initialize_internal_sed(Nqsed)
used=dustem_initialize_internal_sed(Nused)
indpolar=where(str_input_SED.STOKESU_EXT NE la_undef() AND str_input_SED.STOKESQ_EXT NE la_undef(),Npolar)
psi_sed=dustem_initialize_internal_sed(Npolar)
smallp_sed=dustem_initialize_internal_sed(Npolar)
largep_sed=dustem_initialize_internal_sed(Npolar)
syst_var.qext=ptr_new(qsed)
syst_var.uext=ptr_new(used)
(*syst_var.qext).instru_names=dustem_filter2instru(str_input_SED[indq].filter)
(*syst_var.qext).filt_names=strtrim(str_input_SED[indq].filter,2)
(*syst_var.qext).wav=str_input_SED[indq].wavelength
(*syst_var.qext).values=str_input_SED[indq].STOKESQ_EXT
(*syst_var.qext).sigma=la_power(str_input_SED[indq].varianceQQ_EXT,0.5)
;now U
(*syst_var.uext).instru_names=dustem_filter2instru(str_input_SED[indu].filter)
(*syst_var.uext).filt_names=strtrim(str_input_SED[indu].filter,2)
(*syst_var.uext).wav=str_input_SED[indu].wavelength
(*syst_var.uext).values=str_input_SED[indu].STOKESU_EXT
(*syst_var.uext).sigma=la_power(str_input_SED[indu].varianceUU_EXT,0.5)
;=== add up dependent fields
syst_var.psi_ext=ptr_new(psi_sed)
syst_var.fpolext=ptr_new(smallp_sed)
syst_var.polext=ptr_new(largep_sed)
polar_iqu2ippsi,str_input_SED[indpolar].STOKESI_EXT,str_input_SED[indpolar].STOKESQ_EXT,str_input_SED[indpolar].STOKESU_EXT,p,psi,largep=largep
(*syst_var.psi_ext).values=psi
(*syst_var.fpolext).values=p
(*syst_var.polext).values=largep
CIQ_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CIQ_EXT[*]=0. ;because input SEDs don't include covariances, yet
CIU_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CIU_EXT[*]=0. ;because input SEDs don't include covariances, yet
CQU_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CQU_EXT[*]=0. ;because input SEDs don't include covariances, yet
polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI_EXT $
,str_input_SED[indpolar].STOKESQ_EXT $
,str_input_SED[indpolar].STOKESU_EXT $
,str_input_SED[indpolar].VARIANCEII_EXT $
,str_input_SED[indpolar].VARIANCEQQ_EXT $
,str_input_SED[indpolar].VARIANCEUU_EXT $
,CIQ_EXT,CIU_EXT,CQU_EXT,variance_smallp_EXT,variance_psi_EXT,variance_largep_EXT,/lac,snr_largeP=snr_largeP_EXT
(*syst_var.psi_ext).instru_names=dustem_filter2instru(str_input_SED[indpolar].filter)
(*syst_var.psi_ext).filt_names=strtrim(str_input_SED[indpolar].filter,2)
(*syst_var.psi_ext).wav=str_input_SED[indpolar].wavelength
(*syst_var.fpolext).instru_names = (*syst_var.psi_ext).instru_names
(*syst_var.fpolext).filt_names = (*syst_var.psi_ext).filt_names
(*syst_var.fpolext).wav = (*syst_var.psi_ext).wav
(*syst_var.polext).instru_names = (*syst_var.psi_ext).instru_names
(*syst_var.polext).filt_names = (*syst_var.psi_ext).filt_names
(*syst_var.polext).wav = (*syst_var.psi_ext).wav
(*syst_var.psi_ext).sigma=la_power(variance_psi_EXT,0.5)
(*syst_var.fpolext).sigma=la_power(variance_smallp_EXT,0.5)
(*syst_var.polext).sigma=la_power(variance_largep_EXT,0.5)
ENDELSE
ENDIF
END