User: reisswig Date: 2013/04/28 11:56 AM

Modified: /trunk/DataAnalysis/ DataAnalysis.py

Log: Added Hubble constant. Added Luminosity distance caluclator. Added Big Bang Observer and DECIGO noise curves.

File Changes:

Directory: /trunk/DataAnalysis/ ===============================

File [modified]: DataAnalysis.py Delta lines: +214 -1 =================================================================== --- trunk/DataAnalysis/DataAnalysis.py 2012-11-14 14:27:38 UTC (rev 16) +++ trunk/DataAnalysis/DataAnalysis.py 2013-04-28 16:56:56 UTC (rev 17) @@ -19,6 +19,7 @@

constant_c = 299792458 # speed of light in [m/sec] constant_G = 6.67428 * 1e-11 # gravitational constant in [m^3 kg^-1 s^-2] +constant_H0 = 70.0 # Hubble constant

Msol = 1.9891*10**30 # solar mass[kg] ConversionFac = constant_G/constant_c**2 # [m/kg] (mass faktor G/c^2 to transform from geometrized units to SI units) @@ -379,7 +380,144 @@ return res

+def LuminosityDistance(z, Omega_m=0.3175, Omega_lambda=0.6825): + """Computes the luminosity distance in Mpc for a given redshift z + Code from Ed Wright http://www.astro.ucla.edu/~wright/CosmoCalc.html""" + + H0 = constant_H0 + WM = Omega_m # Omega(matter) + WV = Omega_lambda # Omega(vacuum) or lambda + + WR = 0. # Omega(radiation) + WK = 0. # Omega curvaturve = 1-Omega(total) + c = 299792.458 # velocity of light in km/sec + Tyr = 977.8 # coefficent for converting 1/H into Gyr + DTT = 0.5 # time from z to now in units of 1/H0 + DTT_Gyr = 0.0 # value of DTT in Gyr + age = 0.5 # age of Universe in units of 1/H0 + age_Gyr = 0.0 # value of age in Gyr + zage = 0.1 # age of Universe at redshift z in units of 1/H0 + zage_Gyr = 0.0 # value of zage in Gyr + DCMR = 0.0 # comoving radial distance in units of c/H0 + DCMR_Mpc = 0.0 + DCMR_Gyr = 0.0 + DA = 0.0 # angular size distance + DA_Mpc = 0.0 + DA_Gyr = 0.0 + kpc_DA = 0.0 + DL = 0.0 # luminosity distance + DL_Mpc = 0.0 + DL_Gyr = 0.0 # DL in units of billions of light years + V_Gpc = 0.0 + a = 1.0 # 1/(1+z), the scale factor of the Universe + az = 0.5 # 1/(1+z(object)) + + h = H0/100. + WR = 4.165E-5/(h*h) # includes 3 massless neutrino species, T0 = 2.72528 + WK = 1-WM-WR-WV + az = 1.0/(1+1.0*z) + age = 0. + n=1000 # number of points in integrals + for i in range(n): + a = az*(i+0.5)/n + adot = sqrt(WK+(WM/a)+(WR/(a*a))+(WV*a*a)) + age = age + 1./adot

