Skip to content
Snippets Groups Projects

Update quantum code from matgwinc, and vectorize

Merged Update quantum code from matgwinc, and vectorize
fast-quantum-matmul into master
Merged Christopher Wipf requested to merge fast-quantum-matmul into master
Viewing commit 7534900a
Show latest version
1 file
+ 204
− 2
Compare changes
  • Side-by-side
  • Inline
  • 7534900a 
    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
+ 204
− 2
@@ -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))
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -