-
Notifications
You must be signed in to change notification settings - Fork 0
/
data.py
342 lines (299 loc) · 23.5 KB
/
data.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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 27 11:27:27 2020
@author: félix langot
"""
# This is a program to create the for MSci Project
import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from scipy.integrate import quad
import seaborn as sns
sns.set(color_codes=True)
# Compute gamma functions
def Xi(beta):
return (gamma(1.5*beta)/(gamma((1.5*beta)-0.5)))**2 * (gamma((3*beta)-0.5)/gamma(3*beta))
# Compute the cosmological integration factor which depends on z.
def hiya(x):
return ( ((1+x)**2 * (1 + (0.3*x))) - (x*(2+x)*0.7) )**(-0.5)
def FBeta(x,t):
return ( (1 + ((x**2)/(t**2)) )**((1-3*beta[i])/2) )*x
# Equation to determine errors from
def D_a_eq(DT, SX, T_e, z, beta, r_c, L_ee,ff):
return (DT)**2 * (1 / (SX * (10000/3600))) * (9.1093856E-31 * 299792458**2 / (T_e * 1000 * 1.60217662E-19))**2 * ((L_ee/1000000)/(4*np.pi*(ff**2)*(2.72548*6.6524587158E-29)**2 * (1+z)**4)) * 1/(2*np.sqrt(np.pi)) * (Xi(beta)/r_c) / 3.09E+25
def D_a_eq_old(DT, SX, T_e, z, beta, r_c, L_ee):
return (DT)**2 * (1 / (SX * (10000/3600))) * (9.1093856E-31 * 299792458**2 / (T_e * 1000 * 1.60217662E-19))**2 * ((L_ee/1000000)/(4*np.pi*(4)*(2.72548*6.6524587158E-29)**2 * (1+z)**4)) * 1/(2*np.sqrt(np.pi)) * (Xi(beta)/r_c) / 3.09E+25
NumData = 22
TypeData = 1 # 0 # for used data
obsid = np.array([13395,13468,13465,13462,13403,13396,13486,13398,13476,13474,13488,14448,13397,13483,13491,15573,13490,13485,13478,13467,13473])
Exposure_time = np.array([7.16,14.88,10.4,12.23,8.96,6.77,16.45,12.29,20.32,33.44,41.25,30.78,22.63,19.22,17.36,60.28,16.46,29.2,49.35,70.51,24.4,43.95])
RA = np.array([10.204836,16.615226,57.073736,3.3289942,38.678636,62.81544,103.9626,326.46821,30.143591,46.951617,74.116261,327.1812,44.099676,64.34505,8.476725,20.792265,40.861544,58.236718,15.669013,310.82839,25.545159,72.274168])
Dec = np.array([-44.132926,-59.721357,-45.250966,-49.115117,-58.52141,-48.321752,-52.567743,-56.747557,-48.875691,-50.707144,-51.276841,-61.277969,-56.297961,-47.813921,-63.44628,-48.358848,-59.512436,-56.799649,-46.064726,-50.593839,-50.543772,-49.024604])
nH = np.array([0.0263,0.0287,0.0144,0.0226,0.035,0.0165,0.0519,0.0302,0.0235,0.021,0.0186,0.0295,0.0236,0.0182,0.0164,0.0196,0.0357,0.0258,0.0154,0.0296,0.0288,0.0172])
theta_c = np.array([0.50,0.50,0.75,0.75,0.50,1.00,0.50,0.50,0.75,0.50,1.00,0.75,0.50,0.25,0.75,0.25,0.25,0.75,0.25,0.50,0.75,0.50])
XiC = np.array([0.355744155675953,0.325154970486976,0.258658205718938,0.529352757978805,0.26471394536504,0.101543245948168,0.467596538538188,0.301907581418201,0.305097801458771,0.309351236365546,0.260402165571247,0.359404412249901,0.387954025254074,0.242212835449274,0.163895486172248,0.234798765171727,0.400021352893062,0.362323189272605,0.0941227833561085,0.207133426291054,0.367466358197207,0.534404747587648])
# Last Data
z = np.array([0.35,0.348114,0.3582,0.4064,0.415,0.423513,0.4703,0.48,0.498492,0.55,0.561539,0.571,0.58,0.58105,0.597129,0.62,0.635236,0.67,0.721856,0.7234,0.73,0.792])
S_sherpa = np.array([3.5163E-06,2.2282E-06,1.0636E-06,1.5656E-06,1.7906E-05,2.1272E-06,5.0978E-07,1.4727E-06,8.3170E-07,9.8725E-07,7.9385E-07,5.2856E-07,5.5472E-07,1.26E-05,1.19E-06,4.7044E-07,1.4990E-06,5.74E-07,3.91E-07,1.58E-05,4.59E-07,2.75E-07])
S_sherpa_err = np.array([3.1583E-07,2.0103E-07,1.2896E-07,1.3528E-07,1.5549E-06,8.8426E-07,6.2137E-08,1.5683E-07,9.5071E-08,1.0313E-07,6.2915E-08,4.1539E-08,7.0598E-08,9.04E-07,2.10E-07,4.3797E-08,1.9979E-07,8.25E-08,7.89E-08,9.44E-07,5.72E-08,3.11E-08])
r_sherpa = np.array([69.4547,58.1855,69.3704,114.258,21.1456,22.9287,98.8143,62.3842,51.3202,56.0963,52.9654,83.5413,86.1231,19.111,29.4423,52.2606,54.0617,64.9069,28.1084,11.0021,79.7676,117.488])
r_err = np.array([3.84321372,6.07833,5.14755,8.5691148,1.3613886,2.48688288,22.2053376,4.24377552,4.30699752,4.90001004,2.72045496,3.59630844,9.1127256,0.61617588,2.87800812,3.2774826,4.050759,6.717276,4.24356396,0.25610568,8.0932032,14.4964356])
Y = np.array([0.00021948983,0.000115115545,0.00010564338,0.00013506584,0.00016898503,0.00015965472,0.00009097335,0.00014178356,0.000083517225,0.00009526548,0.00009000508,0.00007919483,0.00008716427,0.00017264496,0.00008216682,0.0000911049,0.00009497875,0.000078916,0.00010027521,0.00008734432,0.00011733797,0.00009631345])
Y_unc = np.array([0.000013548807,0.000015860673,0.000012103874,0.0000129467835,0.000014385522,0.000012672225,0.000015437567,0.000017436387,0.0000127068215,0.000013290628,0.000011417068,0.00001196352,0.0000143953,0.000026340524,0.000013893156,0.000014047431,0.000020687343,0.0000122120655,0.000023835268,0.000016675627,0.000012623081,0.000012045249])
beta = np.array([0.757964,0.718976,0.639373,1.01192,0.64635,0.472049,0.914843,0.690375,0.694249,0.699439,0.641377,0.762736,0.800767,0.62068,0.53608,0.61237,0.817283,0.766558,0.464444,0.581962,0.773329,1.02021])
beta_unc = np.array([0.0497747,0.0437258,0.0450767,0.133337,0.0142615,0.0176717,0.449146,0.0530611,0.0548642,0.0779211,0.0356317,0.0467357,0.121666,0.0140217,0.0323047,0.0385831,0.0782782,0.101498,0.0448218,0.0079473,0.111659,0.242764])
L_ee = np.array([2.916970130922960E-15,2.999754261529680E-15,2.778408137182330E-15,2.664333680951090E-15,2.855586409086370E-15,3.082946757535740E-15,2.631978942371750E-15,2.756855924168490E-15,2.792535701460220E-15,2.670882260776740E-15,2.795848614997560E-15,2.716237471147260E-15,2.601849503212530E-15,2.800291096316100E-15,2.870238370954160E-15,2.803084796316230E-15,2.617963995147330E-15,2.468156927059060E-15,2.679436354029690E-15,2.679543257772630E-15,2.489211773767810E-15,2.425355120275520E-15])
L_ee_unc = np.array([1.221379003322260E-16,1.244638432178370E-16,1.845808421541120E-16,8.433782711085340E-17,1.234505730455250E-16,1.613048193286620E-16,1.795253591145390E-16,1.393071969719100E-16,2.174999380511920E-16,1.039761071640040E-16,1.252462964241650E-16,1.380184170110900E-16,1.928251639896100E-16,1.040340549926890E-16,2.825358856325630E-16,1.277011921858670E-16,1.692037208776780E-16,1.124482514946140E-16,1.436202154606410E-16,7.103956195555910E-17,2.298564952807410E-16,9.512855803898500E-17])
Te = np.array([7.78014,6.06551,5.02654,7.6613,5.43676,5.93566,7.90022,8.96459,10.1339,8.07639,7.66745,6.9855,7.99919,5.83822,5.10165666666667,6.57803,5.43698,6.01785,6.10179,5.08212,7.02854,10.1561])
Te_unc = np.array([0.951414,0.579698,1.05384,0.965753,0.708637,0.601316,1.57905,1.70625,2.89344,0.846069,0.92777,1.40405,1.18775,0.487224,0.758313666666667,0.672422,0.72233,0.747407,0.906228,0.247599,0.975405,1.65734])
#theta_c[15] = 0.5
#theta_c[18] = 0.5
# Constants
h = 6.62607004E-34
nu = np.array([150E9, 95E9, 220E9])
k_B = 1.38E-23
T_CMB = 2.72548
cor_s = 0.492
# Initilizations
xx = np.zeros(len(nu))
ff = np.zeros(len(nu))
FacY = np.zeros((NumData,len(nu)))
y0 = np.zeros((NumData,len(nu)))
y0_unc = np.zeros((NumData,len(nu)))
DT = np.zeros((NumData,len(nu)))
DT_unc = np.zeros((NumData,len(nu)))
D_a = np.zeros((NumData,len(nu)))
r_c = np.zeros(NumData)
r_c_unc = np.zeros(NumData)
SX = np.zeros(NumData)
SX_unc = np.zeros(NumData)
for i in np.arange(len(nu)):
xx[i] = h*nu[i]/(k_B*T_CMB)
ff[i] = xx[i]*((np.exp(xx[i])+1)/(np.exp(xx[i])-1)) - 4
for j in np.arange(len(nu)):
for i in np.arange(NumData):
a = theta_c[i]
RES,ERR = quad(FBeta, 0, 0.75, args=(a))
FacY[i,j] = 2*np.pi*RES
y0[i,j] = Y[i]/FacY[i,j]
y0_unc[i,j] = Y_unc[i]/FacY[i,j]
DT[i,j] = T_CMB*(ff[j]*y0[i,j])
DT_unc[i,j] = T_CMB*(ff[j]*y0_unc[i,j])
r_c = cor_s*r_sherpa
r_c_unc = cor_s*r_err
SX = (1/Exposure_time)*3600*100*S_sherpa*1000
SX_unc = (1/Exposure_time)*3600*100*S_sherpa_err*1000
# Initial Data
# Redshift and its unc
z1 = np.array([0.35,0.348114,0.3582,0.4064,0.415,0.423513,0.4703,0.48,0.498492,0.55,0.561539,0.571,0.58,0.58105,0.597129,0.62,0.635236,0.67,0.721856,0.7234,0.73,0.792])
z1_unc = z / 15
# Delta T and its unc
DT1 = np.array([-0.00155430052805673,-0.000815181971670572,-0.000748105556051039,-0.000956458467503602,-0.00119665463025181,-0.00113058275001967,-0.000644220854989454,-0.0010040294904678,-0.000591420873210084,-0.000674615246955078,-0.000637364124669414,-0.000560812161950115,-0.000617247144838919,-0.0012225721460276,-0.000581858082967865,-0.000645152416303552,-0.000672584790280117,-0.000558837648523966,-0.000710091268711629,-0.000618522155212186,-0.000830919905182418,-0.00068203636675998])
DT1_unc = np.array([9.59448457117068E-05,0.000112316148858629,8.57126626310302E-05,9.16816621102044E-05,0.00010186997931053,8.97373969862464E-05,0.000109319955917826,0.000123474447638424,8.99823894682166E-05,9.41165707705255E-05,0.000080849098207692,8.47187564609134E-05,0.000101939221473428,0.00018652841620266,9.83833269503857E-05,9.94758136226199E-05,0.000146495845155973,8.64788125049503E-05,0.000168787636487639,0.000118087183591956,8.93893874900851E-05,8.52975141548692E-05])
# Beta and its unc
beta1 = np.array([0.757964,0.718976,0.813157,0.810782,0.608069,0.557648,1.24219,0.690375,0.644117,0.781599,0.641377,0.943045,0.800767,0.590183,0.529259,0.61237,0.670853,0.766558,0.710763,0.55382,0.773329,1.02021])
beta1_unc = np.array([0.0497747,0.0437258,0.103342,0.133337,0.0142615,0.0176717,0.449146,0.0530611,0.0475828,0.0652756,0.0356317,0.148997,0.121666,0.0121917,0.0323047,0.0385831,0.0460435,0.101498,0.101507,0.00653156,0.111659,0.242764])
# S_X and its unc
SX1 = np.array([176.795195530726,53.9083064516129,33.1936961538461,53.1816843826656,761.323660714286,113.112850812408,9.64443282674772,32.2290820668693,15.0632716535433,9.56680622009569,6.92811927272727,4.25671854545454,8.8244701723376,248.107804370447,24.0627649769585,2.80953483742535,35.4755285540705,7.07238493150685,2.04445714285714,86.272869096582,6.76774918032787,2.25395221843003])
SX1_unc = np.array([15.8798212290503,4.86360483870968,4.02220384615385,3.98195911692559,62.4720535714286,47.0213707533235,1.35983927051672,4.59382262001627,1.73326535433071,0.829114234449761,0.549078109090909,0.554991812865497,1.12307591692444,19.1118418314256,4.3619930875576,0.261560783012608,4.36973511543135,1.01709739726027,0.329428085106383,4.82065579350447,0.843260163934426,0.254435221843003])
# r_c and its unc
r_c1 = np.array([34.1717124,28.627266,46.0875096,41.4140508,9.2523552,22.9045188,66.10266,30.6930264,27.3035892,34.0870884,26.0589768,51.560124,42.3725652,8.4457704,14.2696728,25.7122152,20.8047612,31.9341948,31.289724,4.7809854,39.2456592,57.804096])
r_c1_unc = np.array([3.84321372,6.07833,8.4177264,8.5691148,0.68069676,2.48688288,22.2053376,4.24377552,3.99028236,4.31270964,2.72045496,9.399168,9.1127256,0.59322408,2.87800812,3.2774826,2.98504764,6.717276,7.5943644,0.241941492,8.0932032,14.4964356])
# T_e and its unc
Te1 = np.array([7.78014,6.06551,5.02654,7.6613,6.23959,5.93566,7.90022,8.96459,10.1339,8.07639,7.66745,6.9855,7.99919,6.17429,6.41017,6.57803,7.88763,6.01785,6.10179,5.08212,7.02854,10.1561])
Te1_unc = np.array([0.951414,0.579698,1.05384,0.965753,0.606944,0.601316,1.57905,1.70625,2.89344,0.846069,0.92777,1.40405,1.18775,0.52882,0.960366,0.672422,0.998255,0.747407,0.906228,0.247599,0.975405,1.65734])
# Lambda and its unc
L_ee1 = np.array([2.91697013092296E-15,2.99975426152968E-15,2.77840813718233E-15,2.66433368095109E-15,2.82664776993781E-15,3.08294675753574E-15,2.63197894237175E-15,2.75685592416849E-15,2.79253570146022E-15,2.67088226077674E-15,2.79584861499756E-15,2.71623747114726E-15,2.60184950321253E-15,2.73536837848511E-15,2.75759905538023E-15,2.80308479631623E-15,2.57726979991298E-15,2.46815692705906E-15,2.67943635402969E-15,2.67954325777263E-15,2.48921177376781E-15,2.42535512027552E-15])
L_ee1_unc = np.array([1.22137900332226E-16,1.24463843217837E-16,1.84580842154112E-16,8.43378271108534E-17,1.11112885330488E-16,1.61304819328662E-16,1.79525359114539E-16,1.3930719697191E-16,2.17499938051192E-16,1.03976107164004E-16,1.25246296424165E-16,1.3801841701109E-16,1.9282516398961E-16,9.89977727698239E-17,2.14357029516505E-16,1.27701192185867E-16,1.32102358348748E-16,1.12448251494614E-16,1.43620215460641E-16,7.10395619555591E-17,2.29856495280741E-16,9.5128558038985E-17])
# values computed April 10th
# Redshift and its unc
z0 = np.array([0.35,0.348114,0.3582,0.4064,0.415,0.423513,0.4703,0.48,0.498492,0.55,0.561539,0.571,0.58,0.58105,0.597129,0.62,0.635236,0.67,0.721856,0.7234,0.73,0.792])
z0_unc = z / 15
# Delta T and its unc
DT0 = np.array([-0.00155430052805673,-0.000815181971670572,-0.000748105556051039,-0.000956458467503602,-0.00119665463025181,-0.00113058275001967,-0.000644220854989454,-0.0010040294904678,-0.000591420873210084,-0.000674615246955078,-0.000637364124669414,-0.000560812161950115,-0.000617247144838919,-0.0012225721460276,-0.000581858082967865,-0.000645152416303552,-0.000672584790280117,-0.000558837648523966,-0.000710091268711629,-0.000618522155212186,-0.000830919905182418,-0.00068203636675998])
DT0_unc = np.array([9.59448457117068E-05,0.000112316148858629,8.57126626310302E-05,9.16816621102044E-05,0.00010186997931053,8.97373969862464E-05,0.000109319955917826,0.000123474447638424,8.99823894682166E-05,9.41165707705255E-05,0.000080849098207692,8.47187564609134E-05,0.000101939221473428,0.00018652841620266,9.83833269503857E-05,9.94758136226199E-05,0.000146495845155973,8.64788125049503E-05,0.000168787636487639,0.000118087183591956,8.93893874900851E-05,8.52975141548692E-05])
# Beta and its unc
beta0 = np.array([0.757964,0.718976,0.639373,1.01192,0.64635,0.472049,1.24219,0.690375,0.694249,0.699439,0.641377,0.664896,0.800767,0.62068,0.53608,0.61237,0.817283,0.766558,0.464444,0.581962,0.773329,1.02021])
beta0_unc = np.array([0.0497747,0.0437258,0.0450767,0.133337,0.0142615,0.0176717,0.449146,0.0530611,0.0548642,0.0779211,0.0356317,0.0467357,0.121666,0.0140217,0.0323047,0.0385831,0.0782782,0.101498,0.0448218,0.0079473,0.111659,0.242764])
# S_X and its unc
SX0 = np.array([176.7951955,53.90830645,36.81623077,46.08412101,719.4334821,113.1128508,9.644432827,32.22908207,14.73475394,10.62824641,6.928119273,4.001986909,8.824470172,235.3130073,24.68073733,2.809534837,32.78580802,7.072384932,2.854431611,80.64592256,6.76774918,2.253952218])
SX0_unc = np.array([15.87982123,4.863604839,4.463930769,3.981959117,62.47205357,47.02137075,1.359839271,4.59382262,1.684325197,1.110251196,0.549078109,0.485839766,1.123075917,16.92756712,4.361993088,0.261560783,4.369735115,1.017097397,0.575343465,4.820655794,0.843260164,0.254435222])
# r_c and its unc
r_c0 = np.array([34.1717124,28.627266,34.1302368,56.214936,10.4036352,11.2809204,66.10266,30.6930264,30.0408804,27.5993796,26.0589768,28.5721128,42.3725652,9.402612,14.4856116,25.7122152,26.5983564,31.9341948,13.8293328,5.4130332,39.2456592,57.804096])
r_c0_unc = np.array([3.84321372,6.07833,5.14755,8.5691148,1.3613886,2.48688288,22.2053376,4.24377552,4.30699752,4.90001004,2.72045496,3.59630844,9.1127256,0.61617588,2.87800812,3.2774826,4.050759,6.717276,4.24356396,0.25610568,8.0932032,14.4964356])
# T_e and its unc
Te0 = np.array([7.78014,6.06551,5.02654,7.6613,5.43676,5.93566,7.90022,8.96459,10.1339,8.07639,7.66745,6.9855,7.99919,5.83822,5.101656667,6.57803,5.43698,6.01785,6.10179,5.08212,7.02854,10.1561])
Te0_unc = np.array([0.951414,0.579698,1.05384,0.965753,0.708637,0.601316,1.57905,1.70625,2.89344,0.846069,0.92777,1.40405,1.18775,0.487224,0.758313667,0.672422,0.72233,0.747407,0.906228,0.247599,0.975405,1.65734])
# Lambda and its unc
L_ee0 = np.array([2.91697013092296E-15,2.99975426152968E-15,2.77840813718233E-15,2.66433368095109E-15,2.85558640908637E-15,3.08294675753574E-15,2.63197894237175E-15,2.75685592416849E-15,2.79253570146022E-15,2.67088226077674E-15,2.79584861499756E-15,2.71623747114726E-15,2.60184950321253E-15,2.8002910963161E-15,2.87023837095416E-15,2.80308479631623E-15,2.61796399514733E-15,2.46815692705906E-15,2.67943635402969E-15,2.67954325777263E-15,2.48921177376781E-15,2.42535512027552E-15])
L_ee0_unc = np.array([1.22137900332226E-16,1.24463843217837E-16,1.84580842154112E-16,8.43378271108534E-17,1.23450573045525E-16,1.61304819328662E-16,1.79525359114539E-16,1.3930719697191E-16,2.17499938051192E-16,1.03976107164004E-16,1.25246296424165E-16,1.3801841701109E-16,1.9282516398961E-16,1.04034054992689E-16,2.82535885632563E-16,1.27701192185867E-16,1.69203720877678E-16,1.12448251494614E-16,1.43620215460641E-16,7.10395619555591E-17,2.29856495280741E-16,9.5128558038985E-17])
# DT
DT0 = np.array([-0.00155430052805673,-0.000815181971670572,-0.000748105556051039,-0.000956458467503602,-0.00119665463025181,-0.00113058275001967,-0.000644220854989454,-0.0010040294904678,-0.000591420873210084,-0.000674615246955078,-0.000637364124669414,-0.000560812161950115,-0.000617247144838919,-0.0012225721460276,-0.000581858082967865,-0.000645152416303552,-0.000672584790280117,-0.000558837648523966,-0.000710091268711629,-0.000618522155212186,-0.000830919905182418,-0.00068203636675998])
DT0_unc = np.array([9.59448457117068E-05,0.000112316148858629,8.57126626310302E-05,9.16816621102044E-05,0.00010186997931053,8.97373969862464E-05,0.000109319955917826,0.000123474447638424,8.99823894682166E-05,9.41165707705255E-05,0.000080849098207692,8.47187564609134E-05,0.000101939221473428,0.00018652841620266,9.83833269503857E-05,9.94758136226199E-05,0.000146495845155973,8.64788125049503E-05,0.000168787636487639,0.000118087183591956,8.93893874900851E-05,8.52975141548692E-05])
D_a1 = D_a_eq_old(DT1,SX1,Te1,z1,beta1,r_c1,L_ee1)
D_a0 = D_a_eq_old(DT0,SX0,Te0,z0,beta0,r_c0,L_ee0)
#for j in np.arange(len(nu)):
# DT[15,j] = DT0[15]*ff[j]/2
# DT[18,j] = DT0[18]*ff[j]/2
Vcolor = ['darkblue','red','green']
Vlabel = ['nu=150E9','nu=95E9','nu=220E9']
D_a_unc = np.zeros(NumData)
for j in np.arange(len(nu)):
D_a[:,j] = D_a_eq(DT[:,j],SX,Te,z,beta,r_c,L_ee,ff[j])
plt.errorbar(np.arange(NumData), D_a[:,j], yerr=D_a_unc,fmt='d',color=Vcolor[j],label=Vlabel[j])
plt.xlabel('#')
plt.ylabel('D_a')
plt.legend()
plt.show()
plt.errorbar(np.arange(NumData), D_a[:,0], yerr=D_a_unc,fmt='d',color=Vcolor[0],label=Vlabel[0])
plt.errorbar(np.arange(NumData), D_a0, yerr=D_a_unc,fmt='d',color=Vcolor[1],label='data 0')
#plt.errorbar(np.arange(NumData), D_a1, yerr=D_a_unc,fmt='d',color=Vcolor[2],label='data ini')
plt.xlabel('#')
plt.ylabel('D_a')
plt.legend(loc='upper left')
plt.show()
#====================================================
# Theory
#====================================================
c = 2.99792458 * 10**5 # in km/s
conv = 1000
Ndata = z.shape[0]
# Loop to calculate the integration factor for each z in the list
Ktup = []
Klist = []
errK = []
for i in np.arange(0, Ndata):
Ktup.append(integrate.quad(lambda x: hiya(x), 0, z[i]))
(value, error) = Ktup[i]
Klist.append(value)
errK.append(error)
K = np.asarray(Klist)
errD_a = np.array([0.35683888, 0.76760733, 1.18586014, 0.45285937, 0.22405864,
0.90606152, 2.21080192, 0.64893938, 0.66568483, 0.7993051 ,
0.95424936, 2.28865079, 0.79673148, 0.39894181, 0.83849681,
2.84199489, 0.77619527, 1.25201057, 4.92297751, 0.4311775 ,
1.47255965, 1.48010728])
N = 100
χ_v = np.zeros((N, 1))
a = 40*conv # in km/s/Gpc
b = 110*conv
H0_vec = np.arange(a, b, (b-a)/N)
Dfit_v = np.zeros((N,Ndata))
Dfit = np.zeros(Ndata)
# This loop calculates the χ^2 for different H0 values
for j in np.arange(0, N):
H0 = H0_vec[j]
Dfit[:] = ( 1/H0 )*(c*K)/(1+z)
χ_v[j] = np.sum( (Dfit[:]-D_a[:,0])**2/(errD_a[:]**2) )
minval = np.min(χ_v)
minpos = np.argmin(χ_v)
print('==============================================================================')
print('Estimated value for H =',H0_vec[minpos]/conv)
print('χ =',minval )
print('Reduced χ =',minval/(NumData-1) )
print('==============================================================================')
plt.errorbar(np.arange(NumData), r_c, yerr=r_c_unc,fmt='d',color='darkblue', label='last')
plt.errorbar(np.arange(NumData), r_c0, yerr=r_c_unc,fmt='d',color='red', label='bench')
#plt.errorbar(np.arange(NumData), r_c1, yerr=r_c_unc,fmt='d',color='darkgreen',label='init')
plt.xlabel('#')
plt.ylabel('r_c')
plt.legend()
plt.show()
plt.errorbar(np.arange(NumData), SX, yerr=SX_unc,fmt='d',color='darkblue', label='last')
plt.errorbar(np.arange(NumData), SX0, yerr=SX_unc,fmt='d',color='red', label='bench')
#plt.errorbar(np.arange(NumData), SX1, yerr=SX_unc,fmt='d',color='darkgreen',label='init')
plt.xlabel('#')
plt.ylabel('SX')
plt.legend()
plt.show()
plt.errorbar(np.arange(NumData), beta, yerr=beta_unc,fmt='d',color='darkblue', label='last')
plt.errorbar(np.arange(NumData), beta0, yerr=beta_unc,fmt='d',color='red', label='bench')
#plt.errorbar(np.arange(NumData), beta1, yerr=beta_unc,fmt='d',color='darkgreen',label='init')
plt.xlabel('#')
plt.ylabel('beta')
plt.legend()
plt.show()
plt.errorbar(np.arange(NumData), L_ee, yerr=L_ee_unc,fmt='d',color='darkblue', label='last')
plt.errorbar(np.arange(NumData), L_ee0, yerr=L_ee_unc,fmt='d',color='red', label='bench')
#plt.errorbar(np.arange(NumData), L_ee1, yerr=L_ee_unc,fmt='d',color='darkgreen',label='init')
plt.xlabel('#')
plt.ylabel('L_ee')
plt.legend()
plt.show()
plt.errorbar(np.arange(NumData), Te, yerr=Te_unc,fmt='d',color='darkblue', label='last')
plt.errorbar(np.arange(NumData), Te0, yerr=Te_unc,fmt='d',color='red', label='bench')
#plt.errorbar(np.arange(NumData), Te1, yerr=Te_unc,fmt='d',color='darkgreen',label='init')
plt.xlabel('#')
plt.ylabel('T_e')
plt.legend()
plt.show()
plt.errorbar(np.arange(NumData), DT[:,0], yerr=DT_unc[:,0],fmt='d',color='darkblue', label='last')
plt.errorbar(np.arange(NumData), DT0, yerr=DT_unc[:,0],fmt='d',color='red', label='bench')
#plt.errorbar(np.arange(NumData), DT1, yerr=DT_unc[:,0],fmt='d',color='darkgreen',label='init')
plt.xlabel('nu=150E9')
plt.ylabel('DT')
plt.legend()
plt.show()
plt.errorbar(np.arange(NumData), Y, yerr=Y_unc,fmt='d',color='darkblue')
plt.xlabel('#')
plt.ylabel('Y')
plt.show()
plt.plot(np.arange(NumData),1/Exposure_time,'o')
plt.xlabel('#')
plt.ylabel('1/Exposure_time')
plt.show()
plt.plot(np.arange(NumData),S_sherpa,'o')
plt.xlabel('#')
plt.ylabel('S_sherpa')
plt.show()
plt.errorbar(np.arange(NumData), y0[:,0], yerr=y0_unc[:,0],fmt='d',color='darkblue')
plt.xlabel('nu=150E9')
plt.ylabel('y0')
plt.show()
plt.errorbar(np.arange(NumData), y0[:,1], yerr=y0_unc[:,1],fmt='d',color='red')
plt.xlabel('nu=95E9')
plt.ylabel('y0')
plt.show()
plt.errorbar(np.arange(NumData), y0[:,2], yerr=y0_unc[:,2],fmt='d',color='darkgreen')
plt.xlabel('nu=220E9')
plt.ylabel('y0')
plt.show()
plt.errorbar(np.arange(NumData), DT[:,0], yerr=DT_unc[:,0],fmt='d',color='darkblue')
plt.xlabel('nu=150E9')
plt.ylabel('DT')
plt.show()
plt.errorbar(np.arange(NumData), DT[:,1], yerr=DT_unc[:,1],fmt='d',color='red')
plt.xlabel('nu=95E9')
plt.ylabel('DT')
plt.show()
plt.errorbar(np.arange(NumData), DT[:,2], yerr=DT_unc[:,2],fmt='d',color='darkgreen')
plt.xlabel('nu=220E9')
plt.ylabel('DT')
plt.show()
plt.scatter(XiC,Xi(beta))
plt.xlabel('XiC')
plt.ylabel('Xi(beta)')
plt.show()
fff= -2
beta520 = 0.686
Y520 = 0.73E-10
thetac520 = 80
def FFBeta(x):
return ( (1 + (x**2/thetac520**2) )**((1-3*beta520)/2) )*x
RES,ERR = quad(FFBeta, 0, 3.9*10**(-4))
print('RES =', RES)
FacY520 = RES
y0520 = Y520/FacY520
#y0_unc[i,j] = Y_unc[i]/FacY[i,j]
DT520 = T_CMB*(fff*y0520)
#DT_unc[i,j] = T_CMB*(ff[j]*y0_unc[i,j])
print(DT520)
truedt = -1.44e-3
fracdiff = truedt/DT520
print(fracdiff)
print(1/fracdiff)
DT0 = truedt
Y520 = truedt/(T_CMB*fff) * RES
print(Y520)