diff --git a/gwinc/ifo/CE1/__init__.py b/gwinc/ifo/CE1/__init__.py
index 7baa4055c4d34c6b16ab5b709cad340707552dad..c558f3f551ab27bc92b6c5efe0a5a5b4c6bfac7f 100644
--- a/gwinc/ifo/CE1/__init__.py
+++ b/gwinc/ifo/CE1/__init__.py
@@ -14,10 +14,11 @@ class QuantumVacuum(nb.Budget):
noises = [
QuantumVacuumAS,
QuantumVacuumArm,
- QuantumVacuumSRC,
+ QuantumVacuumSEC,
QuantumVacuumFilterCavity,
QuantumVacuumInjection,
QuantumVacuumReadout,
+ QuantumVacuumQuadraturePhase,
]
diff --git a/gwinc/ifo/CE2/__init__.py b/gwinc/ifo/CE2/__init__.py
index f4e0e69d0008982b6b01d346141f2ac4a10b9fc7..dd608d7f8e27b8dec0f6ee1ffb2ba2d3b9e3434b 100644
--- a/gwinc/ifo/CE2/__init__.py
+++ b/gwinc/ifo/CE2/__init__.py
@@ -14,10 +14,11 @@ class QuantumVacuum(nb.Budget):
noises = [
QuantumVacuumAS,
QuantumVacuumArm,
- QuantumVacuumSRC,
+ QuantumVacuumSEC,
QuantumVacuumFilterCavity,
QuantumVacuumInjection,
QuantumVacuumReadout,
+ QuantumVacuumQuadraturePhase,
]
diff --git a/gwinc/ifo/CE2/ifo.yaml b/gwinc/ifo/CE2/ifo.yaml
index 97f5f0d70a4c2ab7b8ce3c230acd42fcff23f373..40caefe3a6c872a3b54d1be6ec14ded0e74bba42 100644
--- a/gwinc/ifo/CE2/ifo.yaml
+++ b/gwinc/ifo/CE2/ifo.yaml
@@ -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:
diff --git a/gwinc/ifo/noises.py b/gwinc/ifo/noises.py
index c3826f4bbfdeb92cafd289cf4ec86f7d558305c4..7e5b05b9990d3eb2d0d4b04519fd33d6d45f54b5 100644
--- a/gwinc/ifo/noises.py
+++ b/gwinc/ifo/noises.py
@@ -174,8 +174,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']
@@ -184,13 +184,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
@@ -280,18 +286,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):
@@ -336,6 +342,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
diff --git a/gwinc/noise/quantum.py b/gwinc/noise/quantum.py
index 68fe6421455f99335bb0353b9e9f4505f0678f52..260a1ae428f0ce58e3aba5476f4ec44b772a9b1c 100644
--- a/gwinc/noise/quantum.py
+++ b/gwinc/noise/quantum.py
@@ -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
######################################################################