substrate brownian noise simplification

substrate_brownian() function updated to calculate thermal noise of a
single suspended optic.

ifo SubstrateBrownian Noise class modified to handle specific case of
Fabry-Perot Michelson.

gwinc/ifo/noises.py

@@ -202,7 +202,11 @@ class SubstrateBrownian(nb.Noise):
     )
 
     def calc(self):
-        return noise.substratethermal.subbrownian(self.freq, self.ifo)
+        nITM = noise.substratethermal.substrate_brownian(
+            self.freq, self.ifo.Materials, self.ifo.Optics.ITM.BeamRadius)
+        nETM = noise.substratethermal.substrate_brownian(
+            self.freq, self.ifo.Materials, self.ifo.Optics.ETM.BeamRadius)
+        return (nITM + nETM) * 2 * self.ifo.gwinc.dhdl_sqr
 
 
 class SubstrateThermoElastic(nb.Noise):

gwinc/noise/substratethermal.py

@@ -99,91 +99,89 @@ def substrate_thermorefractive(f, materials, wBeam, exact=False):
     else:
         psd = 4*H*beta**2*kappa*kBT*Temp/(pi*r0**4*omega**2*(rho*C)**2)
 
-def subbrownian(f, ifo):
-    """Strain noise from the Brownian thermal noise due to substrate mechanical loss
+    return psd
+
+
+def substrate_brownian(f, materials, wBeam):
+    """Substrate thermal displacement noise spectrum due to substrate mechanical loss
+
+    :f: frequency array in Hz
+    :materials: gwinc optic materials structure
+    :wBeam: beam radius (at 1 / e^2 power)
+
+    :returns: displacement noise power spectrum at :f:, in meters
 
     """
-    wITM = ifo.Optics.ITM.BeamRadius
-    wETM = ifo.Optics.ETM.BeamRadius
-    Y = ifo.Materials.Substrate.MirrorY
-    sigma = ifo.Materials.Substrate.MirrorSigma
+    Y = materials.Substrate.MirrorY
+    sigma = materials.Substrate.MirrorSigma
+    c2 = materials.Substrate.c2
+    n = materials.Substrate.MechanicalLossExponent
+    alphas = materials.Substrate.Alphas
+    kBT = const.kB * materials.Substrate.Temp
 
-    c2 = ifo.Materials.Substrate.c2
-    n = ifo.Materials.Substrate.MechanicalLossExponent
-    alphas = ifo.Materials.Substrate.Alphas
-    L = ifo.Infrastructure.Length
-    kBT = const.kB * ifo.Materials.Substrate.Temp
+    cftm, aftm = substrate_brownian_FiniteCorr(materials, wBeam)
 
     # Bulk substrate contribution
     phibulk = c2 * f**n
-
-    cITM, aITM = subbrownianFiniteCorr(ifo, 'ITM')
-    cETM, aETM = subbrownianFiniteCorr(ifo, 'ETM')
-    cbulk = 8 * kBT * (aITM + aETM) * phibulk / (2 * pi * f)
+    cbulk = 8 * kBT * aftm * phibulk / (2 * pi * f)
 
     # Surface loss contribution
-    # csurfETM = alphas/(Y*pi*wETM^2);
-    # csurfITM = alphas/(Y*pi*wITM^2);
+    # csurf = alphas/(Y*pi*wBeam^2)
+    csurf = alphas*(1-2*sigma)/((1-sigma)*Y*pi*wBeam**2)
+    csurf *= 8 * kBT / (2 * pi * f)
 
-    csurfETM = alphas*(1-2*sigma)/((1-sigma)*Y*pi*wETM**2)
-    csurfITM = alphas*(1-2*sigma)/((1-sigma)*Y*pi*wITM**2)
-    csurf = 8 * kBT * (csurfITM + csurfETM) / (2 * pi * f)
+    return csurf + cbulk
 
-    # account for 2 ITM and 2 ETM, and convert to strain whith 1/L^2
-    n = 2 * (csurf + cbulk) * ifo.gwinc.dhdl_sqr
-
-    return n
 
+def substrate_brownian_FiniteCorr(materials, wBeam):
+    """Substrate brownian noise finite-size test mass correction
 
