Skip to content
Snippets Groups Projects
Commit 0e90ed52 authored by Kevin Kuns's avatar Kevin Kuns
Browse files

add quadrature phase noise 

Fixed the old calculation and added a new term to the quantum sub-budget
parent fd8c9324
No related branches found
No related tags found
No related merge requests found
Pipeline #163333 failed
Stage: dist
Stage: test
Stage: review
Hide whitespace changes
Inline Side-by-side
gwinc/ifo/CE1/__init__.py
+ 2
− 1
View file @ 0e90ed52
......@@ -14,10 +14,11 @@ class QuantumVacuum(nb.Budget):
noises = [
QuantumVacuumAS,
QuantumVacuumArm,
QuantumVacuumSRC,
QuantumVacuumSEC,
QuantumVacuumFilterCavity,
QuantumVacuumInjection,
QuantumVacuumReadout,
QuantumVacuumQuadraturePhase,
]
......
gwinc/ifo/CE2/__init__.py
+ 2
− 1
View file @ 0e90ed52
......@@ -14,10 +14,11 @@ class QuantumVacuum(nb.Budget):
noises = [
QuantumVacuumAS,
QuantumVacuumArm,
QuantumVacuumSRC,
QuantumVacuumSEC,
QuantumVacuumFilterCavity,
QuantumVacuumInjection,
QuantumVacuumReadout,
QuantumVacuumQuadraturePhase,
]
......
gwinc/ifo/CE2/ifo.yaml
+ 1
− 0
View file @ 0e90ed52
......@@ -321,6 +321,7 @@ Squeezer:
AmplitudedB: 15 # SQZ amplitude [dB]
InjectionLoss: 0.02 # power loss to sqz
SQZAngle: 0 # SQZ phase [radians]
LOAngleRMS: 10e-3 # quadrature noise [radians]
# Parameters for frequency dependent squeezing
FilterCavity:
......
gwinc/ifo/noises.py
+ 27
− 8
View file @ 0e90ed52
......@@ -165,8 +165,8 @@ def precomp_quantum(f, ifo):
pbs, parm, _, _, _ = ifo_power(ifo)
power = Struct(parm=parm, pbs=pbs)
noise_dict = noise.quantum.shotrad(f, ifo, mirror_mass, power)
pc.ASvac = noise_dict['asvac']
pc.SRC = noise_dict['src']
pc.ASvac = noise_dict['ASvac']
pc.SEC = noise_dict['SEC']
pc.Arm = noise_dict['arm']
pc.Injection = noise_dict['injection']
pc.PD = noise_dict['pd']
......@@ -175,13 +175,19 @@ def precomp_quantum(f, ifo):
# are noises from the unsqueezed vacuum injected at the back mirror
# Right now there are at most one filter cavity in all the models;
# if there were two, there would also be FC1 and FC1_unsqzd_back, etc.
keys = list(noise_dict.keys())
fc_keys = [key for key in keys if 'FC' in key]
# keys = list(noise_dict.keys())
fc_keys = [key for key in noise_dict.keys() if 'FC' in key]
if fc_keys:
pc.FC = np.zeros_like(pc.ASvac)
for key in fc_keys:
pc.FC += noise_dict[key]
if 'phase' in noise_dict.keys():
pc.Phase = noise_dict['phase']
if 'ofc' in noise_dict.keys():
pc.OFC = noise_dict['OFC']
return pc
......@@ -273,18 +279,18 @@ class QuantumVacuumArm(nb.Noise):
return quantum.Arm
class QuantumVacuumSRC(nb.Noise):
"""Quantum vacuum due to SRC loss
class QuantumVacuumSEC(nb.Noise):
"""Quantum vacuum due to SEC loss
"""
style = dict(
label='SRC Loss',
label='SEC Loss',
color='xkcd:cerulean'
)
@nb.precomp(quantum=precomp_quantum)
def calc(self, quantum):
return quantum.SRC
return quantum.SEC
class QuantumVacuumFilterCavity(nb.Noise):
......@@ -329,6 +335,19 @@ class QuantumVacuumReadout(nb.Noise):
return quantum.PD
class QuantumVacuumQuadraturePhase(nb.Noise):
"""Quantum vacuum noise due to quadrature phase noise
"""
style = dict(
label='Quadrature Phase',
color='xkcd:slate'
)
@nb.precomp(quantum=precomp_quantum)
def calc(self, quantum):
return quantum.Phase
class StandardQuantumLimit(nb.Noise):
"""Standard Quantum Limit
......
gwinc/noise/quantum.py
+ 189
− 182
View file @ 0e90ed52
......@@ -4,6 +4,7 @@
from __future__ import division
import numpy as np
from numpy import pi, sqrt, arctan, sin, cos, exp, log10, conj
from copy import deepcopy
from .. import logger
from .. import const
......@@ -17,71 +18,71 @@ def sqzOptimalSqueezeAngle(Mifo, eta):
return alpha
def getSqzType(ifo):
def getSqzParams(ifo):
"""Determine squeezer type, if any, and extract common parameters
Returns a struct with the following attributes:
SQZ_DB: squeezing in dB
ANTISQZ_DB: anti-squeezing in dB
alpha: freq Indep Squeeze angle [rad]
lambda_in: loss to squeezing before injection [Power]
etaRMS: quadrature noise [rad]
eta: homodyne angle [rad]
"""
params = Struct()
if 'Squeezer' not in ifo:
sqzType = 'None'
sqzType = None
else:
sqzType = ifo.Squeezer.get('Type', 'Freq Independent')
return sqzType
params.sqzType = sqzType
def getSqzParams(ifo):
"""Get squeezer parameters
Returns:
SQZ_DB: squeezing [dB]
alpha: squeezing angle [rad]
lambda_in: loss to squeezing before injection [Power]
ANTISQZ_DB: anti-squeezing [dB]
etaRMS: quadrature phase noise [rad]
"""
# determine squeezer type, if any
# and extract common parameters
sqzType = getSqzType(ifo)
# extract common parameters
if sqzType == 'None':
SQZ_DB = 0 # Squeeing in dB
alpha = 0 # Squeeze angle
lambda_in = 0 # Loss to squeezing before injection [Power]
ANTISQZ_DB = 0 # Anti squeezing in db
etaRMS = 0
else:
SQZ_DB = ifo.Squeezer.AmplitudedB # Squeeing in dB
lambda_in = ifo.Squeezer.InjectionLoss # Loss to squeezing before injection [Power]
alpha = ifo.Squeezer.SQZAngle # Freq Indep Squeeze angle
ANTISQZ_DB = ifo.Squeezer.get('AntiAmplitudedB', SQZ_DB) # Anti squeezing in db
etaRMS = ifo.Squeezer.get('LOAngleRMS', 0) # quadrature noise
# extract squeezer parameters
if sqzType is None:
params.SQZ_DB = 0
params.ANTISQZ_DB = 0
params.alpha = 0
params.lambda_in = 0
params.etaRMS = 0
# switch on squeezing type for other input squeezing modifications
if sqzType == 'None':
pass
else:
params.SQZ_DB = ifo.Squeezer.AmplitudedB
params.ANTISQZ_DB = ifo.Squeezer.get('AntiAmplitudedB', params.SQZ_DB)
params.alpha = ifo.Squeezer.SQZAngle
params.lambda_in = ifo.Squeezer.InjectionLoss
params.etaRMS = ifo.Squeezer.get('LOAngleRMS', 0)
elif sqzType == 'Freq Independent':
logger.debug('You are injecting %g dB of frequency independent squeezing' % SQZ_DB)
# Homodyne Readout phase
eta_orig = ifo.Optics.get('Quadrature', Struct()).get('dc', None)
elif sqzType == 'Optimal':
# compute optimal squeezing angle
alpha = sqzOptimalSqueezeAngle(Mifo, eta)
ifoRead = ifo.get('Squeezer', Struct()).get('Readout', None)
if ifoRead is None:
eta = eta_orig
if eta_orig is None:
raise Exception("must add Quadrature.dc or Readout...")
logger.debug('You are injecting %g dB of squeezing with optimal frequency dependent squeezing angle' % SQZ_DB)
elif ifoRead.Type == 'DC':
eta = np.sign(ifoRead.fringe_side) * \
np.arccos((ifoRead.defect_PWR_W / ifoRead.readout_PWR_W)**.5)
elif sqzType == 'OptimalOptimal':
# compute optimal squeezing angle, assuming optimal readout phase
R = SQZ_DB / (20 * log10(exp(1)))
MnPD = sqzInjectionLoss(Mn, lambda_PD)
MsigPD = Msig * sqrt(1 - lambda_PD)
alpha = sqzOptimalSqueezeAngle(Mifo, [], [R, lambda_in], MsigPD, MnPD)
elif ifoRead.Type == 'Homodyne':
eta = ifoRead.Angle
logger.debug('You are injecting %g dB of squeezing with optimal FD squeezing angle, for optimal readout phase' % SQZ_DB)
else:
raise Exception("Unknown Readout Type")
elif sqzType == 'Freq Dependent':
logger.debug('You are injecting %g dB of squeezing with frequency dependent squeezing angle' % SQZ_DB)
if eta_orig is not None:
# logger.warn((
# 'Quadrature.dc is redundant with '
# 'Squeezer.Readout and is deprecated.'
# ))
if eta_orig != eta:
raise Exception("Quadrature.dc inconsistent with Readout eta")
else:
raise Exception('ifo.Squeezer.Type must be None, Freq Independent, Optimal, or Frequency Dependent, not "%s"' % sqzType)
params.eta = eta
return SQZ_DB, alpha, lambda_in, ANTISQZ_DB, etaRMS
return params
######################################################################
......@@ -90,10 +91,11 @@ def getSqzParams(ifo):
def shotrad(f, ifo, mirror_mass, power):
"""Quantum noise noise strain spectrum
"""Quantum noise strain spectrum
:f: frequency array in Hz
:ifo: gwinc IFO Struct
:mirror_mass: mirror mass in kg
:power: gwinc power Struct
:returns: displacement noise power spectrum at :f:
......@@ -105,10 +107,11 @@ def shotrad(f, ifo, mirror_mass, power):
# Vacuum sources will be individually tracked with the Mnoise dict #
# The keys are the following #
# * arm: arm cavity losses #
# * src: SRC losses #
# * asvac: AS port vacuum #
# * SEC: SEC losses #
# * ASvac: AS port vacuum #
# * pd: readout losses #
# * injection: injection losses #
# * phase: quadrature phase noise #
# #
# If there is a filter cavity there are an additional set of keys #
# starting with 'FC' for the filter cavity losses #
......@@ -116,153 +119,64 @@ def shotrad(f, ifo, mirror_mass, power):
# If there is an output filter cavity the key 'OFC' gives its losses #
######################################################################
#####################################################
# Call IFO Quantum Model
#####################################################
Nfreq = len(f)
# call the main IFO Quantum Model
if 'Type' not in ifo.Optics:
fname = shotradSignalRecycled
else:
namespace = globals()
fname = namespace['shotrad' + ifo.Optics.Type]
coeff, Mifo, Msig, Mn = fname(f, ifo, mirror_mass, power)
Mnoise = dict(arm=Mn[:, :2, :], src=Mn[:, 2:, :])
#####################################################
# Readout
#####################################################
# separate arm and SEC loss
Mnoise = dict(arm=Mn[:, :2, :], SEC=Mn[:, 2:, :])
# useful numbers
# TODO, adjust Struct to allow deep defaulted access
# Homodyne Readout phase
eta_orig = ifo.Optics.get('Quadrature', Struct()).get('dc', None)
sqz_params = getSqzParams(ifo)
ifoRead = ifo.get('Squeezer', Struct()).get('Readout', None)
if ifoRead is None:
eta = eta_orig
if eta_orig is None:
raise Exception("must add Quadrature.dc or Readout...")
elif ifoRead.Type == 'DC':
eta = np.sign(ifoRead.fringe_side) * np.arccos((ifoRead.defect_PWR_W / ifoRead.readout_PWR_W)**.5)
elif ifoRead.Type == 'Homodyne':
eta = ifoRead.Angle
else:
raise Exception("Unknown Readout Type")
if sqz_params.sqzType == 'Optimal':
# compute optimal squeezing angle
sqz_params.alpha = sqzOptimalSqueezeAngle(Mifo, sqz_params.eta)
if eta_orig is not None:
# logger.warn((
# 'Quadrature.dc is redundant with '
# 'Squeezer.Readout and is deprecated.'
# ))
if eta_orig != eta:
raise Exception("Quadrature.dc inconsistent with Readout eta")
vHD = np.array([[sin(sqz_params.eta), cos(sqz_params.eta)]])
def homodyne(signal):
"""Readout the eta quadrature of the signal signal
"""
return np.squeeze(np.sum(abs(getProdTF(vHD, signal))**2, axis=1))
# optomechanical plant
lambda_PD = 1 - ifo.Optics.PhotoDetectorEfficiency # PD losses
Msig = Msig * sqrt(1 - lambda_PD)
plant = homodyne(Msig)
#####################################################
# Input Squeezing
#####################################################
#>>>>>>>> QUANTUM NOISE POWER SPECTRAL DENSITY WITH SQZ [BnC PRD 2004, 62]
#<<<<<<<<<<<<<<<<< Modified to include losses (KM)
#<<<<<<<<<<<<<<<<< Modified to include frequency dependent squeezing angle (LB)
SQZ_DB, alpha, lambda_in, ANTISQZ_DB, etaRMS = getSqzParams(ifo)
sqzType = getSqzType(ifo)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# ------------------------------------------- equation 63 BnC PRD 2004
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Define input matrix of Squeezing
R = SQZ_DB / (20 * log10(exp(1))) # Squeeze factor
R_anti = ANTISQZ_DB / (20 * log10(exp(1))) # Squeeze factor
Msqz = np.array([[exp(-R), 0], [0, exp(R_anti)]])
# expand to Nfreq
Msqz = np.transpose(np.tile(Msqz, (Nfreq,1,1)), axes=(1,2,0))
# add input rotation
MsqzRot = make2x2TF(cos(alpha), -sin(alpha), sin(alpha), cos(alpha))
Msqz = getProdTF(MsqzRot, Msqz)
Msqz = dict(asvac=Msqz)
# Include losses (lambda_in=ifo.Squeezer.InjectionLoss)
Msqz = sqzInjectionLoss(Msqz, lambda_in, 'injection')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Inject squeezed field into the IFO via some filter cavities
if sqzType == 'Freq Dependent' and 'FilterCavity' in ifo.Squeezer:
logger.debug(' Applying %d input filter cavities' % np.atleast_1d(ifo.Squeezer.FilterCavity).size)
Mr, Msqz = sqzFilterCavityChain(
f, np.atleast_1d(ifo.Squeezer.FilterCavity), Msqz)
#####################################################
# IFO Transfer and Output Filter Cavities
#####################################################
# apply the IFO dependent squeezing matrix to get the total noise
# due to quantum fluctuations coming in from the AS port
Msqz = {key: getProdTF(Mifo, nn) for key, nn in Msqz.items()}
# add this to the other noises
Mnoise.update(Msqz)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# pass IFO output through some filter cavities
if 'OutputFilter' in ifo:
if ifo.OutputFilter.Type == 'None':
# do nothing, say nothing
pass
elif ifo.OutputFilter.Type == 'Chain':
logger.debug(' Applying %d output filter cavities' % np.atleast_1d(ifo.OutputFilter.FilterCavity).size)
Mr, Mnoise = sqzFilterCavityChain(
f, np.atleast_1d(ifo.OutputFilter.FilterCavity), Mnoise, key='OFC')
Msig = getProdTF(Mr, Msig)
# Mnoise = getProdTF(Mn, Mnoise);
elif ifo.OutputFilter.Type == 'Optimal':
raise NotImplementedError("Cannot do optimal phase yet")
logger.debug(' Optimal output filtering!')
# get all the losses in the squeezed quadrature
Mnoise = propagate_noise_fc_ifo(
f, ifo, Mifo, Mnoise, sqz_params, sqz_params.alpha)
# compute optimal angle, including upcoming PD losses
MnPD = sqzInjectionLoss(Mnoise, lambda_PD)
zeta = sqzOptimalReadoutPhase(Msig, MnPD)
# add quadrature phase fluctuations if any
if sqz_params.etaRMS:
# get all losses in the anti-squeezed quadrature
Mnoise_antisqz = propagate_noise_fc_ifo(
f, ifo, Mifo, Mnoise, sqz_params, sqz_params.alpha + pi/2)
# rotate by that angle, less the homodyne angle
#zeta_opt = eta;
cs = cos(zeta - eta)
sn = sin(zeta - eta)
Mrot = make2x2TF(cs, -sn, sn, cs)
Mnoise = getProdTF(Mrot, Mnoise)
Msig = getProdTF(Mrot, Msig)
else:
raise Exception('ifo.OutputFilter.Type must be None, Chain or Optimal, not "%s"' % ifo.OutputFilter.Type)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# add PD efficiency
Mnoise = sqzInjectionLoss(Mnoise, lambda_PD, 'pd')
Msig = Msig * sqrt(1 - lambda_PD)
# sum over these since they're not individually tracked
var_antisqz = np.zeros_like(f)
for nn in Mnoise_antisqz.values():
var_antisqz += homodyne(nn)
# and compute the final noise
def HDnoise(eta, Mn):
vHD = np.array([[sin(eta), cos(eta)]])
n = coeff * np.squeeze(np.sum(abs(getProdTF(vHD, Mn))**2, axis=1)) / \
np.squeeze(np.sum(abs(getProdTF(vHD, Msig))**2, axis=1))
return n
# The total variance is
# V(alpha) * cos(etaRMS)^2 + V(alpha + pi/2) * sin(etaRMS)^2
Mnoise = {key: cos(sqz_params.etaRMS)**2 * homodyne(nn)
for key, nn in Mnoise.items()}
Mnoise['phase'] = sin(sqz_params.etaRMS)**2 * var_antisqz
# include quadrature noise (average of +- the RMS)
n = {key: (HDnoise(eta+etaRMS, Mn) + HDnoise(eta-etaRMS, Mn)) / 2
for key, Mn in Mnoise.items()}
else:
Mnoise = {key: homodyne(nn) for key, nn in Mnoise.items()}
n = {key: nn * ifo.Infrastructure.Length**2 for key, nn in n.items()}
# calibrate into displacement
# coeff can be removed if shotradSignalRecycledBnC isn't used
Mnoise = {key: coeff*nn/plant for key, nn in Mnoise.items()}
psd = {key: nn * ifo.Infrastructure.Length**2 for key, nn in Mnoise.items()}
return n
return psd
######################################################################
......@@ -625,6 +539,99 @@ def shotradSignalRecycledBnC(f, ifo, mirror_mass, power):
return coeff, Mifo, Msig, Mnoise
def propagate_noise_fc_ifo(f, ifo, Mifo, Mnoise, sqz_params, alpha):
"""Propagate quantum noise through the filter cavities and IFO
f: frequency vector [Hz]
ifo: ifo Struct
Mifo: IFO input-output relation for the AS port
Mnoise: dictionary of existing noise sources
sqz_params: Struct of squeeze parameters
alpha: squeeze quadrature
Returns a new Mnoise dict with the noises due to AS port vacuum, injection
loss, and any input or output filter cavities added.
"""
Mnoise = deepcopy(Mnoise)
#>>>>>>>> QUANTUM NOISE POWER SPECTRAL DENSITY WITH SQZ [BnC PRD 2004, 62]
#<<<<<<<<<<<<<<<<< Modified to include losses (KM)
#<<<<<<<<<<<<<<<<< Modified to include frequency dependent squeezing angle (LB)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# ------------------------------------------- equation 63 BnC PRD 2004
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Define input matrix of Squeezing
R = sqz_params.SQZ_DB / (20 * log10(exp(1))) # squeeze factor
R_anti = sqz_params.ANTISQZ_DB / (20 * log10(exp(1))) # anti-squeeze factor
Msqz = np.array([[exp(-R), 0], [0, exp(R_anti)]])
# expand to Nfreq
Msqz = np.transpose(np.tile(Msqz, (len(f), 1, 1)), axes=(1, 2, 0))
# add input rotation
MsqzRot = make2x2TF(cos(alpha), -sin(alpha), sin(alpha), cos(alpha))
Msqz = getProdTF(MsqzRot, Msqz)
Msqz = dict(ASvac=Msqz)
# Include losses (lambda_in=ifo.Squeezer.InjectionLoss)
Msqz = sqzInjectionLoss(Msqz, sqz_params.lambda_in, 'injection')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Inject squeezed field into the IFO via some filter cavities
if sqz_params.sqzType == 'Freq Dependent' and 'FilterCavity' in ifo.Squeezer:
logger.debug(' Applying %d input filter cavities' % np.atleast_1d(ifo.Squeezer.FilterCavity).size)
Mr, Msqz = sqzFilterCavityChain(
f, np.atleast_1d(ifo.Squeezer.FilterCavity), Msqz)
# apply the IFO dependent squeezing matrix to get the total noise
# due to quantum fluctuations coming in from the AS port
Msqz = {key: getProdTF(Mifo, nn) for key, nn in Msqz.items()}
# add this to the other noises
Mnoise.update(Msqz)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# pass IFO output through some filter cavities
if 'OutputFilter' in ifo:
if ifo.OutputFilter.Type == 'None':
# do nothing, say nothing
pass
elif ifo.OutputFilter.Type == 'Chain':
logger.debug(' Applying %d output filter cavities' % np.atleast_1d(ifo.OutputFilter.FilterCavity).size)
Mr, Mnoise = sqzFilterCavityChain(
f, np.atleast_1d(ifo.OutputFilter.FilterCavity), Mnoise, key='OFC')
Msig = getProdTF(Mr, Msig)
# Mnoise = getProdTF(Mn, Mnoise);
elif ifo.OutputFilter.Type == 'Optimal':
raise NotImplementedError("Cannot do optimal phase yet")
logger.debug(' Optimal output filtering!')
# compute optimal angle, including upcoming PD losses
MnPD = sqzInjectionLoss(Mnoise, lambda_PD)
zeta = sqzOptimalReadoutPhase(Msig, MnPD)
# rotate by that angle, less the homodyne angle
#zeta_opt = eta;
cs = cos(zeta - eta)
sn = sin(zeta - eta)
Mrot = make2x2TF(cs, -sn, sn, cs)
Mnoise = getProdTF(Mrot, Mnoise)
Msig = getProdTF(Mrot, Msig)
else:
raise Exception('ifo.OutputFilter.Type must be None, Chain or Optimal, not "%s"' % ifo.OutputFilter.Type)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# add PD efficiency
lambda_PD = 1 - ifo.Optics.PhotoDetectorEfficiency # PD losses
Mnoise = sqzInjectionLoss(Mnoise, lambda_PD, 'pd')
return Mnoise
######################################################################
# Auxiliary functions
######################################################################
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment