-
Notifications
You must be signed in to change notification settings - Fork 0
/
cc_calc_quijote.py
419 lines (386 loc) · 22.3 KB
/
cc_calc_quijote.py
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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Calculate colour corrections for QUIJOTE
# Version history:
#
# 16-Jun-2019 M. Peel Started
# 17-Jun-2019 M. Peel Tidied up
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
from astrocode.astroutils import *
from fastcc import *
from cc_calc_functions import *
outdir = 'plots_2021_01_15/'
print(outdir)
ensure_dir(outdir)
# alphas = [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0]
alphas = [-0.31, 2.0]
alphas = np.asarray(alphas)
# QUIJOTE MFI
pol1 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/Pol1_bandpass.dat')
plot_bandpass([pol1[0], pol1[1]],outdir+'mfi_pol1_ch1.png')
plot_bandpass([pol1[0], pol1[2]],outdir+'mfi_pol1_ch2.png')
plot_bandpass([pol1[0], pol1[3]],outdir+'mfi_pol1_ch3.png')
plot_bandpass([pol1[0], pol1[4]],outdir+'mfi_pol1_ch4.png')
plot_bandpass([pol1[0], pol1[5]],outdir+'mfi_pol1_ch5.png')
plot_bandpass([pol1[0], pol1[6]],outdir+'mfi_pol1_ch6.png')
plot_bandpass([pol1[0], pol1[7]],outdir+'mfi_pol1_ch7.png')
plot_bandpass([pol1[0], pol1[8]],outdir+'mfi_pol1_ch8.png')
# pol2 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/Pol2_bandpass.dat')
# pol2 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020_05_21_Pol2_bandpass.dat')
pol2 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020-10_Horn2_bandpass_intensity.dat')
pol2_pol = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020-10_Horn2_bandpass_polarization.dat')
# pol2 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020_05_22_0deg_Pol2_bandpass.dat')
# pol2 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020_05_22_22p5deg_Pol2_bandpass.dat')
# minlength = np.min([len(pol2a[0]),len(pol2b[0])])
# pol2a = pol2a[:,:minlength]
# pol2b = pol2b[:,:minlength]
# pol2a[pol2a < 0.0] = 0.0
# pol2b[pol2b < 0.0] = 0.0
# pol2 = np.sqrt(pol2a**2+pol2b**2)
# for i in range(0,len(pol2)):
# pol2[i] = pol2[i] / np.sum(pol2[i])
# pol2[~np.isfinite(pol2)] = 0.0
plot_bandpass([pol2[0], pol2[1]],outdir+'mfi_pol2_ch1.png')
plot_bandpass([pol2[0], pol2[2]],outdir+'mfi_pol2_ch2.png')
plot_bandpass([pol2[0], pol2[3]],outdir+'mfi_pol2_ch3.png')
plot_bandpass([pol2[0], pol2[4]],outdir+'mfi_pol2_ch4.png')
plot_bandpass([pol2[0], pol2[5]],outdir+'mfi_pol2_ch5.png')
plot_bandpass([pol2[0], pol2[6]],outdir+'mfi_pol2_ch6.png')
plot_bandpass([pol2[0], pol2[7]],outdir+'mfi_pol2_ch7.png')
plot_bandpass([pol2[0], pol2[8]],outdir+'mfi_pol2_ch8.png')
# pol3_orig = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/Pol3_bandpass.dat')
# pol3 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/Horn3_bandpass_intensity.dat')
# pol3_pol = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/Horn3_bandpass_polarization.dat')
pol3 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020-10_Horn3_bandpass_intensity.dat')
pol3_pol = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020-10_Horn3_bandpass_polarization.dat')
# pol3 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020_05_22_0deg_Pol3_bandpass.dat')
# pol3 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020_05_22_22p5deg_Pol3_bandpass.dat')
# minlength = np.min([len(pol3a[0]),len(pol3b[0])])
# pol3a = pol3a[:,:minlength]
# pol3b = pol3b[:,:minlength]
# pol3 = np.sqrt(pol3a**2+pol3b**2)
plot_bandpass([pol3[0], pol3[1]],outdir+'mfi_pol3_ch1.png')
plot_bandpass([pol3[0], pol3[2]],outdir+'mfi_pol3_ch2.png')
plot_bandpass([pol3[0], pol3[3]],outdir+'mfi_pol3_ch3.png')
plot_bandpass([pol3[0], pol3[4]],outdir+'mfi_pol3_ch4.png')
plot_bandpass([pol3[0], pol3[5]],outdir+'mfi_pol3_ch5.png')
plot_bandpass([pol3[0], pol3[6]],outdir+'mfi_pol3_ch6.png')
plot_bandpass([pol3[0], pol3[7]],outdir+'mfi_pol3_ch7.png')
plot_bandpass([pol3[0], pol3[8]],outdir+'mfi_pol3_ch8.png')
# pol4 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/Pol4_bandpass.dat')
# pol4 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020_05_21_Pol4_bandpass.dat')
pol4 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020-10_Horn4_bandpass_intensity.dat')
pol4_pol = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020-10_Horn4_bandpass_polarization.dat')
# pol4 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020_05_22_0deg_Pol4_bandpass.dat')
# pol4 = read_quijote_mfi_bandpass('/Users/mpeel/Documents/maps/quijote_mfi/2020_05_22_22p5deg_Pol4_bandpass.dat')
# minlength = np.min([len(pol4a[0]),len(pol4b[0])])
# pol4a = pol4a[:,:minlength]
# pol4b = pol4b[:,:minlength]
# pol4 = np.sqrt(pol4a**2+pol4b**2)
for i in range(1,9):
pol4[i][pol4[0]<14] = 0.0
plot_bandpass([pol4[0], pol4[1]],outdir+'mfi_pol4_ch1.png')
plot_bandpass([pol4[0], pol4[2]],outdir+'mfi_pol4_ch2.png')
plot_bandpass([pol4[0], pol4[3]],outdir+'mfi_pol4_ch3.png')
plot_bandpass([pol4[0], pol4[4]],outdir+'mfi_pol4_ch4.png')
plot_bandpass([pol4[0], pol4[5]],outdir+'mfi_pol4_ch5.png')
plot_bandpass([pol4[0], pol4[6]],outdir+'mfi_pol4_ch6.png')
plot_bandpass([pol4[0], pol4[7]],outdir+'mfi_pol4_ch7.png')
plot_bandpass([pol4[0], pol4[8]],outdir+'mfi_pol4_ch8.png')
plot_bandpass_all([[pol1[0], pol1[1]], [pol1[0], pol1[2]], [pol1[0], pol1[3]], [pol1[0], pol1[4]], [pol1[0], pol1[5]],[pol1[0], pol1[6]], [pol1[0], pol1[7]], [pol1[0], pol1[8]], [pol2[0], pol2[1]], [pol2[0], pol2[2]], [pol2[0], pol2[3]], [pol2[0], pol2[4]], [pol2[0], pol2[5]], [pol2[0], pol2[6]], [pol2[0], pol2[7]], [pol2[0], pol2[8]], [pol3[0], pol3[1]], [pol3[0], pol3[2]], [pol3[0], pol3[3]], [pol3[0], pol3[4]], [pol3[0], pol3[5]], [pol3[0], pol3[6]], [pol3[0], pol3[7]], [pol3[0], pol3[8]], [pol4[0], pol4[1]], [pol4[0], pol4[2]], [pol4[0], pol4[3]], [pol4[0], pol4[4]], [pol4[0], pol4[5]], [pol4[0], pol4[6]], [pol4[0], pol4[7]],[pol4[0], pol4[8]]], outdir+'mfi_comb.pdf')
mfi111 = combine_bandpasses([pol1[0], pol1[2]], [pol1[0], pol1[4]],[pol1[0], pol1[6]],[pol1[0], pol1[8]])
plot_bandpass(mfi111,outdir+'mfi_111.png')
mfi113 = combine_bandpasses([pol1[0], pol1[1]], [pol1[0], pol1[3]],[pol1[0], pol1[5]],[pol1[0], pol1[7]])
plot_bandpass(mfi113,outdir+'mfi_113.png')
mfi217 = combine_bandpasses([pol2[0], pol2[2]], [pol2[0], pol2[4]],[pol2[0], pol2[6]],[pol2[0], pol2[8]])
plot_bandpass(mfi217,outdir+'mfi_217.png')
mfi219 = combine_bandpasses([pol2[0], pol2[1]], [pol2[0], pol2[3]],[pol2[0], pol2[5]],[pol2[0], pol2[7]])
plot_bandpass(mfi219,outdir+'mfi_219.png')
mfi217_pol = combine_bandpasses([pol2_pol[0], pol2_pol[2]], [pol2_pol[0], pol2_pol[4]],[pol2_pol[0], pol2_pol[6]],[pol2_pol[0], pol2_pol[8]])
plot_bandpass(mfi217_pol,outdir+'mfi_217_pol.png')
mfi219_pol = combine_bandpasses([pol2_pol[0], pol2_pol[1]], [pol2_pol[0], pol2_pol[3]],[pol2_pol[0], pol2_pol[5]],[pol2_pol[0], pol2_pol[7]])
plot_bandpass(mfi219_pol,outdir+'mfi_219_pol.png')
mfi311 = combine_bandpasses([pol3[0], pol3[2]], [pol3[0], pol3[4]],[pol3[0], pol3[6]],[pol3[0], pol3[8]])
plot_bandpass(mfi311,outdir+'mfi_311.png')
mfi313 = combine_bandpasses([pol3[0], pol3[1]], [pol3[0], pol3[3]],[pol3[0], pol3[5]],[pol3[0], pol3[7]])
plot_bandpass(mfi313,outdir+'mfi_313.png')
mfi311pol = combine_bandpasses([pol3_pol[0], pol3_pol[2]], [pol3_pol[0], pol3_pol[4]],[pol3_pol[0], pol3_pol[6]],[pol3_pol[0], pol3_pol[8]])
plot_bandpass(mfi311pol,outdir+'mfi_311_pol.png')
mfi313pol = combine_bandpasses([pol3_pol[0], pol3_pol[1]], [pol3_pol[0], pol3_pol[3]],[pol3_pol[0], pol3_pol[5]],[pol3_pol[0], pol3_pol[7]])
plot_bandpass(mfi313pol,outdir+'mfi_313_pol.png')
# mfi311orig = combine_bandpasses([pol3_orig[0], pol3_orig[2]], [pol3_orig[0], pol3_orig[4]],[pol3_orig[0], pol3_orig[6]],[pol3_orig[0], pol3_orig[8]])
# plot_bandpass(mfi311orig,outdir+'mfi_311_orig.png')
# mfi313orig = combine_bandpasses([pol3_orig[0], pol3_orig[1]], [pol3_orig[0], pol3_orig[3]],[pol3_orig[0], pol3_orig[5]],[pol3_orig[0], pol3_orig[7]])
# plot_bandpass(mfi313orig,outdir+'mfi_313_orig.png')
mfi417 = combine_bandpasses([pol4[0], pol4[2]], [pol4[0], pol4[4]],[pol4[0], pol4[6]],[pol4[0], pol4[8]])
plot_bandpass(mfi417,outdir+'mfi_417.png')
mfi419 = combine_bandpasses([pol4[0], pol4[1]], [pol4[0], pol4[3]],[pol4[0], pol4[5]],[pol4[0], pol4[7]])
plot_bandpass(mfi419,outdir+'mfi_419.png')
mfi417_pol = combine_bandpasses([pol4_pol[0], pol4_pol[2]], [pol4_pol[0], pol4_pol[4]],[pol4_pol[0], pol4_pol[6]],[pol4_pol[0], pol4_pol[8]])
plot_bandpass(mfi417_pol,outdir+'mfi_417_pol.png')
mfi419_pol = combine_bandpasses([pol4_pol[0], pol4_pol[1]], [pol4_pol[0], pol4_pol[3]],[pol4_pol[0], pol4_pol[5]],[pol4_pol[0], pol4_pol[7]])
plot_bandpass(mfi419_pol,outdir+'mfi_419_pol.png')
plot_bandpass_all([mfi111, mfi113, mfi217, mfi219, mfi311, mfi313, mfi417, mfi419], outdir+'mfi_comb_band.png')
# And the weighted combination
W417 = np.asarray([0.3616885 , 0.73196787, 0.73196787])
W419 = np.asarray([0.41912782, 0.7880423 , 0.7880423 ])
W217 = 1.0-W417
W219 = 1.0-W419
mfi17 = np.asarray(mfi217)*W217[0] + np.asarray(mfi417)*W417[0]
mfi19 = np.asarray(mfi219)*W219[0] + np.asarray(mfi419)*W419[0]
mfi17_pol = np.asarray(mfi217_pol)*W217[1] + np.asarray(mfi417_pol)*W417[1]
mfi19_pol = np.asarray(mfi219_pol)*W219[1] + np.asarray(mfi419_pol)*W419[1]
plot_bandpass(mfi17,outdir+'mfi_17.png')
plot_bandpass(mfi19,outdir+'mfi_19.png')
plot_bandpass(mfi17_pol,outdir+'mfi_17_pol.png')
plot_bandpass(mfi19_pol,outdir+'mfi_19_pol.png')
mfi111_corrections = np.ones(len(alphas))
mfi113_corrections = np.ones(len(alphas))
mfi217_corrections = np.ones(len(alphas))
mfi219_corrections = np.ones(len(alphas))
mfi217_pol_corrections = np.ones(len(alphas))
mfi219_pol_corrections = np.ones(len(alphas))
mfi311_corrections = np.ones(len(alphas))
mfi311pol_corrections = np.ones(len(alphas))
# mfi311orig_corrections = np.ones(len(alphas))
mfi313_corrections = np.ones(len(alphas))
mfi313pol_corrections = np.ones(len(alphas))
# mfi313orig_corrections = np.ones(len(alphas))
mfi417_corrections = np.ones(len(alphas))
mfi419_corrections = np.ones(len(alphas))
mfi417_pol_corrections = np.ones(len(alphas))
mfi419_pol_corrections = np.ones(len(alphas))
mfi17_corrections = np.ones(len(alphas))
mfi19_corrections = np.ones(len(alphas))
mfi17_pol_corrections = np.ones(len(alphas))
mfi19_pol_corrections = np.ones(len(alphas))
# Trim the bandpasses
# mfi111 = trim_bandpass(mfi111, 11.2, 4.5)
# mfi113 = trim_bandpass(mfi111, 12.8, 4.5)
mfi217 = trim_bandpass(mfi217, 16.7, 4.5)
mfi217_pol = trim_bandpass(mfi217_pol, 16.7, 4.5)
mfi219 = trim_bandpass(mfi219, 18.7, 4.5)
mfi219_pol = trim_bandpass(mfi219_pol, 18.7, 4.5)
mfi311 = trim_bandpass(mfi311, 11.1, 4.5)
mfi311pol = trim_bandpass(mfi311pol, 11.1, 4.5)
mfi313 = trim_bandpass(mfi313, 12.9, 4.5)
mfi313pol = trim_bandpass(mfi313pol, 12.9, 4.5)
mfi417 = trim_bandpass(mfi417, 17.0, 4.5)
mfi417_pol = trim_bandpass(mfi417_pol, 17.0, 4.5)
mfi419 = trim_bandpass(mfi419, 19.0, 4.5)
mfi419_pol = trim_bandpass(mfi419_pol, 19.0, 4.5)
mfi17 = trim_bandpass(mfi17, 16.8, 4.5)
mfi17_pol = trim_bandpass(mfi17_pol, 16.8, 4.5)
mfi19 = trim_bandpass(mfi19, 18.8, 4.5)
mfi19_pol = trim_bandpass(mfi19_pol, 18.8, 4.5)
reference_freqs = [11.2, 12.8, 16.7, 18.7, 11.1, 12.9, 17.0, 19.0, 16.8, 18.8]
reference_freqs = [11.2, 12.8, 16.807675, 18.831303, 10.987673, 12.901321, 16.938613, 18.901596, 16.8, 18.8]
calindex = 2.0#-0.3
# calindex = -0.31
usecorr = False
for i in range(0,len(alphas)):
mfi111_corrections[i] = calc_correction(mfi111, reference_freqs[0], alphas[i],calindex=calindex,usecorr=usecorr)
mfi113_corrections[i] = calc_correction(mfi113, reference_freqs[1], alphas[i],calindex=calindex,usecorr=usecorr)
mfi217_corrections[i] = calc_correction(mfi217, reference_freqs[2], alphas[i],calindex=calindex,usecorr=usecorr)
mfi219_corrections[i] = calc_correction(mfi219, reference_freqs[3], alphas[i],calindex=calindex,usecorr=usecorr)
mfi217_pol_corrections[i] = calc_correction(mfi217_pol, reference_freqs[2], alphas[i],calindex=calindex,usecorr=usecorr)
mfi219_pol_corrections[i] = calc_correction(mfi219_pol, reference_freqs[3], alphas[i],calindex=calindex,usecorr=usecorr)
mfi311_corrections[i] = calc_correction(mfi311, reference_freqs[4], alphas[i],calindex=calindex,usecorr=usecorr)
mfi313_corrections[i] = calc_correction(mfi313, reference_freqs[5], alphas[i],calindex=calindex,usecorr=usecorr)
mfi311pol_corrections[i] = calc_correction(mfi311pol, reference_freqs[4], alphas[i],calindex=calindex,usecorr=usecorr)
mfi313pol_corrections[i] = calc_correction(mfi313pol, reference_freqs[5], alphas[i],calindex=calindex,usecorr=usecorr)
# mfi311orig_corrections[i] = calc_correction(mfi311orig, 11.1, alphas[i],calindex=calindex,usecorr=usecorr)
# mfi313orig_corrections[i] = calc_correction(mfi313orig, 12.9, alphas[i],calindex=calindex,usecorr=usecorr)
mfi417_corrections[i] = calc_correction(mfi417, reference_freqs[6], alphas[i],calindex=calindex,usecorr=usecorr)
mfi419_corrections[i] = calc_correction(mfi419, reference_freqs[7], alphas[i],calindex=calindex,usecorr=usecorr)
mfi417_pol_corrections[i] = calc_correction(mfi417_pol, reference_freqs[6], alphas[i],calindex=calindex,usecorr=usecorr)
mfi419_pol_corrections[i] = calc_correction(mfi419_pol, reference_freqs[7], alphas[i],calindex=calindex,usecorr=usecorr)
mfi17_corrections[i] = calc_correction(mfi17, reference_freqs[8], alphas[i],calindex=calindex,usecorr=usecorr)
mfi19_corrections[i] = calc_correction(mfi19, reference_freqs[9], alphas[i],calindex=calindex,usecorr=usecorr)
mfi17_pol_corrections[i] = calc_correction(mfi17_pol, reference_freqs[8], alphas[i],calindex=calindex,usecorr=usecorr)
mfi19_pol_corrections[i] = calc_correction(mfi19_pol, reference_freqs[9], alphas[i],calindex=calindex,usecorr=usecorr)
print(alphas)
print('From bandpasses:')
print(mfi111_corrections)
print(mfi113_corrections)
print(mfi217_corrections)
print(mfi219_corrections)
print(mfi217_pol_corrections)
print(mfi219_pol_corrections)
print(mfi311_corrections)
print(mfi313_corrections)
print(mfi311pol_corrections)
print(mfi313pol_corrections)
# print(mfi311orig_corrections)
# print(mfi313orig_corrections)
print(mfi417_corrections)
print(mfi419_corrections)
print(mfi417_pol_corrections)
print(mfi419_pol_corrections)
print(mfi17_corrections)
print(mfi19_corrections)
print(mfi17_pol_corrections)
print(mfi19_pol_corrections)
exit()
# print('From fastcc:')
# print(fastcc('Q111',alpha=alphas))
# print(fastcc('Q113',alpha=alphas))
# print(fastcc('Q217',alpha=alphas))
# print(fastcc('Q219',alpha=alphas))
# print(fastcc('Q311',alpha=alphas))
# print(fastcc('Q313',alpha=alphas))
plt.plot(alphas,mfi111_corrections,label='111')
plt.plot(alphas,mfi113_corrections,label='113')
plt.plot(alphas,mfi217_corrections,label='217')
plt.plot(alphas,mfi219_corrections,label='219')
plt.plot(alphas,mfi311_corrections,label='311')
plt.plot(alphas,mfi313_corrections,label='313')
# plt.plot(alphas,mfi311orig_corrections,'-g',label='311orig')
# plt.plot(alphas,mfi313orig_corrections,'-.g',label='313orig')
plt.plot(alphas,mfi417_corrections,label='417')
plt.plot(alphas,mfi419_corrections,label='419')
l = plt.legend(prop={'size':11})
l.set_zorder(20)
plt.xlabel('Spectral index')
plt.ylabel('Colour correction')
plt.tight_layout()
plt.savefig(outdir+'mfi_corrections.png')
plt.clf()
plt.close()
plt.plot(alphas,mfi311_corrections,'r-',label='11')
plt.plot(alphas,mfi313_corrections,'g-',label='13')
plt.plot(alphas,mfi17_corrections,'b-',label='17')
plt.plot(alphas,mfi19_corrections,'m-',label='19')
plt.plot(alphas,mfi311pol_corrections,'r-.',label='11p')
plt.plot(alphas,mfi313pol_corrections,'g-.',label='13p')
plt.plot(alphas,mfi17_pol_corrections,'b-.',label='17p')
plt.plot(alphas,mfi19_pol_corrections,'m-.',label='19p')
l = plt.legend(prop={'size':11})
l.set_zorder(20)
plt.xlabel('Spectral index')
plt.ylabel('Colour correction')
plt.tight_layout()
plt.savefig(outdir+'mfi_band_corrections.png')
plt.clf()
plt.close()
plt.plot(alphas,mfi111_corrections,label='111')
plt.plot(alphas,mfi113_corrections,label='113')
plt.plot(alphas,mfi217_pol_corrections,label='217pol')
plt.plot(alphas,mfi219_pol_corrections,label='219pol')
plt.plot(alphas,mfi311pol_corrections,label='311pol')
plt.plot(alphas,mfi313pol_corrections,label='313pol')
plt.plot(alphas,mfi417_pol_corrections,label='417pol')
plt.plot(alphas,mfi419_pol_corrections,label='419pol')
l = plt.legend(prop={'size':11})
l.set_zorder(20)
plt.xlabel('Spectral index')
plt.ylabel('Colour correction')
plt.tight_layout()
plt.savefig(outdir+'mfi_corrections_pol.png')
plt.clf()
plt.close()
plt.plot(alphas,mfi217_corrections,'b--',label='217')
plt.plot(alphas,mfi219_corrections,'g--',label='219')
plt.plot(alphas,mfi417_corrections,'b-.',label='417')
plt.plot(alphas,mfi419_corrections,'g-.',label='419')
plt.plot(alphas,mfi217_pol_corrections,'r--',label='217pol')
plt.plot(alphas,mfi219_pol_corrections,'m--',label='219pol')
plt.plot(alphas,mfi417_pol_corrections,'r-.',label='417pol')
plt.plot(alphas,mfi419_pol_corrections,'m-.',label='419pol')
plt.plot(alphas,mfi17_corrections,'b-',label='17')
plt.plot(alphas,mfi19_corrections,'g-',label='19')
plt.plot(alphas,mfi17_pol_corrections,'r-',label='17pol')
plt.plot(alphas,mfi19_pol_corrections,'m-',label='19pol')
l = plt.legend(prop={'size':11})
l.set_zorder(20)
plt.xlabel('Spectral index')
plt.ylabel('Colour correction')
plt.savefig(outdir+'mfi_corrections_17_19.pdf')
plt.clf()
plt.close()
print('Horn 1:')
params = np.polyfit(alphas,mfi111_corrections,2)
print("'Q111': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 11.2],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='111_fit')
params = np.polyfit(alphas,mfi113_corrections,2)
print("'Q113': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 12.8],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='113_fit')
print('Horn 2:')
params = np.polyfit(alphas,mfi217_corrections,2)
print("'Q217': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 16.7],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='217_fit')
params = np.polyfit(alphas,mfi219_corrections,2)
print("'Q219': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 18.7],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='219_fit')
params = np.polyfit(alphas,mfi217_pol_corrections,2)
print("'Q217p': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 16.7],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='217pol_fit')
params = np.polyfit(alphas,mfi219_pol_corrections,2)
print("'Q219p': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 18.7],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='219pol_fit')
print('Horn 3:')
params = np.polyfit(alphas,mfi311_corrections,2)
print("'Q311': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 11.1],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='311_fit')
params = np.polyfit(alphas,mfi313_corrections,2)
print("'Q313': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 12.9],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='313_fit')
params = np.polyfit(alphas,mfi311pol_corrections,2)
print("'Q311p': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 11.1],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='311pol_fit')
params = np.polyfit(alphas,mfi313pol_corrections,2)
print("'Q313p': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 12.9],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='313pol_fit')
print('Horn 4:')
params = np.polyfit(alphas,mfi417_corrections,2)
print("'Q417': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 17.0],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='417_fit')
params = np.polyfit(alphas,mfi419_corrections,2)
print("'Q419': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 19.0],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='419_fit')
params = np.polyfit(alphas,mfi417_pol_corrections,2)
print("'Q417p': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 17.0],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='417pol_fit')
params = np.polyfit(alphas,mfi419_pol_corrections,2)
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='419pol_fit')
print("'Q419p': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 19.0],')
l = plt.legend(prop={'size':11})
params = np.polyfit(alphas,mfi17_corrections,2)
print("'Q17': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 16.8],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='17_fit')
params = np.polyfit(alphas,mfi19_corrections,2)
print("'Q19': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 18.8],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='19_fit')
params = np.polyfit(alphas,mfi17_pol_corrections,2)
print("'Q17p': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 16.8],')
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='17pol_fit')
params = np.polyfit(alphas,mfi19_pol_corrections,2)
plt.plot(alphas,params[2]+params[1]*alphas+params[0]*alphas**2,'-',label='19pol_fit')
print("'Q19p': ["+str(params[2]) + ', ' + str(params[1]) + ', ' + str(params[0]) + ', 18.8],')
l = plt.legend(prop={'size':11})
l.set_zorder(20)
plt.xlabel('Spectral index')
plt.ylabel('Colour correction')
plt.savefig(outdir+'mfi_corrections_fit.png')
plt.clf()
plt.close()
# Compare with Roke's values
roke_11 = [0.9589, 0.01601, 0.002185]
roke_13 = [1.022, -0.01436, 0.001688]
roke_17 = [1.029, -0.01618, 0.0008349]
roke_19 = [1.029, -0.01588, 0.0007932]
plt.plot(alphas,1.0/(roke_11[0]+roke_11[1]*alphas+roke_11[2]*alphas**2),'-',label='11_roke')
plt.plot(alphas,1.0/(roke_13[0]+roke_13[1]*alphas+roke_13[2]*alphas**2),'-',label='13_roke')
plt.plot(alphas,1.0/(roke_17[0]+roke_17[1]*alphas+roke_17[2]*alphas**2),'-',label='17_roke')
plt.plot(alphas,1.0/(roke_19[0]+roke_19[1]*alphas+roke_19[2]*alphas**2),'-',label='19_roke')
l = plt.legend(prop={'size':11})
l.set_zorder(20)
plt.savefig(outdir+'mfi_corrections_roke.pdf')
plt.clf()
plt.close()
print(mfi219_corrections[alphas == -0.5] - mfi419_corrections[alphas == -0.5])
print(mfi217_corrections[alphas == -0.5] - mfi417_corrections[alphas == -0.5])
exit()