quantum.py functions updated from matgwinc

The matgwinc updates include a new shotradSignalRecycled(), and minor
corrections to the previous code, which is still available as
shotradSignalRecycledBnC().

shotradSignalRecycled() has been vectorized for good performance,
with its two TF inversions pre-computed by sympy.

gwinc/noise/quantum.py

@@ -254,7 +254,209 @@ def shotrad(f, ifo):
     return n * ifo.gwinc.sinc_sqr
 
 
+def compile_ARM_RES_TF():
+    import sympy as sp
+    ID = sp.eye(2)
+    rITM, tArm, exp_2jOmegaL_c, K = sp.symbols('rITM tArm exp_2jOmegaL_c K')
+    ARM = tArm * exp_2jOmegaL_c * sp.Matrix([[1, 0], [-K, 1]])
+    ARM_RES = (ID - rITM*ARM)**-1
+
+    subexprs, ARM_RES_expr = sp.cse(ARM_RES)
+    for expr in subexprs:
+        print(str(expr[0]), '=', str(expr[1]))
+    print('RES', '=', str(ARM_RES_expr[0]).replace('Matrix', 'np.array').replace(', 0]', ', np.zeros(nf)]'))
+
+
+def compile_SEC_RES_TF():
+    import sympy as sp
+    ID = sp.eye(2)
+    phi, exp_1jOmegal_c, tArm, exp_2jOmegaL_c, K, r00, r10, r11, R, T, rITM, tSR, rho = sp.symbols('phi exp_1jOmegal_c tArm exp_2jOmegaL_c K r00 r10 r11 R T rITM tSR rho')
+    SEr = sp.Matrix([[sp.cos(phi), sp.sin(phi)], [-sp.sin(phi), sp.cos(phi)]])
+    SE = SEr * exp_1jOmegal_c
+    ARM = tArm * exp_2jOmegaL_c * sp.Matrix([[1, 0], [-K, 1]])
+    ARM_RES = sp.Matrix([[r00, 0], [r10, r11]])
+    rho_ARM = ARM_RES * ((R + T) * ARM - rITM * ID)
+    SEC = tSR * SE * rho_ARM * SE
+    SEC_RES = (ID + rho*SEC)**-1
+
+    subexprs, SEC_RES_expr = sp.cse(SEC_RES)
+    for expr in subexprs:
+        print(str(expr[0]), '=', str(expr[1]))
+    print('RES', '=', str(SEC_RES_expr[0]).replace('Matrix', 'np.array'))
+
+
 def shotradSignalRecycled(f, ifo):
