Merge branch 'quantum-susceptibility' into 'master'

Radiation pressure calculations use suspension instead of free mass susceptibility

Closes #71

See merge request gwinc/pygwinc!114

gwinc/ifo/noises.py

@@ -132,19 +132,21 @@ def precomp_suspension(f, ifo):
         pc.VHCoupling.theta = ifo.Suspension.VHCoupling.theta
     else:
         pc.VHCoupling.theta = ifo.Infrastructure.Length / const.R_earth
-    hForce, vForce, hTable, vTable = suspension.suspQuad(f, ifo.Suspension)
+    hForce, vForce, hTable, vTable, tst_suscept = suspension.suspQuad(
+        f, ifo.Suspension)
     pc.hForce = hForce
     pc.vForce = vForce
     pc.hTable = hTable
     pc.vTable = vTable
+    pc.tst_suscept = tst_suscept
     return pc
 
 
-def precomp_quantum(f, ifo):
+@nb.precomp(sustf=precomp_suspension)
+def precomp_quantum(f, ifo, sustf):
     pc = Struct()
-    mirror_mass = mirror_struct(ifo, 'ETM').MirrorMass
     power = ifo_power(ifo)
-    noise_dict = noise.quantum.shotrad(f, ifo, mirror_mass, power)
+    noise_dict = noise.quantum.shotrad(f, ifo, sustf, power)
     pc.ASvac = noise_dict['ASvac']
     pc.SEC = noise_dict['SEC']
     pc.Arm = noise_dict['arm']

gwinc/noise/quantum.py

@@ -90,12 +90,12 @@ def getSqzParams(ifo):
 ######################################################################
 
 
-def shotrad(f, ifo, mirror_mass, power):
+def shotrad(f, ifo, sustf, power):
     """Quantum noise strain spectrum
 
     :f: frequency array in Hz
     :ifo: gwinc IFO Struct
-    :mirror_mass: mirror mass in kg
+    :sustf: suspension transfer function struct
     :power: gwinc power Struct
 
     :returns: displacement noise power spectrum at :f:
@@ -120,12 +120,12 @@ def shotrad(f, ifo, mirror_mass, power):
     ######################################################################
 
     # call the main IFO Quantum Model
-    if 'Type' not in ifo.Optics:
-        fname = shotradSignalRecycled
+    if 'Type' not in ifo.Optics or ifo.Optics.Type == 'SignalRecycled':
+        coeff, Mifo, Msig, Mn = shotradSignalRecycled(f, ifo, sustf, power)
+    elif ifo.Optics.Type == 'SignalRecycledBnC':
+        coeff, Mifo, Msig, Mn = shotradSignalRecycledBnC(f, ifo, power)
     else:
-        namespace = globals()
-        fname = namespace['shotrad' + ifo.Optics.Type]
-    coeff, Mifo, Msig, Mn = fname(f, ifo, mirror_mass, power)
+        raise ValueError('Unrecognized IFO type ' + ifo.Optics.Type)
 
     # separate arm and SEC loss
     Mnoise = dict(arm=Mn[:, :2, :], SEC=Mn[:, 2:, :])
@@ -221,7 +221,7 @@ def compile_SEC_RES_TF():
 # Main IFO quantum models
 ######################################################################
 
-def shotradSignalRecycled(f, ifo, mirror_mass, power):
+def shotradSignalRecycled(f, ifo, sustf, power):
     """Quantum noise model for signal recycled IFO (see shotrad for more info)
 
     New version July 2016 by JH based on transfer function formalism
@@ -242,7 +242,7 @@ def shotradSignalRecycled(f, ifo, mirror_mass, power):
     L = ifo.Infrastructure.Length                # Length of arm cavities [m]
     l = ifo.Optics.SRM.CavityLength              # SRC Length [m]
     T = ifo.Optics.ITM.Transmittance             # ITM Transmittance [Power]
-    m = mirror_mass                              # Mirror mass [kg]
+    m = ifo.Suspension.Stage[0].Mass             # Mirror mass [kg]
 
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     bsloss = ifo.Optics.BSLoss                   # BS Loss [Power]
@@ -298,7 +298,7 @@ def shotradSignalRecycled(f, ifo, mirror_mass, power):
 
     Omega = 2*pi*f                               # Signal angular frequency [rad/s]
     h_SQL = sqrt(8 * hbar / (m * (Omega * L)**2)) # SQL Strain
-    K = 16 * P * omega_0 / (m * c**2 * Omega**2)
+    K = - 16 * P * omega_0 / c**2 * sustf.tst_suscept
 
     # arm cavity
     exp_2jOmegaL_c = exp(2j*Omega*L/c)
@@ -391,7 +391,7 @@ def shotradSignalRecycled(f, ifo, mirror_mass, power):
     return coeff, Mifo, Msig, Mnoise
 
 
-def shotradSignalRecycledBnC(f, ifo, mirror_mass, power):
+def shotradSignalRecycledBnC(f, ifo, power):
     """Quantum noise model for signal recycled IFO
 
     See shotrad for more info.
@@ -423,7 +423,7 @@ def shotradSignalRecycledBnC(f, ifo, mirror_mass, power):
     L = ifo.Infrastructure.Length                # Length of arm cavities [m]
     l = ifo.Optics.SRM.CavityLength              # SRC Length [m]
     T = ifo.Optics.ITM.Transmittance             # ITM Transmittance [Power]
-    m = mirror_mass                              # Mirror mass [kg]
+    m = ifo.Suspension.Stage[0].Mass             # Mirror mass [kg]
 
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     bsloss = ifo.Optics.BSLoss                   # BS Loss [Power]

gwinc/squeeze.py

@@ -9,8 +9,7 @@ def sql(ifo):
     c = const.c
     power = ifo_power(ifo)
     w0 = 2 * pi * c / ifo.Laser.Wavelength
-    rho = ifo.Materials.Substrate.MassDensity
-    m = pi * ifo.Materials.MassRadius**2 * ifo.Materials.MassThickness * rho
+    m = ifo.Suspension.Stage[0].Mass
     Titm = ifo.Optics.ITM.Transmittance
     Tsrm = ifo.Optics.SRM.Transmittance
     tSR = sqrt(Tsrm)

gwinc/suspension.py

@@ -587,9 +587,13 @@ def suspQuad(f, sus):
 
     # Calculate horizontal/vertical TFs with all losses on,
     # for the vibration isolation calculation
-    k = khr + 1j*(khi[:,0,:] + khi[:,1,:])
-    hTable = top_displ_to_tst_displ(k, masses, f)
     k = kvr + 1j*(kvi[:,0,:] + kvi[:,1,:])
     vTable = top_displ_to_tst_displ(k, masses, f)
 
-    return hForce, vForce, hTable, vTable
+    k = khr + 1j*(khi[:,0,:] + khi[:,1,:])
+    hTable = top_displ_to_tst_displ(k, masses, f)
+
+    # test mass susceptibility for radiation pressure calculations
+    tst_suscept = tst_force_to_tst_displ(k, masses, f)
+
+    return hForce, vForce, hTable, vTable, tst_suscept