From 7534900a22097a2871d64410ffc8b5eaafcf44ad Mon Sep 17 00:00:00 2001
From: Christopher Wipf <wipf@ligo.mit.edu>
Date: Wed, 22 Aug 2018 18:47:56 -0700
Subject: [PATCH] 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 | 206 ++++++++++++++++++++++++++++++++++++++++-
1 file changed, 204 insertions(+), 2 deletions(-)
diff --git a/gwinc/noise/quantum.py b/gwinc/noise/quantum.py
index 6c0257e9..73f0112d 100644
--- a/gwinc/noise/quantum.py
+++ b/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))
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
GitLab