newtonian noise simplifications

Simplify newtonian noise calculations to return displacement noise for
simplified configuration: either single test mass or single arm.

gwinc/ifo/noises.py

@@ -3,6 +3,7 @@ import numpy as np
 from .. import nb
 from .. import noise
 
+
 class QuantumVacuum(nb.Noise):
     """Quantum Vacuum
 
@@ -50,7 +51,8 @@ class Newtonian(nb.Noise):
     )
 
     def calc(self):
-        return noise.newtonian.gravg(self.freq, self.ifo)
+        n = noise.newtonian.gravg(self.freq, self.ifo.Seismic)
+        return n * 4 * self.ifo.gwinc.dhdl_sqr
 
 
 class NewtonianRayleigh(nb.Noise):
@@ -63,7 +65,8 @@ class NewtonianRayleigh(nb.Noise):
     )
 
     def calc(self):
-        return noise.newtonian.gravg_rayleigh(self.freq, self.ifo)
+        n = noise.newtonian.gravg_rayleigh(self.freq, self.ifo.Seismic)
+        return n * 2 * self.ifo.gwinc.dhdl_sqr
 
 
 class NewtonianBody(nb.Noise):
@@ -76,8 +79,9 @@ class NewtonianBody(nb.Noise):
     )
 
     def calc(self):
-        return (noise.newtonian.gravg_pwave(self.freq, self.ifo)
-               + noise.newtonian.gravg_swave(self.freq, self.ifo))
+        np = noise.newtonian.gravg_pwave(self.freq, self.ifo.Seismic)
+        ns = noise.newtonian.gravg_swave(self.freq, self.ifo.Seismic)
+        return (np + ns) * 4 * self.ifo.gwinc.dhdl_sqr
 
 
 class NewtonianInfrasound(nb.Noise):
@@ -90,7 +94,8 @@ class NewtonianInfrasound(nb.Noise):
     )
 
     def calc(self):
-        return noise.newtonian.atmois(self.freq, self.ifo)
+        n = noise.newtonian.atmois(self.freq, self.ifo.Atmospheric, self.ifo.Seismic)
+        return n * 2 * self.ifo.gwinc.dhdl_sqr
 
 
 class SuspensionThermal(nb.Noise):

gwinc/noise/newtonian.py

-from __future__ import division
+'''Functions to calculate Newtonian noise
 
+'''
+from __future__ import division
 import numpy as np
+from numpy import pi, sqrt, exp, log10
 import scipy.integrate as scint
 
-from numpy import pi, sqrt, exp
 from .seismic import seismic_ground_NLNM
 from .. import const
 
 
-def gravg(f, ifo):
-    """Return estimate of newtonian noise contribribution to |h(f)|^2
+def gravg(f, seismic):
+    """Gravity gradient noise for single test mass
+
+    :f: frequency array in Hz
+    :seismic: gwinc Seismic struct
 
-    N = GRAVG(F, IFO) returns the gravity gradient (Newtonian) noise
-    contribution in strain^2/Hz for four mirrors.
+    :returns: displacement noise power spectrum at :f:, in meters
 
     References:
 
@@ -21,7 +25,7 @@ def gravg(f, ifo):
 
      Driggers and Harms 2011, ``Results of Phase 1 Newtonian Noise
      Measurements at the LIGO Sites,'' February-March 2011.  T1100237.
-     https://dcc.ligo.org/cgi-bin/private/DocDB/ShowDocument?docid=60064
+     https://dcc.ligo.org/LIGO-T1100237
 
     Written by Enrico Camagna (?)
 
@@ -31,20 +35,18 @@ def gravg(f, ifo):
 
     """
 
-    fk = ifo.Seismic.KneeFrequency
-    a = ifo.Seismic.LowFrequencyLevel
-    L = ifo.Infrastructure.Length
-    gamma = ifo.Seismic.Gamma
-    ggcst = const.G
-    rho = ifo.Seismic.Rho
+    fk = seismic.KneeFrequency
+    a = seismic.LowFrequencyLevel
+    gamma = seismic.Gamma
+    rho = seismic.Rho
     # factor to account for correlation between masses
     # and the height of the mirror above the ground