+ zage = az*age/n + zage_Gyr = (Tyr/H0)*zage + DTT = 0.0 + DCMR = 0.0 + + def DCMT(): + ratio = 1.00 + x = sqrt(abs(WK))*DCMR + if (x > 0.1): + if (WK > 0): ratio = 0.5*(exp(x)-exp(-x))/x + else: ratio = sin(x)/x; + y = ratio*DCMR + return y + y = x*x + # statement below fixed 13-Aug-03 to correct sign error in expansion + if (WK < 0): y = -y + ratio = 1 + y/6 + y*y/120 + y= ratio*DCMR + return y + + + def VCM(): + ratio = 1.00 + x = sqrt(abs(WK))*DCMR; + if (x > 0.1): + if (WK > 0): ratio = (0.125*(exp(2*x)-exp(-2*x))-x/2)/(x*x*x/3) + else: ratio = (x/2 - sin(2*x)/4)/(x*x*x/3) + y = ratio*DCMR*DCMR*DCMR/3 + return y + y = x*x + # statement below fixed 13-Aug-03 to correct sign error in expansion + if (WK < 0): y = -y + ratio = 1 + y/5 + (2/105)*y*y + y= ratio*DCMR*DCMR*DCMR/3 + return y + + + h = H0/100 + WR = 4.165E-5/(h*h) # includes 3 massless neutrino species, T0 = 2.72528 + WK = 1-WM-WR-WV + az = 1.0/(1+1.0*z) + age = 0; + for i in range (0, n+1): + a = az*(i+0.5)/n + adot = sqrt(WK+(WM/a)+(WR/(a*a))+(WV*a*a)) + age = age + 1/adot + zage = az*age/n + # correction for annihilations of particles not present now like e+/e- + # added 13-Aug-03 based on T_vs_t.f + lpz = log((1+1.0*z))/log(10.0) + dzage = 0 + if (lpz > 7.500): dzage = 0.002 * (lpz - 7.500) + if (lpz > 8.000): dzage = 0.014 * (lpz - 8.000) + 0.001 + if (lpz > 8.500): dzage = 0.040 * (lpz - 8.500) + 0.008 + if (lpz > 9.000): dzage = 0.020 * (lpz - 9.000) + 0.028 + if (lpz > 9.500): dzage = 0.019 * (lpz - 9.500) + 0.039 + if (lpz > 10.000): dzage = 0.048 + if (lpz > 10.775): dzage = 0.035 * (lpz - 10.775) + 0.048 + if (lpz > 11.851): dzage = 0.069 * (lpz - 11.851) + 0.086 + if (lpz > 12.258): dzage = 0.461 * (lpz - 12.258) + 0.114 + if (lpz > 12.382): dzage = 0.024 * (lpz - 12.382) + 0.171 + if (lpz > 13.055): dzage = 0.013 * (lpz - 13.055) + 0.188 + if (lpz > 14.081): dzage = 0.013 * (lpz - 14.081) + 0.201 + if (lpz > 15.107): dzage = 0.214 + zage = zage*10.0**dzage + + zage_Gyr = (Tyr/H0)*zage + DTT = 0.0 + DCMR = 0.0 + # do integral over a=1/(1+z) from az to 1 in n steps, midpoint rule + for i in range (0,n+1): + a = az+(1-az)*(i+0.5)/n + adot = sqrt(WK+(WM/a)+(WR/(a*a))+(WV*a*a)) + DTT = DTT + 1/adot + DCMR = DCMR + 1/(a*adot) + + DTT = (1-az)*DTT/n + DCMR = (1-az)*DCMR/n + age = DTT+zage + age_Gyr = age*(Tyr/H0) + DTT_Gyr = (Tyr/H0)*DTT + DCMR_Gyr = (Tyr/H0)*DCMR + DCMR_Mpc = (c/H0)*DCMR + DA = az*DCMT() + DA_Mpc = (c/H0)*DA + kpc_DA = DA_Mpc/206.264806 + DA_Gyr = (Tyr/H0)*DA + DL = DA/(az*az) + DL_Mpc = (c/H0)*DL + DL_Gyr = (Tyr/H0)*DL + V_Gpc = 4*pi*pow(0.001*c/H0,3)*VCM() + + return DL_Mpc + + def ShiftToZero(DF, trange): """Given a range [.,.] shift the wave according to that intervall to zero""" s = 0 @@ -427,7 +565,7 @@ res[-1].append([]) return res

-def LoadModeArray(basename, lmax, col_time=1, col_re=2, col_im=3, filenameformat="cce_conversionscript"): +def LoadModeArray(basename, lmax, col_time=1, col_re=2, col_im=3, filenameformat="cce_conversion_script"): """Loads an array of modes up to lmax from disk. We assume that files begin with 'basename' and end with _l*m*.dat""" res = InitModeArray(lmax) @@ -551,7 +689,71 @@ return S