-def subbrownianFiniteCorr(ifo, opticName):
-    """Amplitude coefficient of mirror thermal noise
+    :materials: gwinc optic materials structure
+    :wBeam: beam radius (at 1 / e^2 power)
 
-    Contribution for finite-size test masses.
-    
-    [cftm, aftm] = subbrownianFiniteCorr(ifo, opticName)
+    :returns: correction factors tuple:
     cftm = finite mirror correction factor
     aftm = amplitude coefficient for thermal noise:
            thermal noise contribution to displacement noise is
            S_x(f) = (8 * kB * T / (2*pi*f)) * Phi(f) * aftm
-    
+
     Equation references to Bondu, et al. Physics Letters A 246 (1998)
     227-236 (hereafter BHV) and Liu and Thorne gr-qc/0002055 (hereafter LT)
 
     """
-    # get some numbers
-    a = ifo.Materials.MassRadius
-    h = ifo.Materials.MassThickness
-    w = ifo.Optics[opticName].BeamRadius
-    Y = ifo.Materials.Substrate.MirrorY
-    sigma = ifo.Materials.Substrate.MirrorSigma
+    a = materials.MassRadius
+    h = materials.MassThickness
+    Y = materials.Substrate.MirrorY
+    sigma = materials.Substrate.MirrorSigma
 
-    # do the work
-    r0 = w / sqrt(2)                     # LT uses e-folding of power
+    # LT uses e-folding of power
+    r0 = wBeam / sqrt(2)
     km = zeta/a
 
-    Qm = exp(-2*km*h)                    # LT eq. 35a
+    Qm = exp(-2*km*h) # LT eq. 35a
 
     Um = (1-Qm)*(1+Qm)+4*h*km*Qm
-    Um = Um/((1-Qm)**2-4*(km*h)**2*Qm)   # LT 53 (BHV eq. btwn 29 & 30)
+    Um = Um/((1-Qm)**2-4*(km*h)**2*Qm) # LT 53 (BHV eq. btwn 29 & 30)
 
     x = exp(-(zeta*r0/a)**2/4)
-    s = sum(x/(zeta**2*j0m))             # LT 57
+    s = sum(x/(zeta**2*j0m)) # LT 57
 
     x2 = x*x
     U0 = sum(Um*x2/(zeta*j0m**2))
     U0 = U0*(1-sigma)*(1+sigma)/(pi*a*Y) # LT 56 (BHV eq. 3)
 
-    p0 = 1/(pi*a**2)                     # LT 28
+    p0 = 1/(pi*a**2) # LT 28
     DeltaU = (pi*h**2*p0)**2
     DeltaU = DeltaU + 12*pi*h**2*p0*sigma*s
     DeltaU = DeltaU + 72*(1-sigma)*s**2
-    DeltaU = DeltaU*a**2/(6*pi*h**3*Y)   # LT 54
+    DeltaU = DeltaU*a**2/(6*pi*h**3*Y) # LT 54
 
-    aftm = DeltaU + U0                   # LT 58 (eq. following BHV 31)
+    # LT 58 (eq. following BHV 31)
+    aftm = DeltaU + U0
 
-    # amplitude coef for infinite TM
-    #   factored out: (8 * kB * T * Phi) / (2 * pi * f)
-    aitm = (1 - sigma**2) / (2 * sqrt(2 * pi) * Y * r0) # LT 59
+    # amplitude coef for infinite TM, LT 59
+    # factored out: (8 * kB * T * Phi) / (2 * pi * f)
+    aitm = (1 - sigma**2) / (2 * sqrt(2 * pi) * Y * r0)
 
     # finite mirror correction
     cftm = aftm / aitm