-    beta = ifo.Seismic.Beta
-    h = ifo.Seismic.TestMassHeight
-    c_rayleigh = ifo.Seismic.RayleighWaveSpeed
+    beta = seismic.Beta
+    h = seismic.TestMassHeight
+    c_rayleigh = seismic.RayleighWaveSpeed
 
-    if 'Omicron' in ifo.Seismic:
-        omicron = ifo.Seismic.Omicron
+    if 'Omicron' in seismic:
+        omicron = seismic.Omicron
     else:
         omicron = 1
 
@@ -53,14 +55,14 @@ def gravg(f, ifo):
 
     # modelization of seismic noise (vertical)
     ground = a*coeff + a*(1-coeff)*(fk/f)**2
-    if 'Site' in ifo.Seismic and ifo.Seismic.Site == 'LLO':
+    if 'Site' in seismic and seismic.Site == 'LLO':
         ground = a*coeff*(fk/f) + a*(1-coeff)*(fk/f)**2
 
     # effective GG spring frequency, with G gravitational
-    fgg = sqrt(ggcst * rho) / (2*pi)
+    fgg = sqrt(const.G * rho) / (2*pi)
 
     # fixed numerical factors, 5/9/06, PF
-    n = (beta*4*pi/L*(fgg**2/f**2)*ground)**2
+    n = (beta*2*pi*(fgg**2/f**2)*ground)**2
 
     # The following two terms are corrections due to Jan Harms
     # https://git.ligo.org/rana-adhikari/CryogenicLIGO/issues/45
@@ -72,25 +74,29 @@ def gravg(f, ifo):
     # Feedforward cancellation
     n /= (omicron**2)
 
-    return n * ifo.gwinc.sinc_sqr
+    return n
+
+
+def gravg_rayleigh(f, seismic):
+    """Gravity gradient noise for single arm cavity from seismic Rayleigh waves
 
+    :f: frequency array in Hz
+    :seismic: gwinc Seismic structure
+
+    :returns: displacement noise power spectrum at :f:, in meters
 
-def gravg_rayleigh(f, ifo):
-    """Gravity gradient noise for seismic Rayleigh waves
     Following Harms LRR: https://doi.org/10.1007/lrr-2015-3
 
     """
-    fk = ifo.Seismic.KneeFrequency
-    a = ifo.Seismic.LowFrequencyLevel
-    L = ifo.Infrastructure.Length
-    gamma = ifo.Seismic.Gamma
-    ggcst = const.G
-    rho = ifo.Seismic.Rho
-    h = ifo.Seismic.TestMassHeight
-    c_rayleigh = ifo.Seismic.RayleighWaveSpeed
-
-    if 'Omicron' in ifo.Seismic:
-        omicron = ifo.Seismic.Omicron
+    fk = seismic.KneeFrequency
+    a = seismic.LowFrequencyLevel
+    gamma = seismic.Gamma
+    rho = seismic.Rho
+    h = seismic.TestMassHeight
+    c_rayleigh = seismic.RayleighWaveSpeed
+
+    if 'Omicron' in seismic:
+        omicron = seismic.Omicron
     else:
         omicron = 1
 
@@ -99,41 +105,47 @@ def gravg_rayleigh(f, ifo):
 
     # modelization of seismic noise (vertical)
     ground = a*coeff + a*(1-coeff)*(fk/f)**2
-    if 'Site' in ifo.Seismic and ifo.Seismic.Site == 'LLO':
+    if 'Site' in seismic and seismic.Site == 'LLO':
         ground = a*coeff*(fk/f) + a*(1-coeff)*(fk/f)**2
 
     # Harms LRR eqs. 35, 96, and 98
     w = 2 * pi * f
     k = w / c_rayleigh
-    kP = w / ifo.Seismic.pWaveSpeed
-    kS = w / ifo.Seismic.sWaveSpeed
+    kP = w / seismic.pWaveSpeed
+    kS = w / seismic.sWaveSpeed
     qzP = sqrt(k**2 - kP**2)
     qzS = sqrt(k**2 - kS**2)
     zeta = sqrt(qzP / qzS)
 
     gnu = k * (1 - zeta) / (qzP - k * zeta)
 
-    n = 2 * (2 * pi * ggcst * rho * exp(-h * k) * gnu)**2 * ground**2 / w**4
+    n = (2 * pi * const.G * rho * exp(-h * k) * gnu)**2 * ground**2 / w**4
 
     n /= omicron**2
 
-    return n * ifo.gwinc.dhdl_sqr
+    return n
+
+
+def gravg_pwave(f, seismic):
+    """Gravity gradient noise for single test mass from seismic p-waves
+
+    :f: frequency array in Hz
+    :seismic: gwinc Seismic structure
 
+    :returns: displacement noise power spectrum at :f:, in meters
 
-def gravg_pwave(f, ifo):
-    """Gravity gradient noise for seismic p-waves
     Following Harms LRR: https://doi.org/10.1007/lrr-2015-3
 
     """
-    ggcst = const.G
-    cP = ifo.Seismic.pWaveSpeed
-    levelP = ifo.Seismic.pWaveLevel
+    cP = seismic.pWaveSpeed
+    levelP = seismic.pWaveLevel
+    tmheight = seismic.TestMassHeight
+    rho_ground = seismic.Rho
+
     kP = (2 * pi * f) / cP
 
-    rho_ground = ifo.Seismic.Rho
     psd_ground_pwave = (levelP * seismic_ground_NLNM(f))**2
 
-    tmheight = ifo.Seismic.TestMassHeight 
     xP = np.abs(kP * tmheight)
 
     if tmheight >= 0:
@@ -148,74 +160,81 @@ def gravg_pwave(f, ifo):
         height_supp_power = (3 / 4) * np.array([scint.quad(lambda th, x: np.sin(th)**3
                 * (2 - np.exp(-x * np.sin(th)))**2, 0, pi, args=(x,))[0]
                 for x in xP])
-    psd_gravg_pwave = ((2 * pi * ggcst * rho_ground)**2
+    psd_gravg_pwave = ((2 * pi * const.G * rho_ground)**2
             * psd_ground_pwave * height_supp_power)
-    psd_gravg_pwave *= 4 / (2 * pi * f)**4
-    return psd_gravg_pwave * ifo.gwinc.dhdl_sqr
+    psd_gravg_pwave /= (2 * pi * f)**4
+    return psd_gravg_pwave
+
+
+def gravg_swave(f, seismic):
+    """Gravity gradient noise for single test mass from seismic s-waves
 
+    :f: frequency array in Hz
+    :seismic: gwinc Seismic structure
+
+    :returns: displacement noise power spectrum at :f:, in meters
 
-def gravg_swave(f, ifo):
-    """Gravity gradient noise for seismic s-waves
     Following Harms LRR: https://doi.org/10.1007/lrr-2015-3
 
     """
-    ggcst = const.G
-    cS = ifo.Seismic.sWaveSpeed
-    levelS = ifo.Seismic.sWaveLevel
+    cS = seismic.sWaveSpeed
+    levelS = seismic.sWaveLevel
+    tmheight = seismic.TestMassHeight
+    rho_ground = seismic.Rho
+
     kS = (2 * pi * f) / cS
 
-    rho_ground = ifo.Seismic.Rho
     psd_ground_swave = (levelS * seismic_ground_NLNM(f))**2
 
-    tmheight = ifo.Seismic.TestMassHeight 
     xS = np.abs(kS * tmheight)
 
     # For both surface and underground facilities
     height_supp_power = (3 / 2) * np.array([scint.quad(lambda th, x: np.sin(th)**3
             * np.exp(-2 * x * np.sin(th)), 0, pi / 2, args=(x,))[0]
             for x in xS])
-    psd_gravg_swave = ((2 * pi * ggcst * rho_ground)**2
+    psd_gravg_swave = ((2 * pi * const.G * rho_ground)**2
             * psd_ground_swave * height_supp_power)
-    psd_gravg_swave *= 4 / (2 * pi * f)**4
+    psd_gravg_swave /= (2 * pi * f)**4
 
-    return psd_gravg_swave * ifo.gwinc.dhdl_sqr
+    return psd_gravg_swave
 
 
-def atmois(f, ifo):
-    """
-    Atmospheric infrasound calculation following Harms LRR.
+def atmois(f, atmos, seismic):
+    """Atmospheric infrasound newtonian noise for single arm cavity
 
-    """
-    p_air = ifo.Atmospheric.AirPressure
-    rho_air = ifo.Atmospheric.AirDensity
-    ai_air = ifo.Atmospheric.AdiabaticIndex
-    c_sound = ifo.Atmospheric.SoundSpeed
+    :f: frequency array in Hz
+    :atmos: gwinc Atmospheric structure
+    :seismic: gwinc Seismic structure
 
-    L = ifo.Infrastructure.Length
-    ggcst = const.G
-    h = ifo.Seismic.TestMassHeight
+    :returns: displacement noise power spectrum at :f:, in meters
+
+    """
+    p_air = atmos.AirPressure
+    rho_air = atmos.AirDensity
+    ai_air = atmos.AdiabaticIndex
+    c_sound = atmos.SoundSpeed
+    h = seismic.TestMassHeight
 
     w = 2 * pi * f
     k = w / c_sound
 
     # Pressure spectrum
     try:
-        a_if = ifo.Atmospheric.InfrasoundLevel1Hz
-        e_if = ifo.Atmospheric.InfrasoundExponent
+        a_if = atmos.InfrasoundLevel1Hz
+        e_if = atmos.InfrasoundExponent
         psd_if = (a_if * f**e_if)**2
     except AttributeError:
         psd_if = atmoBowman(f)**2
 
     # Harms LRR (2015), eq. 172
     # https://doi.org/10.1007/lrr-2015-3
-    # With an extra factor 2 for two arms
     # And with the Bessel terms ignored... for 4 km this amounts to a 10%
     # correction at 10 Hz and a 30% correction at 1 Hz
-    coupling_if = 4./3 * (4 * pi / (k * w**2) * ggcst * rho_air / (ai_air * p_air))**2
+    coupling_if = 2/3 * (4 * pi / (k * w**2) * const.G * rho_air / (ai_air * p_air))**2
 
     n_if = coupling_if * psd_if
 
-    return n_if * ifo.gwinc.dhdl_sqr
+    return n_if
 
 
 def atmoBowman(f):
@@ -223,12 +242,15 @@ def atmoBowman(f):
    
     """ 
     freq = np.array([
-                0.01, 0.0155, 0.0239, 0.0367, 0.0567,
-                0.0874, 0.1345, 0.2075, 0.32, 0.5,
-                0.76, 1.17, 1.8, 2.79, 4.3,
-                6.64, 10, 100])
-    pressure_asd = np.sqrt([22.8, 4, 0.7, 0.14, 0.027, 0.004,
-                            0.0029, 0.0039, 7e-4, 1.44e-4, 0.37e-4,
-                            0.12e-4, 0.56e-5, 0.35e-5, 0.26e-5, 0.24e-5,
-                            2e-6, 2e-6])
-    return 10**(np.interp(np.log10(f), np.log10(freq), np.log10(pressure_asd)))
+        0.01, 0.0155, 0.0239, 0.0367, 0.0567,
+        0.0874, 0.1345, 0.2075, 0.32, 0.5,
+        0.76, 1.17, 1.8, 2.79, 4.3,
+        6.64, 10, 100,
+    ])
+    pressure_asd = sqrt([
+        22.8, 4, 0.7, 0.14, 0.027, 0.004,
+        0.0029, 0.0039, 7e-4, 1.44e-4, 0.37e-4,
+        0.12e-4, 0.56e-5, 0.35e-5, 0.26e-5, 0.24e-5,
+        2e-6, 2e-6,
+    ])
+    return 10**(np.interp(log10(f), log10(freq), log10(pressure_asd)))