+    """Quantum noise model for signal recycled IFO (see shotrad for more info)
+
+    New version July 2016 by JH based on transfer function formalism
+
+    coeff = frequency dependent overall noise coefficient (Nx1)
+            (not required anymore, but kept for compatibility with shotrad.m)
+    Mifo = IFO input-output relation for the AS port
+    Msig = signal transfer to the AS port
+    Mnoise = noise fields produced by losses in the IFO at the AS port
+
+    """
+    # f                                          # Signal Freq. [Hz]
+    lambda_ = ifo.Laser.Wavelength               # Laser Wavelength [m]
+    hbar    = scipy.constants.hbar               # Plancks Constant [Js]
+    c       = scipy.constants.c                  # SOL [m/s]
+    omega_0 = 2*pi*c/lambda_                     # Laser angular frequency [rads/s]
+
+    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+    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       = ifo.Materials.MirrorMass           # Mirror mass [kg]
+
+    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+    bsloss  = ifo.Optics.BSLoss                  # BS Loss [Power]
+    mismatch = 1 - ifo.Optics.coupling           # Mismatch
+    mismatch = mismatch + ifo.TCS.SRCloss        # Mismatch
+
+    # BSloss + mismatch has been incorporated into a SRC Loss
+    lambda_SR = 1 - (1 - mismatch) * (1 - bsloss)  # SR cavity loss [Power]
+
+    tau     = sqrt(ifo.Optics.SRM.Transmittance)  # SRM Transmittance [amplitude]
+    rho     = sqrt(1 - tau**2)                   # SRM Reflectivity [amplitude]
+
+    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+    ds      = ifo.Optics.SRM.Tunephase           # SRC Detunning
+    phi     = ds/2
+
+    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+    lambda_arm = 1 - (1 - ifo.Optics.Loss)**2 * (1 - ifo.Optics.ETM.Transmittance)
+
+    R = 1 - T - ifo.Optics.Loss                  # ITM Reflectivity [Power]
+
+    P = ifo.gwinc.parm                           # use precomputed value
+
+    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+    nf = len(f)
+
+    ID = np.array([[np.ones(nf), np.zeros(nf)], [np.zeros(nf), np.ones(nf)]])
+
+    # transfer matrices for dark port input and signal field
+    Mifo = zeros((2,2,nf), dtype=complex)
+    Msig = zeros((2,1,nf), dtype=complex)
+
+    # transfer matrices for SEC and arm loss fields
+    Mp = zeros((2,2,nf), dtype=complex)
+    Mn = zeros((2,2,nf), dtype=complex)
+
+    # SRC rotation matrix
+    SEr = np.array([[np.tile(cos(phi), nf), np.tile(sin(phi), nf)],
+                    [np.tile(-sin(phi), nf), np.tile(cos(phi), nf)]])
+
+    # some precomputed parameters
+    tITM = sqrt(T)
+    rITM = sqrt(R)
+
+    tArm = sqrt(1 - lambda_arm)
+    tSR = sqrt(1 - lambda_SR)
+    tSig = sqrt((1 - lambda_arm / 2) * (1 - lambda_SR / 2))
+    RT_SRM = rho**2 + tau**2
+
+    lossArm = sqrt(lambda_arm)
+    lossSR = sqrt(lambda_SR)
+
+    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)
+
+    # arm cavity
+    exp_2jOmegaL_c = exp(2j*Omega*L/c)
+    ARM = tArm * exp_2jOmegaL_c * np.array([[np.ones(nf), np.zeros(nf)], [-K, np.ones(nf)]])
+
+    # the following code is generated by compile_ARM_RES_TF()
+    # and is equivalent to the following:
+    # RES = zeros((2,2,nf), dtype=complex)
+    # RES = np.linalg.pinv(ID.transpose((2,0,1)) - rITM * ARM.transpose((2,0,1))).transpose((1,2,0))
+    x0 = exp_2jOmegaL_c*rITM*tArm
+    x1 = -x0 + 1
+    x2 = 1/x1
+    RES = np.array([[x2, np.zeros(nf)], [-K*x0/x1**2, x2]])
+    # end of generated code
+
+    rho_ARM = getProdTF(RES, (R + T) * ARM - rITM * ID)
+    tau_ARM = tITM * RES
+
+    # signal-extraction cavity
+    SE = SEr * exp(1j * Omega * l / c)
+    SEC = getProdTF(tSR * SE, getProdTF(rho_ARM, SE))
+
+    exp_1jOmegal_c = exp(1j*Omega*l/c)
+    r00 = RES[0,0,:]
+    r10 = RES[1,0,:]
+    r11 = RES[1,1,:]
+
+    # the following code is generated by compile_SEC_RES_TF()
+    # and is equivalent to the following:
+    # RES = zeros((2,2,nf), dtype=complex)
+    # RES = np.linalg.pinv(ID.transpose((2,0,1)) + rho * SEC.transpose((2,0,1))).transpose((1,2,0))
+    x0 = cos(phi)
+    x1 = exp_2jOmegaL_c*tArm*(R + T)
+    x2 = -rITM + x1
+    x3 = exp_1jOmegal_c**2*r11*tSR*x2
+    x4 = sin(phi)
+    x5 = exp_1jOmegal_c*r00*tSR*x2
+    x6 = exp_1jOmegal_c*tSR*(-K*r11*x1 + r10*x2)
+    x7 = exp_1jOmegal_c*(x0*x6 - x4*x5)
+    x8 = rho*(x0**2*x3 + x4*x7) + 1
+    x9 = x0*x3*x4
+    x10 = exp_1jOmegal_c*(x0*x5 + x4*x6)
+    x11 = x10*x4 + x9
+    x12 = x0*x7 - x9
+    x13 = rho*(x0*x10 - x3*x4**2) + 1
+    x14 = 1/(-rho**2*x11*x12 + x13*x8)
+    x15 = rho*x14
+    RES = np.array([[x14*x8, -x11*x15], [-x12*x15, x13*x14]])
+    # end of generated code
+
+    rho_SEC = getProdTF(RES, RT_SRM * SEC + rho * ID)
+    tau_SEC = tau * getProdTF(RES, SE)
+    tau_SEC_ARM = getProdTF(tau_SEC, tau_ARM)
+
+    # signal field
+    Msig = tSig * exp(1j * Omega * L / c) * \
+           getProdTF(tau_SEC_ARM, np.array([[np.zeros(nf)], [sqrt(2 * K) / h_SQL]]))
+
+    # dark-port input field
+    Mifo = rho_SEC
+
+    # loss field from arm cavity
+    Mn = lossArm * tau_SEC_ARM
+
+    # loss field from signal-extraction cavity
+    Mp = lossSR * tau_SEC
+
+    # adapt to GWINC phase convention
+    Msig = Msig[[1, 0], :, :]
+    Msig[1, 0, :] = -Msig[1, 0, :]
+
+    def adapt_to_gwinc(Mx):
+        My = zeros(Mx.shape, dtype=complex)
+        My[0, 0, :] = Mx[1, 1, :]
+        My[1, 1, :] = Mx[0, 0, :]
+        My[0, 1, :] = -Mx[1, 0, :]
+        My[1, 0, :] = -Mx[0, 1, :]
+        return My
+
+    Mifo = adapt_to_gwinc(Mifo)
+    Mn = adapt_to_gwinc(Mn)
+    Mp = adapt_to_gwinc(Mp)
+
+    # overall coefficient
+    coeff = 1
+
+    # put all loss fields together
+    Mnoise = np.hstack([Mn, Mp])
+
+    return coeff, Mifo, Msig, Mnoise
+
+
+def shotradSignalRecycledBnC(f, ifo):
     """Quantum noise model for signal recycled IFO
 
     See shotrad for more info.
@@ -342,7 +544,7 @@ def shotradSignalRecycled(f, ifo):
                          1/2*lambda_SR * ( sin(2*phi)+Kappa*sin(phi)**2) )
 
     C21_L   = tau**2 * ( (sin(2*phi)-Kappa*cos(phi)**2 ) +
-                         1/2*epsilon*cos(phi) * ( (3+exp_2jbeta )*Kappa*sin(phi) - 4*cos_beta**2*sin(phi) ) +
+                         1/2*epsilon*cos(phi) * ( (3+exp_2jbeta )*Kappa*cos(phi) - 4*cos_beta**2*sin(phi) ) +
                          1/2*lambda_SR * ( -sin(2*phi) + Kappa*cos(phi)**2) )
 
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -381,7 +583,7 @@ def shotradSignalRecycled(f, ifo):
                                     2*cos_beta*(invexp_1jbeta*cos(phi)-rho*exp_1jbeta*(cos(phi)+Kappa*sin(phi))))
     N22     = -sqrt(2*epsilon)*tau*(-invexp_1jbeta+rho*exp_1jbeta)*cos_beta*cos(phi)
     N12     = -sqrt(2*epsilon)*tau*(invexp_1jbeta+rho*exp_1jbeta)*cos_beta*sin(phi);
-    N21     = sqrt(2*epsilon)*tau*(-Kappa*(1+rho)*cos(phi)+
+    N21     = sqrt(epsilon/2)*tau*(-Kappa*(1+rho)*cos(phi)+
                                    2*cos_beta*(invexp_1jbeta+rho*exp_1jbeta)*cos_beta*sin(phi))
 
     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -