[Commits] [svn:einsteintoolkit] pyGWAnalysis/trunk/DataAnalysis/ (Rev. 17)

reisswig at tapir.caltech.edu reisswig at tapir.caltech.edu
Sun Apr 28 11:56:56 CDT 2013


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!"""



More information about the Commits mailing list