newtonian: update Newtonian noise calculation

Implements the changes recommended by Jan Harms in https://git.ligo.org/rana-adhikari/CryogenicLIGO/issues/45
Fixes #18

gwinc/ifo/A+.yaml

@@ -58,11 +58,13 @@ Seismic:
   Site: 'LHO'                     # LHO or LLO (only used for Newtonian noise)
   KneeFrequency: 10               # Hz; freq where 'flat' noise rolls off
   LowFrequencyLevel: 1e-9         # m/rtHz; seismic noise level below f_knee
-  Gamma: .8                       # abruptness of change at f_knee
+  Gamma: 0.8                      # abruptness of change at f_knee
   Rho: 1.8e3                      # kg/m^3; density of the ground nearby
-  Beta: 0.6                       # quiet times beta: 0.35-0.60
-  # noisy times beta: 0.15-1.4
+  Beta: 0.8                       # quiet times beta: 0.35-0.60
+                                  # noisy times beta: 0.15-1.4
   Omicron: 1                      # Feedforward cancellation factor
+  TestMassHeight: 1.5             # m
+  RayleighWaveSpeed: 250          # m/s
 
 Suspension:
   Type: 'Quad'

gwinc/ifo/CEcryo.yaml

@@ -72,10 +72,12 @@ Seismic:
   Site: 'LHO'                      # LHO or LLO (only used for Newtonian noise)
   KneeFrequency: 10                # Hz; freq where 'flat' noise rolls off
   LowFrequencyLevel: 1e-9          # m/rtHz; seismic noise level below f_knee
-  Gamma: .8                        # abruptness of change at f_knee
+  Gamma: 0.8                       # abruptness of change at f_knee
   Rho: 1.8e3                       # kg/m^3; density of the ground nearby
-  Beta: 1                          # quiet times beta = 0.35-0.60; noisy times beta = 0.15-1.4
+  Beta: 0.8                        # quiet times beta = 0.35-0.60; noisy times beta = 0.15-1.4
   Omicron: 10                      # Feedforward cancellation factor
+  TestMassHeight: 1.5              # m
+  RayleighWaveSpeed: 250           # m/s
   #darmSeiSusFile: 'CryogenicLIGO/Sensitivity/GWINC/seismic.mat'
 
 Suspension:

gwinc/ifo/Voyager.yaml

@@ -70,10 +70,12 @@ Seismic:
   Site: 'LHO'                      # LHO or LLO (only used for Newtonian noise)
   KneeFrequency: 10                # Hz; freq where 'flat' noise rolls off
   LowFrequencyLevel: 1e-9          # m/rtHz; seismic noise level below f_knee
-  Gamma: .8                        # abruptness of change at f_knee
+  Gamma: 0.8                       # abruptness of change at f_knee
   Rho: 1.8e3                       # kg/m^3; density of the ground nearby
-  Beta: 1                          # quiet times beta = 0.35-0.60; noisy times beta = 0.15-1.4
+  Beta: 0.8                        # quiet times beta = 0.35-0.60; noisy times beta = 0.15-1.4
   Omicron: 10                      # Feedforward cancellation factor
+  TestMassHeight: 1.5              # m
+  RayleighWaveSpeed: 250           # m/s
   #darmSeiSusFile: 'CryogenicLIGO/Sensitivity/GWINC/seismic.mat'
 
 Suspension:

gwinc/ifo/aLIGO.yaml

@@ -58,11 +58,13 @@ Seismic:
   Site: 'LHO'                     # LHO or LLO (only used for Newtonian noise)
   KneeFrequency: 10               # Hz; freq where 'flat' noise rolls off
   LowFrequencyLevel: 1e-9         # m/rtHz; seismic noise level below f_knee
-  Gamma: .8                       # abruptness of change at f_knee
+  Gamma: 0.8                      # abruptness of change at f_knee
   Rho: 1.8e3                      # kg/m^3; density of the ground nearby
-  Beta: 0.5                       # quiet times beta: 0.35-0.60
-  # noisy times beta: 0.15-1.4
+  Beta: 0.8                       # quiet times beta: 0.35-0.60
+                                  # noisy times beta: 0.15-1.4
   Omicron: 1                      # Feedforward cancellation factor
+  TestMassHeight: 1.5             # m
+  RayleighWaveSpeed: 250          # m/s
 
 Suspension:
   Type: 'Quad'

gwinc/noise/newtonian.py

 from __future__ import division
 import scipy.constants
-from numpy import pi, sqrt
-import logging
+from numpy import pi, sqrt, exp
 
 
 def gravg(f, ifo):
@@ -9,100 +8,63 @@ def gravg(f, ifo):
 
     N = GRAVG(F, IFO) returns the gravity gradient (Newtonian) noise
     contribution in strain^2/Hz for four mirrors.
-    
+
     References:
-    
+
      Saulson 1984,           http://dx.doi.org/10.1103/PhysRevD.30.732
      Hughes and Thorne 1998, http://dx.doi.org/10.1103/PhysRevD.58.122002
