More accurate displacement to strain conversion

Use function from R. Schilling CQG 14 1997.  Particularly needed for longer
arms.

gwinc/gwinc.py

@@ -7,6 +7,7 @@ import logging
 
 from .precomp import precompIFO
 from . import suspension
+from . import util
 from . import noise
 from . import plot
 
@@ -40,6 +41,13 @@ def noise_calc(ifo, f):
     ifo.Suspension.hTable = hTable
     ifo.Suspension.vTable = vTable
 
+    ##############################
+    # strain to length conversion
+    #
+    # used for converting displacement to strain in all noise calculations
+
+    ifo.gwinc.dhdl_sqr, ifo.gwinc.sinc_sqr = util.dhdl(f, ifo.Infrastructure.Length)
+
     ##############################
 
     noises = OrderedDict()

gwinc/noise/coatingthermal.py

@@ -26,7 +26,7 @@ def coatbrownian(f, ifo):
     SbrITM  = getCoatBrownian(f, ifo, wBeam_ITM, dOpt_ITM)
     SbrETM  = getCoatBrownian(f, ifo, wBeam_ETM, dOpt_ETM)
 
-    n = 2 * (SbrITM + SbrETM) / Length**2
+    n = 2 * (SbrITM + SbrETM) * ifo.gwinc.dhdl_sqr
 
     return n
 
@@ -135,7 +135,8 @@ def thermooptic(f, ifo):
     StoITM, junk1, junk2, junk3  = getCoatThermoOptic(f, ifo, wBeam_ITM, dOpt_ITM[:])
     StoETM, junk1, junk2, junk3  = getCoatThermoOptic(f, ifo, wBeam_ETM, dOpt_ETM[:])
   
-    n = 2 * (StoITM + StoETM) / Length**2
+    n = 2 * (StoITM + StoETM) * ifo.gwinc.dhdl_sqr
+
     return n
 
 

gwinc/noise/seismic.py

@@ -28,9 +28,9 @@ def seismic(f, ifo):
     n = nv + nh
 
     # Convert into Strain PSD (4 TMs)
-    nh *= 4 / ifo.Infrastructure.Length**2
-    nv *= 4 / ifo.Infrastructure.Length**2
-    n  *= 4 / ifo.Infrastructure.Length**2
+    nh *= 4 * ifo.gwinc.dhdl_sqr
+    nv *= 4 * ifo.gwinc.dhdl_sqr
+    n  *= 4 * ifo.gwinc.dhdl_sqr
 
     return n, nh, nv
 

gwinc/noise/substratethermal.py

@@ -133,7 +133,8 @@ def subbrownian(f, ifo):
     csurf = 8 * kBT * (csurfITM + csurfETM) / (2 * pi * f)
 
     # account for 2 ITM and 2 ETM, and convert to strain whith 1/L^2
-    n = 2 * (csurf + cbulk) / L**2
+    n = 2 * (csurf + cbulk) * ifo.gwinc.dhdl_sqr
+
     return n
 
 
@@ -222,7 +223,8 @@ def subtherm(f, ifo):
     SETM = SETM * subthermFiniteCorr(ifo, 'ETM')
 
     # 2 mirrors each type, factor omega^2, dimensionless strain
-    n = 2*(SITM + SETM)/(2*pi*f*L)**2
+    n = 2 * (SITM + SETM)/(2*pi*f)**2 * ifo.gwinc.dhdl_sqr
+
     return n
 
 

gwinc/noise/suspensionthermal.py

@@ -68,5 +68,6 @@ def susptherm(f, ifo):
             noise += 4 * kB * Temp[ii] * abs(imag(dxdF[ii,:])) / w
 
     # 4 masses, turn into gravitational wave strain
-    noise *= 4 / ifo.Infrastructure.Length**2
+    noise *= 4 * ifo.gwinc.dhdl_sqr
+
     return np.squeeze(noise)

gwinc/util.py

+from numpy import pi, sin, exp, abs
+import scipy.constants
+
+
+def dhdl(f, armlen):
+    """Strain to length conversion for noise power spetra
+
+    This takes into account the GW wavelength and is only important
+    when this is comparable to the detector arm length.
+
+    From R. Schilling, CQG 14 (1997) 1513-1519, equation 5,
+    with n = 1, nu = 0.05, ignoring overall phase and cos(nu)^2.
+    A small value of nu is used instead of zero to avoid infinities.
+
+    Returns the square of the dh/dL function, and the same divided by
+    the arm length squared.
+
+    """
+    c = scipy.constants.c
+    nu_small = 0.05
+    omega_arm = pi * f * armlen / c
+    omega_arm_f = (1 - sin(nu_small)) * pi * f * armlen / c
+    omega_arm_b = (1 + sin(nu_small)) * pi * f * armlen / c
+    sinc_sqr = 4 / abs(sin(omega_arm_f) * exp(-1j * omega_arm) / omega_arm_f
+                       + sin(omega_arm_b) * exp(1j * omega_arm) / omega_arm_b)**2
+    # keep DC value equal to 1
+    sinc_sqr /= sinc_sqr[0]
+    dhdl_sqr = sinc_sqr / armlen**2
+    return dhdl_sqr, sinc_sqr