+def SensitivityeLISA(f): + # Noise curve for eLISA-NGO. + # Function from eLISA yellow book, Sec. 4.3

+ f0 = 150.0e-3 + a = 0.41 + T = sqrt(1.0+(f/(a*f0))**2) + + Da0 = 3.0 + fL = 0.1e-3 + fH = 8.0e-3 + da = Da0*1.0e-15 * sqrt(1.0+(f/fH)**4) * sqrt(1.0+fL/f) + + dx_DRS = 2.0*da / (2.0*pi*f)**2 + + Dx0 = 12.0 + f0 = 2.8e-3 + dx_IMS = Dx0*1e-12 * sqrt(1.0+(f0/f)**4) + + # Is this correct??? No... + S_IMS = dx_IMS + S_DRS = dx_DRS + + # arm length + L = 1e9 + + #h = sqrt(5.0)*2.0/sqrt(3.0)*T*sqrt(S_IMS + S_DRS) / L + # Use this in terms of dx_IMS, dx_DRS: + h = sqrt(5.0)*2.0/sqrt(3.0)*T*(dx_IMS + dx_DRS) / L + + S = h**2 + + return S + + +def SensitivityDECIGO(f): + """ From arXiv:1110.2865 + For this sensitivity curve, you don't need to select the + or x polarization component since the detector is a triangle. + The 60 degree abgle between detector arms is incoproated into the noise curve. + The antenna pattern is incorporated in the noise curve. + For SNRs, use an additional factor of sqrt(4) since BBO consists of 4 detectors!""" + + # White-dwarf confusion noise + SWD = 4.2e-47 * f**(-7.0/3.0) * exp(-2.0 * (f/5e-2)**2) + + fc = 7.69 + S = 3.30 * 1e-50*1.0/f**4 + 3.09*1e-47 * (1.0 + (f/fc)**2) + SWD + return S + + +def SensitivityBBO(f): + """ From arXiv:1110.2865 + For this sensitivity curve, you don't need to select the + or x polarization component since the detector is a triangle. + The 60 degree abgle between detector arms is incoproated into the noise curve. + The antenna pattern is incorporated in the noise curve. + For SNRs, use an additional factor of sqrt(4) since BBO consists of 4 detectors!""" + + # White-dwarf confusion noise + SWD = 4.2e-47 * f**(-7.0/3.0) * exp(-2.0 * (f/5e-2)**2) + + S = 6.15*1e-51*(1.0/f)**4 + 1.95*1e-48 + 1.20*1e-48*f**2 + SWD + return S + + + def SNR(fhSI, d, Sn, fLow, fHi): """Given a reconstructed h in SI units in Fourier domain and distance d (in Mpc) this computes the SNR for sensitivity curve Sn""" """Remember to select either h_+ or h_\cross! Otherise the SNR will be twice as large!""" @@ -561,6 +763,17 @@ return rho

+def SNRavg(fhSIModes, lmax, d, Sn, fLow, fHi): + """Given a modes array in SI units in Fourier domain and distance d (in Mpc) this computes the angle averaged SNR for sensitivity curve Sn""" + S = map(Sn, fhSIModes[2][2].x) # create table of values for noise curve + + rho = 0.0 + for ll in range(2, lmax+1): + for mm in range(-ll, ll+1): + rho += 1.0/pi * WaveFunction(fhSIModes[ll][mm].x, abs(fhSIModes[ll][mm].f/(d*Mpc))**2 / S).Integrate(fLow, fHi, False).real + return sqrt(rho) + + def IntCharFreq(fhSI, d, Sn, fLow, fHi): """Given a reconstructed h in SI units in Fourier domain and distance d (in Mpc) this computes the integrated characteristic signal frequency for sensitivity curve Sn""" """Remember to select either h_+ or h_\cross! Otherise the SNR will be twice as large!"""