-    
+
      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
-    
-    See also: ground.m
-    
+
     Written by Enrico Camagna (?)
-    
+
     added to Bench by Gregg Harry 8/27/03
     seismic spectrum modified by Jan Harms 05/11/2010
     Calculates gravity gradient noise for four mirrors
 
     """
-    L     = ifo.Infrastructure.Length
-    G     = scipy.constants.G
-    rho   = ifo.Seismic.Rho
+
+    fk = ifo.Seismic.KneeFrequency
+    a = ifo.Seismic.LowFrequencyLevel
+    L = ifo.Infrastructure.Length
+    gamma = ifo.Seismic.Gamma
+    ggcst = scipy.constants.G
+    rho = ifo.Seismic.Rho
     # factor to account for correlation between masses
     # and the height of the mirror above the ground
-    beta  = ifo.Seismic.Beta
+    beta = ifo.Seismic.Beta
+    h = ifo.Seismic.TestMassHeight
+    c_rayleigh = ifo.Seismic.RayleighWaveSpeed
 
-    if 'Spectrum' in ifo.Seismic:
-        seismic = ifo.Seismic.Spectrum
+    if 'Omicron' in ifo.Seismic:
+        omicron = ifo.Seismic.Omicron
     else:
-        seismic = ground(ifo.Seismic, f)
-
-    # effective GG spring frequency, with G gravitational
-    fgg = sqrt(G * rho) / (2*pi)
-
-    # fixed numerical factors, 5/9/06, PF
-    n = (beta*4*pi/L*(fgg**2/f**2)*seismic)**2
-
-    # Feedforward cancellation
-    n /= (ifo.Seismic.Omicron)**2
-
-    return n
-
-
-def ground(Seismic, f):
-    """Estimate of seismic displacement spectrum
-
-    N = GROUND(SEISMIC, F) returns the estimated ground (seismic)
-    displacement spectrum in meters / rtHz at frequencies F using
-    parameters SEISMIC (usually a field of the IFOMODEL).
-    
-    Example usage:
-    
-        ifo = IFOModel();
-        f = logspace(log10(1), log10(60), 101);
-        n = ground(ifo.Seismic, f);
-        loglog(f, n)
-        ylabel('meters / rtHz');
-        xlabel('Hz');
-    
-    This function is currently only used by gravg.m (Newtonian Noise).
-    
-    Jan says: The seismic NN is 90th percentile. But this was just to
-    lift the spectra above the aLIGO GWINC sensitivity. 90th, 95th,
-    whatever percentile does not matter since you want to guarantee
-    that NN is not limiting at any time. So in the spectral plot the
-    only important information is that NN can be high enough to be
-    seen in the detector.
-    
-    References:
-      Waldman and Fritschel, 'Reference Seismic Data for LLO', T0900312
-      https://dcc.ligo.org/cgi-bin/private/DocDB/ShowDocument?docid=3315
-
-    """
-    fk = Seismic.KneeFrequency
-    a1 = Seismic.LowFrequencyLevel
-    a2 = a1*100
-    gamma = Seismic.Gamma
+        omicron = 1
 
     # a sort of theta function (Fermi distr.)
     coeff = 1/(1 + 3**(gamma*(f-fk)))
 
-    # modelization of seismic noise (velocity)
-    if 'Site' not in Seismic:
-        logging.info('defaulting to Livingston site')
-        Seismic.Site = 'LLO'
+    # modelization of seismic noise (vertical)
+    ground = a*coeff + a*(1-coeff)*(fk/f)**2
+    if 'Site' in ifo.Seismic and ifo.Seismic.Site == 'LLO':
+        ground = a*coeff*(fk/f) + a*(1-coeff)*(fk/f)**2
 
-    if Seismic.Site == 'LHO':
-        n = (2*pi*f)**(4/3)*(a1*coeff + a1*(1-coeff)*(fk/f)**(9/3))
-    elif Seismic.Site == 'LLO':
-        n = a2*coeff + a2*(1-coeff)*(fk/f)**2
-    else:
-        raise 'Unknown seismic.site'
+    # effective GG spring frequency, with G gravitational
+    fgg = sqrt(ggcst * rho) / (2*pi)
+
+    # fixed numerical factors, 5/9/06, PF
+    n = (beta*4*pi/L*(fgg**2/f**2)*ground)**2
 
-    # convert into displacement
-    n /= 2*pi*f
+    # The following two terms are corrections due to Jan Harms
+    # https://git.ligo.org/rana-adhikari/CryogenicLIGO/issues/45
+    # (1) projection of NN force onto the direction of the arm
+    n = n * 1/2
+    # (2) exponential cutoff at frequency (seismic speed)/(test mass height)
+    n = n * exp(-4*pi*f*h/c_rayleigh)
+
+    # Feedforward cancellation
+    n /= (omicron**2)
 
-    return n
+    return n * ifo.gwinc.sinc_sqr