Select Git revision
doxygen_index.h
Forked from
einsteinathome / graphicsframework
Source project has a limited visibility.
ifo.py 43.21 KiB
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
import numpy as np
import math
import copy
import warnings
import cmath
from pykat import finesse
import pykat.components
import pykat.exceptions as pkex
import pykat.external.peakdetect as peak
import matplotlib.pyplot as plt
import pkg_resources
from scipy.optimize import fmin
global nsilica, clight
nsilica = 1.44963098985906
clight = 299792458.0
class aLIGO(object):
"""
Object storing the Finesse input file and some auxilliary information
about Avanced LIGO interferometers.
References:
[1] C. Bond `How to stay in shape: Overcoming beam and mirror distortions
in advanced gravitational wave interferometers', PhD thesis, University
of Birmingham, 2014, http://etheses.bham.ac.uk/5223/2/Bond14PhD.pdf
[2] D. Martynov, `Lock Acquisition and Sensitivity Analysis of Advanced
LIGO Interferometers', PhD thesis, Caltech, 2015
http://thesis.library.caltech.edu/8899/1/DenisMartynovThesis.pdf
[3] A. Staley, `Locking the Advanced LIGO Gravitational Wave Detector:
with a focus on the Arm Length Stabilization Technique', PhD thesis,
Columbia University, 2015
https://academiccommons.columbia.edu/catalog/ac:189457
"""
def __init__(self, _name="default", katfile=None, verbose = False, debug=False):
names = ['default', 'LLO', 'LHO']
if debug:
self.kat = finesse.kat(tempdir=".",tempname="test")
else:
self.kat = finesse.kat()
self.kat.verbose=verbose
self._data_path=pkg_resources.resource_filename('pykat.gw_detectors','finesse_files/')
if katfile:
self.kat.loadKatFile(katfile)
else:
if _name not in names: # TODO different files not yet implemented
printf("aLIGO name `{}' not recognised, must be 'default', 'LLO' or 'LHO'",_name)
#print(data_path)
self.kat.loadKatFile(self._data_path+"aLIGO.kat")
# ----------------------------------------------------------------------
# set variables to zero first
self.DCoffset = 0.0
self.DCoffsetW = 0.0
# ----------------------------------------------------------------------
# get and derive parameters from the kat file
# get main frequencies
if "f1" in self.kat.constants.keys():
self.f1 = float(self.kat.constants["f1"].value)
else:
self.f1 = 9099471.0
if "f2" in self.kat.constants.keys():
self.f2 = float(self.kat.constants["f2"].value)
else:
self.f2 = 5.0 * self.f1
if "f3" in self.kat.constants.keys():
self.f3 = float(self.kat.constants["f3"].value)
# TODO add else here!
# check modultion frequencies
if (5 * self.f1 != self.f2):
print(" ** Warning: modulation frequencies do not match: 5*f1!=f2")
# defining a dicotionary for the main mirror positions (tunings)
self.tunings = {}
self.tunings = self.get_tunings(self.kat)
self.compute_derived_lengths(self.kat)
# ----------------------------------------------------------------------
# define ports and signals
# useful ports
self.POP_f1 = port("POP_f1", "nPOP", self.f1, phase=101)
self.POP_f2 = port("POP_f2", "nPOP", self.f2, phase=13)
self.REFL_f1 = port("REFL_f1", "nREFL", self.f1, phase=101)
self.REFL_f2 = port("REFL_f2", "nREFL", self.f2, phase=14)
self.AS_DC = port("AS_DC", "nSRM2")
self.POW_BS = port("PowBS", "nPRBS*")
self.POW_X = port("PowX", "nITMX2")
self.POW_Y = port("PowY", "nITMY2")
# pretune DOF
self.preARMX = DOF("ARMX", self.POW_X, "", "ETMX", 1, 1.0)
self.preARMY = DOF("ARMY", self.POW_Y, "", "ETMY", 1, 1.0)
self.preMICH = DOF("AS" , self.AS_DC, "", ["ITMX", "ETMX", "ITMY", "ETMY"], [1,1,-1,-1], 6.0)
self.prePRCL = DOF("PRCL", self.POW_BS, "", "PRM", 1, 10.0)
self.preSRCL = DOF("SRCL", self.AS_DC, "", "SRM", 1, 10.0)
# control scheme as in [1] Table C.1
self.PRCL = DOF("PRCL", self.POP_f1, "I", "PRM", 1, 100.0)
self.MICH = DOF("MICH", self.POP_f2, "Q", ["ITMX", "ETMX", "ITMY", "ETMY"], [1,1,-1,-1], 100.0)
self.CARM = DOF("CARM", self.REFL_f1, "I", ["ETMX", "ETMY"], [1, 1], 1.5)
self.DARM = DOF("DARM", self.AS_DC, "", ["ETMX", "ETMY"], [1,-1], 1.0)
self.SRCL = DOF("SRCL", self.REFL_f2, "I", "SRM", 1, 1e2)
def adjust_PRC_length(self, kat, verbose=False):
"""
Adjust PRC length so that it fulfils the requirement
lPRC = (N+1/2) * c/(2*f1), see [1] equation C.1
In the current design N=3
"""
vprint(verbose, "-- adjusting PRC length")
ltmp = 0.5 * clight / self.f1
delta_l = 3.5 * ltmp - self.lPRC
vprint(verbose, " adusting kat.lp1.L by {:.4g}m".format(delta_l))
kat.lp1.L += delta_l
self.compute_derived_lengths(kat)
def check_f1_PRC_resonance(self, _kat):
"""
Plot the sideband amplitudes for mdulation frequecy
f1 (~ 9MHz) in the PRC, to check the resonance
condition.
"""
kat = _kat.deepcopy()
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
startf = self.f1-4000.0
stopf = self.f1+4000.0
if _kat.maxtem == "off":
nmStr = None
else:
nmStr = "0 0"
code = """
ad f1p {0} {1} nPRM2
ad f1m {0} -{1} nPRM2
xaxis mod1 f lin {2} {3} 200
put f1p f $x1
put f1m f $mx1
""".format(nmStr, self.f1, startf, stopf)
kat.parseCommands(code)
out = kat.run()
ax.plot(out.x-self.f1,np.abs(out["f1p"]), label=" f1")
ax.plot(out.x-self.f1,np.abs(out["f1m"]), label="-f1")
ax.set_xlim([np.min(out.x-self.f1), np.max(out.x-self.f1)])
ax.set_xlabel("delta_f1 [Hz]")
ax.set_ylabel('sqrt(W) ')
ax.grid()
ax.legend()
plt.tight_layout()
plt.show(block=0)
def compute_derived_lengths(self, kat, verbose=False):
"""
Compute derived length from individual space components.
Design values are currently:
lPRC = 57.656, lSRC = 56.008, lSchnupp = 0.08
and the individual lengths:
PRC: Lp1 16.6107, Lp2 16.1647, Lp3 19.5381
SRC: Ls1 15.7586, Ls2 15.4435, Ls3 19.3661
"""
# distances between HR surfaces:
self.lpr = kat.lp1.L + kat.lp2.L + kat.lp3.L
self.lx = kat.lx1.L + kat.BSsub1.L * kat.BSsub1.n + kat.ITMXsub.L * kat.ITMXsub.n
self.ly = kat.ly1.L + kat.ITMYsub.L * kat.ITMYsub.n
self.lsr = kat.ls1.L + kat.ls2.L + kat.ls3.L + kat.BSsub2.L * kat.BSsub2.n
# resulting combined distances (single, not roundtrip)
self.lMI = 0.5 * (self.lx + self.ly)
self.lPRC = self.lpr + self.lMI
self.lSRC = self.lsr + self.lMI
self.lSchnupp = self.lx - self.ly
self.compute_derived_resonances(kat)
def compute_derived_resonances(self, kat):
self.fsrX = 0.5 * clight / float(kat.LX.L)
self.fsrY = 0.5 * clight / float(kat.LY.L)
self.fsrPRC = 0.5 * clight / self.lPRC
self.fsrSRC = 0.5 * clight / self.lSRC
self.f1_PRC = 3.5 * self.fsrPRC
def lengths_status(self, kat):
print(" .--------------------------------------------------.")
print("| - arm length: |")
print("| Lx = {:11.7}m, Ly = {:11.7}m |".format(float(kat.LX.L), float(kat.LY.L)))
print("| - small MI and recycling lengths: | ")
print("| lx = {:11.7}m, ly = {:11.7}m |".format(self.lx, self.ly))
print("| lpr = {:11.7}m, lsr = {:11.7}m |".format(self.lpr, self.lsr))
print("| lMI = {:11.7}m, lSchnupp = {:11.5}m |".format(self.lMI, self.lSchnupp))
print("| lPRC = {:11.7}m, lSRC = {:11.7}m |".format(self.lPRC, self.lSRC))
print("+---------------------------------------------------+")
print("| - associated cavity frequencies [Hz]: |")
print("| fsrx = {:11.8}, fsry = {:11.8} |".format(self.fsrX, self.fsrY))
print("| fsrPRC = {:11.8}, fsrSRC = {:11.8} |".format(self.fsrPRC, self.fsrSRC))
print("| f1_PRC = {:11.8} |".format(self.f1_PRC))
print("| f1 = {:11.8}, f2 = {:11.9} |".format(self.f1, self.f2))
print(" `--------------------------------------------------'")
def get_tunings(self, kat):
self.tunings["maxtem"] = kat.maxtem
self.tunings["PRM"] = kat.PRM.phi
self.tunings["ITMX"] = kat.ITMX.phi
self.tunings["ETMX"] = kat.ETMX.phi
self.tunings["ITMY"] = kat.ITMY.phi
self.tunings["ETMY"] = kat.ETMY.phi
self.tunings["BS"] = kat.BS.phi
self.tunings["SRM"] = kat.SRM.phi
return self.tunings
def set_tunings(self, kat, tunings):
kat.maxtem = tunings["maxtem"]
kat.PRM.phi = tunings["PRM"]
kat.ITMX.phi = tunings["ITMX"]
kat.ETMX.phi = tunings["ETMX"]
kat.ITMY.phi = tunings["ITMY"]
kat.ETMY.phi = tunings["ETMY"]
kat.BS.phi = tunings["BS"]
kat.SRM.phi = tunings["SRM"]
def apply_lock_feedback(self, kat, out):
tuning = self.get_tunings(kat)
if "ETMX_lock" in out.ylabels:
tuning["ETMX"] += float(out["ETMX_lock"])
else:
print(" ** Warning: could not find ETMX lock")
if "ETMY_lock" in out.ylabels:
tuning["ETMY"] += float(out["ETMY_lock"])
else:
print(" ** Warning: could not find ETMY lock")
if "PRCL_lock" in out.ylabels:
tuning["PRM"] += float(out["PRCL_lock"])
else:
print(" ** Warning: could not find PRCL lock")
if ("MICH_lock" in out.ylabels) and ("ITMY_lock" in out.ylabels):
tuning["ITMX"] += float(out["MICH_lock"])
tuning["ITMY"] += float(out["ITMY_lock"])
else:
print(" ** Warning: could not find MICH (ITMY) lock")
if "SRCL_lock" in out.ylabels:
tuning["SRM"] += float(out["SRCL_lock"])
else:
print(" ** Warning: could not find SRCL lock")
self.set_tunings(kat, tuning)
def pretune(self, _kat, pretune_precision=1.0e-4, verbose=False):
print("-- pretuning interferometer to precision {0:2g} deg = {1:2g} m".format(pretune_precision, pretune_precision*_kat.lambda0/360.0))
kat=_kat.deepcopy()
#_kat.ITMX.phi=0.0
#_kat.ITMY.phi=0.0
vprint(verbose, " scanning X arm (maximising power)")
kat1 = kat.deepcopy()
make_transparent(kat1,["PRM","SRM"])
make_transparent(kat1,["ITMY", "ETMY"])
kat1.BS.setRTL(0.0,1.0,0.0) # set BS refl. for X arm
phi, precision = self.scan_to_precision(kat1, self.preARMX, pretune_precision)
phi=round(phi/pretune_precision)*pretune_precision
phi=round_to_n(phi,5)
vprint(verbose, " found max/min at: {} (precision = {:2g})".format(phi, precision))
self.preARMX.apply_tuning(_kat,phi)
vprint(verbose, " scanning Y arm (maximising power)")
kat = _kat.deepcopy()
make_transparent(kat,["PRM","SRM"])
make_transparent(kat,["ITMX", "ETMX"])
kat.BS.setRTL(1.0,0.0,0.0) # set BS refl. for Y arm
phi, precision = self.scan_to_precision(kat, self.preARMY, pretune_precision)
phi=round(phi/pretune_precision)*pretune_precision
phi=round_to_n(phi,5)
vprint(verbose, " found max/min at: {} (precision = {:2g})".format(phi, precision))
self.preARMY.apply_tuning(_kat,phi)
vprint(verbose, " scanning MICH (minimising power)")
kat = _kat.deepcopy()
make_transparent(kat,["PRM","SRM"])
phi, precision = self.scan_to_precision(kat, self.preMICH, pretune_precision, minmax="min", precision=30.0)
phi=round(phi/pretune_precision)*pretune_precision
phi=round_to_n(phi,5)
vprint(verbose, " found max/min at: {} (precision = {:2g})".format(phi, precision))
self.preMICH.apply_tuning(_kat,phi, add=True)
vprint(verbose, " scanning PRCL (maximising power)")
kat = _kat.deepcopy()
make_transparent(kat,["SRM"])
phi, precision = self.scan_to_precision(kat, self.prePRCL, pretune_precision)
phi=round(phi/pretune_precision)*pretune_precision
phi=round_to_n(phi,5)
vprint(verbose, " found max/min at: {} (precision = {:2g})".format(phi, precision))
self.prePRCL.apply_tuning(_kat,phi)
vprint(verbose, " scanning SRCL (maximising carrier power, then adding 90 deg)")
kat = _kat.deepcopy()
phi, precision = self.scan_to_precision(kat, self.preSRCL, pretune_precision, phi=0)
phi=round(phi/pretune_precision)*pretune_precision
phi=round_to_n(phi,4)-90.0
vprint(verbose, " found max/min at: {} (precision = {:2g})".format(phi, precision))
self.preSRCL.apply_tuning(_kat,phi)
print(" ... done")
def scan_to_precision(self, kat, DOF, pretune_precision, minmax="max", phi=0.0, precision=60.0):
while precision>pretune_precision*DOF.scale:
out = scan_DOF(kat, DOF, xlimits = [phi-1.5*precision, phi+1.5*precision])
phi, precision = find_peak(out, DOF.port.portName, minmax=minmax)
#print("** phi= {}".format(phi))
return phi, precision
def pretune_status(self, _kat):
kat = _kat.deepcopy()
kat.verbose = False
kat.noxaxis = True
pretune_DOFs = [self.preARMX, self.preARMY, self.prePRCL, self.preMICH, self.preSRCL]
_detStr=""
for p in pretune_DOFs:
_sigStr = p.port.signal(kat)
_detStr = "\n".join([_detStr, _sigStr])
kat.parseCommands(_detStr)
out = kat.run()
Pin = float(kat.L0.P)
tunings = self.get_tunings(kat)
_maxtemStr = "{:3}".format(tunings["maxtem"])
if tunings["maxtem"] == -1:
_maxtemStr="off"
print(" .--------------------------------------------------.")
print(" | pretuned for maxtem = {}, phase = {:2} |".format(_maxtemStr, int(kat.phase)))
keys_t = list(tunings.keys())
keys_t.remove("maxtem")
print(" .--------------------------------------------------.")
print(" | port power[W] pow. ratio | optics tunings |")
print(" +----------------------------|---------------------+")
idx_p = 0
idx_t = 0
while (idx_p < len(pretune_DOFs) or idx_t < len(keys_t)):
if idx_p < len(pretune_DOFs):
p = pretune_DOFs[idx_p]
print(" | {:5}: {:9.4g} {:9.4g} |".format(p.name, float(out[p.port.name]), float(out[p.port.name])/Pin),end="")
idx_p +=1
else:
print(" | |", end="")
if idx_t < len(keys_t):
t=keys_t[idx_t]
print(" {:5}: {:9.3g} |".format(t, float(self.tunings[t])))
idx_t +=1
else:
print(" |")
print(" `--------------------------------------------------'")
# probably extra and can be removed
def power_ratios(self, _kat):
kat = _kat.deepcopy()
kat.verbose = False
kat.noxaxis = True
ports = [self.POW_X, self.POW_Y, self.AS_DC, self.POW_BS]
_detStr = ""
for p in ports:
_sigStr = p.signal(kat)
_detStr = "\n".join([_detStr, _sigStr])
kat.parseCommands(_detStr)
out = kat.run()
Pin = float(kat.L0.P)
print("-- power ratios (Pin = {0:.3g} W)".format(Pin))
for p in ports:
print(" {0:6} = {1:8.3g} W ({0:6}/Pin = {2:8.2g})" .format(p.name, float(out[p.name]), float(out[p.name])/Pin))
return
def plot_pretuning_powers(self, _kat, xlimits=[-10,10]):
kat = _kat.deepcopy()
kat.verbose = False
kat.noxaxis = True
dofs = [self.preARMX, self.preARMY, self.preMICH, self.prePRCL]
idx=1
fig = plt.figure()
for d in dofs:
ax = fig.add_subplot(2,2,idx)
idx+=1
out = scan_DOF(kat, d, xlimits = np.multiply(d.scale, xlimits), relative = True)
ax.semilogy(out.x,out[d.signal_name(kat)])
ax.set_xlim([np.min(out.x), np.max(out.x)])
ax.set_xlabel("phi [deg] {}".format(d.optics[0]))
ax.set_ylabel('{} [W] '.format(d.signal_name(kat)))
ax.grid()
plt.tight_layout()
plt.show(block=0)
def plot_error_signals(self, _kat, xlimits=[-1,1]):
kat = _kat.deepcopy()
kat.verbose = False
kat.noxaxis = True
dofs = [self.DARM, self.CARM, self.PRCL, self.SRCL, self.MICH]
idx=1
fig = plt.figure()
for d in dofs:
ax = fig.add_subplot(2,3,idx)
idx+=1
out = scan_DOF(kat, d, xlimits = np.multiply(d.scale,xlimits), relative = True)
ax.plot(out.x,out[d.signal_name(kat)])
ax.set_xlim([np.min(out.x), np.max(out.x)])
ax.set_xlabel("{} [deg]".format(d.name))
ax.set_ylabel('{} [W] '.format(d.port.name))
ax.grid()
plt.tight_layout()
plt.show(block=0)
def set_DC_offset(self, _kat, DCoffset=None, verbose=False):
if DCoffset:
self.DCoffset=DCoffset
print("-- applying user-defined DC offset:")
pretuning = self.get_tunings(_kat)
pretuning["ETMY"] += self.DCoffset
pretuning["ETMX"] -= self.DCoffset
self.set_tunings(_kat, pretuning)
kat = _kat.deepcopy()
sigStr = self.AS_DC.signal(kat)
signame = self.AS_DC.signal_name(kat)
kat.parseCommands(sigStr)
kat.noxaxis=True
out = kat.run()
self.DCoffsetW=float(out[signame])
else:
# Finding light power in AS port (mostly due to RF sidebands now
kat = _kat.deepcopy()
sigStr = self.AS_DC.signal(kat)
signame = self.AS_DC.signal_name(kat)
kat.parseCommands(sigStr)
kat.noxaxis=True
out = kat.run()
print("-- adjusting DCoffset based on light in dark port:")
waste_light = round(float(out[signame]),1)
print(" waste light in AS port of {:2} W".format(waste_light))
kat_lock = _kat.deepcopy()
self.find_DC_offset(kat_lock, 2*waste_light)
pretuning = self.get_tunings(kat_lock)
pretuning["ETMY"] += self.DC_offset
pretuning["ETMX"] -= self.DC_offset
self.set_tunings(_kat, pretuning)
self.DCoffset_meter = self.DCoffset / 360.0 * _kat.lambda0
vprint(verbose, " DCoffset = {:6.4} deg ({:6.4}m)".format(self.DCoffset, self.DCoffset_meter))
vprint(verbose, " at dark port power: {:6.4}W".format(self.DCoffsetW))
def find_DC_offset(self, _kat, AS_power, precision=1e-4, verbose=False):
"""
Returns the DC offset of DARM that corrponds to the
specified power in the AS power.
"""
vprint(verbose, " finding DC offset for AS power of {:3g} W".format(AS_power))
kat = _kat.deepcopy()
kat.verbose = False
kat.noxaxis = True
_sigStr = self.AS_DC.signal(kat)
kat.parseCommands(_sigStr)
Xphi = float(kat.ETMX.phi)
Yphi = float(kat.ETMY.phi)
def powerDiff(phi, kat, Xphi, Yphi, AS_power):
kat.ETMY.phi = Yphi + phi
kat.ETMX.phi = Xphi - phi
out=kat.run()
#print(out[self.AS_DC.name]-AS_power)
return np.abs(out[self.AS_DC.name]-AS_power)
vprint(verbose, " starting peak search...")
out=fmin(powerDiff,0,xtol=precision,ftol=1e-3,args=(kat, Xphi, Yphi, AS_power), disp=verbose)
vprint(verbose, " ... done")
vprint(verbose, " DC offset for AS_DC={} W is: {}".format(AS_power, out[0]))
self.DCoffset = round(out[0],6)
self.DCoffsetW = AS_power
return self.DCoffset
def generate_tuning_block(self, kat):
code1 = """
%%% FTblock tunings
###########################################################################"""
code2 = """
const phi_ITMX {:.8}
const phi_ITMY {:.8}
const phi_ETMX {:.8}
const phi_ETMY {:.8}""".format(float(kat.ITMX.phi), float(kat.ITMY.phi), float(kat.ETMX.phi),float(kat.ETMY.phi))
code3 = """
const phi_BS {:.8}
const phi_PRM {:.8}
const phi_SRM {:.8}""".format(float(kat.BS.phi), float(kat.PRM.phi), float(kat.SRM.phi))
code4 = """
###########################################################################
%%% FTend tuning"""
return "".join([code1, code2, code3, code4])
def generate_errsig_block(self, kat, noplot=False):
code1 = """
%%% FTblock errsigs
###########################################################################
"""
sigDARM = self.DARM.signal(kat)
sigCARM = self.CARM.signal(kat)
sigPRCL = self.PRCL.signal(kat)
sigMICH = self.MICH.signal(kat)
sigSRCL = self.SRCL.signal(kat)
code2 = "\n".join([sigDARM, sigCARM, sigPRCL, sigMICH, sigSRCL])
code3= ""
if noplot:
nameDARM = self.DARM.signal_name(kat)
nameCARM = self.CARM.signal_name(kat)
namePRCL = self.PRCL.signal_name(kat)
nameMICH = self.MICH.signal_name(kat)
nameSRCL = self.SRCL.signal_name(kat)
code3 = """
noplot {}
noplot {}
noplot {}
noplot {}
noplot {}""".format(nameDARM, nameCARM, namePRCL, nameMICH, nameSRCL)
code4 = """
###########################################################################
%%% FTend errsigs"""
return "".join([code1, code2, code3, code4])
def generate_locks(self, _kat, gainsAdjustment = [0.5, 0.005, 1.0, 0.5, 0.025], gains=None, accuracies=None, verbose=True):
"""
gainsAdjustment: factors to apply to loop gains computed from optical gains
gains: override loop gain [W per deg]
accuracies: overwrite error signal threshold [W]
"""
kat = _kat.deepcopy()
# optical gains in W/rad
ogDARM = optical_gain(kat, self.DARM, self.DARM)
ogCARM = optical_gain(kat, self.CARM, self.CARM)
ogPRCL = optical_gain(kat, self.PRCL, self.PRCL)
ogMICH = optical_gain(kat, self.MICH, self.MICH)
ogSRCL = optical_gain(kat, self.SRCL, self.SRCL)
if gains == None:
# manually tuning relative gains
factor = -1.0 * 180 / math.pi # convert from rad/W to -1 * deg/W
gainDARM = round_to_n(gainsAdjustment[0] * factor / ogDARM, 2) # manually tuned
gainCARM = round_to_n(gainsAdjustment[1] * factor / ogCARM, 2) # factor 0.005 for better gain hirarchy with DARM
gainPRCL = round_to_n(gainsAdjustment[2] * factor / ogPRCL, 2) # manually tuned
gainMICH = round_to_n(gainsAdjustment[3] * factor / ogMICH, 2) # manually tuned
gainSRCL = round_to_n(gainsAdjustment[4] * factor / ogSRCL, 2) # gain hirarchy with MICH
gains = [ gainDARM, gainCARM, gainPRCL, gainMICH, gainSRCL]
self.lockGains = gains.copy()
# rms: loop accuracies in meters (manually tuned for the loops to work
# with the default file)
# to compute accuracies from rms, we convert
# rms to radians as rms_rad = rms * 2 pi/lambda
# and then multiply by the optical gain.
rms = [1e-13, 1e-13, 1e-12, 1e-11, 50e-11] # default accuracies in meters
factor = 2.0 * math.pi / kat.lambda0 # convert from m to radians
if accuracies == None:
accDARM = round_to_n(np.abs(factor * rms[0] * ogDARM),2)
accCARM = round_to_n(np.abs(factor * rms[1] * ogCARM),2)
accPRCL = round_to_n(np.abs(factor * rms[2] * ogPRCL),2)
accMICH = round_to_n(np.abs(factor * rms[3] * ogMICH),2)
accSRCL = round_to_n(np.abs(factor * rms[4] * ogSRCL),2)
accuracies = [accDARM, accCARM, accPRCL, accMICH, accSRCL]
self.lockAccuracies = accuracies.copy()
nameDARM = self.DARM.signal_name(kat)
nameCARM = self.CARM.signal_name(kat)
namePRCL = self.PRCL.signal_name(kat)
nameMICH = self.MICH.signal_name(kat)
nameSRCL = self.SRCL.signal_name(kat)
self.lockNames = [nameDARM, nameCARM, namePRCL, nameMICH, nameSRCL]
factor1 = 2.0 * math.pi / 360.0
factor2 = 2.0 * math.pi / kat.lambda0
factor3 = 360.0 / kat.lambda0
factor4 = -1.0 * 180 / math.pi
if verbose:
print(" .--------------------------------------------------.")
print(" | Parameters for locks: |")
print(" +--------------------------------------------------+")
print(" | -- optical gains [W/rad], [W/deg] and [W/m]: |")
print(" | DARM: {:12.5}, {:12.5}, {:12.5} |".format(ogDARM, ogDARM*factor1, ogDARM*factor2))
print(" | CARM: {:12.5}, {:12.5}, {:12.5} |".format(ogCARM, ogCARM*factor1, ogCARM*factor2))
print(" | PRCL: {:12.5}, {:12.5}, {:12.5} |".format(ogPRCL, ogPRCL*factor1, ogPRCL*factor2))
print(" | MICH: {:12.5}, {:12.5}, {:12.5} |".format(ogMICH, ogMICH*factor1, ogMICH*factor2))
print(" | SRCL: {:12.5}, {:12.5}, {:12.5} |".format(ogSRCL, ogSRCL*factor1, ogSRCL*factor2))
print(" +--------------------------------------------------+")
print(" | -- defult loop accuracies [deg], [m] and [W]: |")
print(" | DARM: {:12.6}, {:12.6}, {:12.6} |".format(factor3*rms[0], rms[0], np.abs(rms[0]*ogDARM*factor2)))
print(" | CARM: {:12.6}, {:12.6}, {:12.6} |".format(factor3*rms[1], rms[1], np.abs(rms[1]*ogCARM*factor2)))
print(" | PRCL: {:12.6}, {:12.6}, {:12.6} |".format(factor3*rms[2], rms[2], np.abs(rms[2]*ogPRCL*factor2)))
print(" | MICH: {:12.6}, {:12.6}, {:12.6} |".format(factor3*rms[3], rms[3], np.abs(rms[3]*ogMICH*factor2)))
print(" | SRCL: {:12.6}, {:12.6}, {:12.6} |".format(factor3*rms[4], rms[4], np.abs(rms[4]*ogSRCL*factor2)))
print(" +--------------------------------------------------+")
print(" | -- extra gain factors (factor * 1/optical_gain): |")
print(" | DARM: {:5.4} * {:12.6} = {:12.6} |".format(gainsAdjustment[0],factor4/ogDARM, gainsAdjustment[0]*factor4/ogDARM))
print(" | CARM: {:5.4} * {:12.6} = {:12.6} |".format(gainsAdjustment[1],factor4/ogCARM, gainsAdjustment[1]*factor4/ogCARM))
print(" | PRCL: {:5.4} * {:12.6} = {:12.6} |".format(gainsAdjustment[2],factor4/ogPRCL, gainsAdjustment[2]*factor4/ogPRCL))
print(" | MICH: {:5.4} * {:12.6} = {:12.6} |".format(gainsAdjustment[3],factor4/ogMICH, gainsAdjustment[3]*factor4/ogMICH))
print(" | SRCL: {:5.4} * {:12.6} = {:12.6} |".format(gainsAdjustment[4],factor4/ogSRCL, gainsAdjustment[4]*factor4/ogSRCL))
print(" `--------------------------------------------------'")
def generate_lock_block(self, kat, verbose=False):
if self.lockNames == None or self.lockAccuracies == None or self.lockGains == None:
raise pkex.BasePyKatException("run gerate_locks before generate_lock_block")
code1 = """
%%% FTblock locks
###########################################################################
set AS_f2_I_re {} re
set CARM_err {} re
set PRCL_err {} re
set MICH_err {} re
set SRCL_err {} re
func DARM_err = $AS_f2_I_re - {}
""".format(self.lockNames[0], self.lockNames[1], self.lockNames[2], self.lockNames[3], self.lockNames[4], self.DCoffsetW)
code2 = """
lock DARM_lock $DARM_err {:8.2} {:8.2}
lock CARM_lock $CARM_err {:8.2g} {:8.2g}
lock PRCL_lock $PRCL_err {:8.2g} {:8.2g}
lock MICH_lock $MICH_err {:8.2g} {:8.2g}
lock SRCL_lock $SRCL_err {:8.2g} {:8.2g}""".format(self.lockGains[0], self.lockAccuracies[0], self.lockGains[1], self.lockAccuracies[1], self.lockGains[2], self.lockAccuracies[2], self.lockGains[3], self.lockAccuracies[3], self.lockGains[4], self.lockAccuracies[4])
code3 = """
noplot ITMY_lock
func ITMY_lock = (-1.0) * $MICH_lock
func ETMX_lock = $CARM_lock + $MICH_lock + $DARM_lock
func ETMY_lock = $CARM_lock - $MICH_lock - $DARM_lock
put* PRM phi $PRCL_lock
put* ITMX phi $MICH_lock
put* ITMY phi $ITMY_lock
put* ETMX phi $ETMX_lock
put* ETMY phi $ETMY_lock
put* SRM phi $SRCL_lock
noplot PRCL_lock
noplot SRCL_lock
noplot MICH_lock
noplot DARM_lock
noplot CARM_lock
noplot ETMX_lock
noplot ETMY_lock
###########################################################################
%%% FTend locks"""
if verbose:
print(" .--------------------------------------------------.")
print(" | Lock commands used: |")
print(" +--------------------------------------------------+")
for l in code2.splitlines():
print (" | {:49}|".format(l))
print(" `--------------------------------------------------'")
return "".join([code1, code2, code3])
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
class DOF(object):
"""
Defining a degree of freedom for the interferometer, includes the
objects and how to move them, and the default output port to read
out the DOF signal.
"""
def __init__(self, _DOFName, _port, _quad, _optics, _factors, _scale, sigtype="z"):
self.name = _DOFName
self.port = _port
self.quad = _quad
self.sigtype = sigtype
self.optics=make_list_copy(_optics)
self.factors=make_list_copy(_factors)
# scaling factor, to compensate for lower sensitivity compared
# to DARM (in tuning plots for example)
# Thus DARM has a scale of 1, all other DOFs a scale >1
self.scale = _scale
def apply_tuning(self, kat, phi, add=False):
for idx, o in enumerate(self.optics):
if add:
kat.components[o].phi += phi * self.factors[idx]
else:
kat.components[o].phi = phi * self.factors[idx]
def signal(self, kat):
return self.port.signal(kat, self.quad, sigtype=self.sigtype)
def signal_name(self, kat):
return self.port.signal_name(kat, self.quad, sigtype=self.sigtype)
def transfer(self, kat, fsig, phase2=None):
return self.port.transfer(kat, self.quad, fsig=fsig, phase2=phase2, sigtype=self.sigtype)
def transfer_name(self, kat):
return self.port.transfer_name(kat, self.quad)
def fsig(self, _fsigName, fsig=1.0):
_fsigStr= ""
for idx, o in enumerate(self.optics):
phase = 0.0
if self.factors[idx] == -1:
phase = 180.0
_fsigStr = "\n".join([_fsigStr, "fsig {} {} {} {} ".format(_fsigName, o, fsig, phase)])
return _fsigStr
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
class port(object):
"""
Defining an output port for the interferometer, can be either a
pd or a pd1 detector (for error signal generation).
"""
def __init__(self, _portName, _nodeNames, f=None, phase=None):
self.portName = _portName
self.nodeNames = make_list_copy(_nodeNames)
self.f=f # demodulation frequency, float
self.phase = phase # demodulation frequency for I quadrature, float
self.name = self.portName
def check_nodeName(self, kat):
self.nodeName = None
for node in self.nodeNames:
_node = node
if _node[-1] == "*":
_node = _node[:-1]
if _node in kat.nodes:
self.nodeName=node
break
if self.nodeName==None:
raise pkex.BasePyKatException("port {}: cannot find any of these nodes: '{}'".format(self.name,self.nodeNames))
def amplitude_name(self, kat, f, n=None, m=None, sigtype="z"):
name = self.name + "_ad"
return name
def amplitude(self, kat, f, n=None, m=None, sigtype="z"):
self.check_nodeName(kat)
name = self.amplitude_name(kat, f, n=n, m=m, sigtype=sigtype)
if n==None and m==None:
return "ad {} {} {}".format(name, f, self.nodeName)
else:
return "ad {} {} {} {} {}".format(name, f, n, m, self.nodeName)
def signal_name(self, kat, quad="I", sigtype="z"):
name = self.name
if self.f!=None:
name = self.name+"_"+quad
return name
def signal(self, kat, quad="I", sigtype="z"):
self.check_nodeName(kat)
name = self.signal_name(kat, quad=quad, sigtype=sigtype)
if sigtype != "z":
raise pkex.BasePyKatException("alignment signals are not implemented yet")
if self.f==None:
return "pd {} {}".format(name, self.nodeName)
else:
if quad !="I" and quad != "Q":
raise pkex.BasePyKatException("quadrature must be 'I' or 'Q'")
phase = self.IQ_phase(quad, self.phase)
return "pd1 {} {} {} {}".format(name, self.f, phase, self.nodeName)
def IQ_phase(self, quad, phase):
if quad== "Q":
phase = phase + 90.0
if phase >=360.0 :
phase -= 360.0
return phase
def transfer_name(self, kat, quad="I"):
name = self.name
if self.f!=None:
name = self.name+"_"+quad
return name
def transfer(self, kat, quad="I", fsig=1.0, phase2=None, sigtype="z"):
self.check_nodeName(kat)
name = self.transfer_name(kat, quad=quad)
if sigtype!="z":
raise pkex.BasePyKatException("alignment signals are not implemented yet")
if self.f==None:
if phase2 == None:
return "pd1 {} {} {}".format(name, fsig, self.nodeName)
else:
return "pd1 {} {} {} {}".format(name, fsig, phase2, self.nodeName)
else:
if quad !="I" and quad != "Q":
raise pkex.BasePyKatException("quadrature must be 'I' or 'Q'")
phase = self.IQ_phase(quad, self.phase)
if phase2 == None:
return "pd2 {} {} {} {} {}".format(name, self.f , phase, fsig, self.nodeName)
else:
return "pd2 {} {} {} {} {} {}".format(name, self.f, phase, fsig, phase2, self.nodeName)
def scan_optics_string(_optics, _factors, _varName, linlog="lin", xlimits=[-100, 100], steps=200, axis=1,relative=False):
optics=make_list_copy(_optics)
factors=make_list_copy(_factors)
if len(optics) != len(factors):
raise pkex.BasePyKatException("you must provide a factor for each optics")
if linlog not in ["lin", "log"]:
raise pkex.BasePyKatException("linlog must be 'lin' or 'log'")
_tuneStr = "var {} 0\n".format(_varName)
if axis==1:
_tuneStr += "xaxis {} phi {} {} {} {}".format(_varName, linlog, xlimits[0], xlimits[1], steps)
elif (axis==2 or axis==3):
_tuneStr += "x{}axis {} phi {} {} {} {}".format(axis, _varName, linlog, xlimits[0], xlimits[1], steps)
else:
raise pkex.BasePyKatException("axis must be 1, 2 or 3")
_putStr = ""
for idx, o in enumerate(optics):
if factors[idx] == 1:
_xStr="$x"
elif factors[idx] == -1:
_xStr="$mx"
else:
raise pkex.BasePyKatException("optics factors must be 1 or -1")
if (relative):
_putCmd = "put*"
else:
_putCmd = "put"
_putStr = "\n".join([_putStr, "{} {} phi {}{}".format(_putCmd, o, _xStr, axis)])
_tuneStr += _putStr
return _tuneStr
def make_transparent(kat, _components):
"""
Function to make certain mirror or beamsplitter objects transparent
Usage:
make_transparent(kat, "PRM")
or
make_transparent(kat, ["PRM", "SRM"])
The function modifies the kat object. Use something like
kat = copy.deepcopy(basekat)
before if you want to preseve the original kat object
"""
components=make_list_copy(_components)
for o in kat.components.values():
if o.name in components:
if isinstance(o, pykat.components.AbstractMirrorComponent):
o.setRTL(0.0,1.0,0.0)
components = [c for c in components if c != o.name]
else:
raise pkex.BasePyKatException("{} is not a mirror or beamsplitter".format(o))
#print(' made {} transparent'.format(o))
if len(components) != 0:
raise pkex.BasePyKatException("Cannot find component {}".format(components))
return kat
def reconnect_nodes(kat, component1, idx1, node_name, verbose=False):
c_string = component1.getFinesseText()
c = c_string[0].split()
new_string = " ".join(c[:-2])
nodes = {}
nodes[0] = c[-2]
nodes[1] = c[-1]
nodes[idx1]=node_name
new_string = new_string + " " + nodes[0] + " " + nodes[1]
vprint(verbose, " new string ='{}'".format(new_string))
kat.parseCommands(new_string)
def remove_commands(kat, _commands, verbose=False):
commands=make_list_copy(_commands)
# removing commands
for o in kat.commands.values():
if o.name in commands:
o.remove()
commands = [c for c in commands if c != o.name]
vprint(verbose, ' {} removed'.format(o))
if len(commands) != 0:
raise pkex.BasePyKatException("Cannot find command(s) {}".format(commands))
return kat
def remove_components(kat, _components, component_in=None, component_out=None, verbose=False):
components=make_list_copy(_components)
if kat.components[components[-1]].nodes[1]:
node_in = kat.components[components[-1]].nodes[1].name
else:
node_in = None
node_out = kat.components[components[0]].nodes[0].name
# removing components
for o in kat.components.values():
if o.name in components:
o.remove()
components = [c for c in components if c != o.name]
vprint(verbose, ' {} removed'.format(o))
if len(components) != 0:
raise pkex.BasePyKatException("Cannot find component(s) {}".format(components))
# reconnecting nodes if requested
if component_in:
reconnect_nodes(kat, kat.components[component_in],1, node_in)
if component_out:
reconnect_nodes(kat, kat.components[component_out],0, node_out)
return kat
def BS_optical_path(thickness, n=nsilica, angle=45.0):
"""
Compute optical path length in BS substrate, default
parameters assume angle of incidence of 45 deg and fused
silica substrate.
thickness: substrate thickness [m]
n: index of refraction of substrate
angle: angle of incidence (in vacuum) [deg]
"""
angle_subst = math.asin(math.sin(math.radians(angle))/n)
L = thickness / math.cos(angle_subst)
return math.degrees(angle_subst), L
def scan_DOF(_kat, DOF, xlimits=[-100, 100], steps=200, relative=False):
kat = _kat.deepcopy()
scan_string = scan_optics_string(DOF.optics, DOF.factors, "scan", linlog="lin", xlimits=xlimits, steps=steps, axis=1,relative=relative)
kat.parseCommands(scan_string)
sigStr = DOF.signal(kat)
kat.parseCommands(sigStr)
out = kat.run()
return out
def scan_optics(_kat, _optics, _factors, xlimits=[-100, 100], steps=200,relative=False):
"""
Scans one or more optics (by changing its tuning).
Parameters:
optics: list of names of components to be tuned
factors: list of scaling factors for the tuning for each optics, first element must be 1.0
xlimits: limits of scan, defaults to [-100, 100]
steps: number of steps to use in scan, default is 200
Usage:
scan_optis(kat, "PRM", 1)
scan_optis(kat, ["ETMX", "ETMY", [1, -1])
"""
kat = _kat.deepcopy()
optics=make_list_copy(_optics)
factors=make_list_copy(_factors)
scan_string = scan_optics_string(optics, factors, "scan", linlog="lin", xlimits=xlimits, steps=steps, axis=1,relative=relative)
kat.parseCommands(scan_string)
out = kat.run()
return out
def optical_gain(_kat, DOF_sig, DOF_det, f=10.0):
kat = _kat.deepcopy()
_fsigStr = DOF_sig.fsig("sig1", fsig=f)
_detStr = DOF_det.transfer(kat, fsig=f)
_detName = DOF_det.transfer_name(kat)
kat.parseCommands(_fsigStr)
kat.parseCommands(_detStr)
kat.noxaxis = True
kat.parseCommands("yaxis lin abs:deg")
out = kat.run()
return float(np.real(out[_detName]))
def find_peak(out, detector, minmax='max', debug=False):
"""
Expects an output of a kat scan and find the max/min output on
a sepcific detector. Useful for pre-tuning cavities or interferometers.
Returns the tuning of the maximum (or minimum) and the precision
Parameters:
detector: name of detetor to use
minmax, string to indicate maximum or minum detection, default 'max'
Usage:
find_peak(out, "pdout")
find_peak(out, "pdout", minmax='min')
"""
stepsize = out.x[1]-out.x[0]
if debug:
print(" stepsize (precision) of scan: {0:g}".format(stepsize))
_max, _min = peak.peakdetect( out[detector],out.x, 1)
if debug:
plt.figure()
plt.plot(out.x,out[detector])
print("max: ")
print(_max)
print("min: ")
print(_min)
if minmax == 'max':
X = [p[0] for p in _max]
Y = [p[1] for p in _max]
X_out = X[np.argmax(Y)]
Y_out = np.max(Y)
elif minmax == 'min':
X = [p[0] for p in _min]
Y = [p[1] for p in _min]
X_out = X[np.argmin(Y)]
Y_out = np.min(Y)
else:
raise pkex.BasePyKatException("maxmin must be 'max' or 'min'")
if debug:
plt.plot(X_out,Y_out,'o')
plt.xlabel('tuning [deg]')
plt.ylabel('{0} output'.format(detector))
plt.show()
return X_out, stepsize
def make_list_copy(_l):
"""
Utility function, takes a list of strings or single string
and returns a copy of a list, e.g.
"string" copy to ["string"]
["string1", "string2"] copy to ["string1", "string2"]
"""
if not isinstance(_l, list):
_l = [_l]
return _l[:] # copy the list, just to be save
def round_to_n(x, n):
if not x: return 0
power = -int(math.floor(math.log10(abs(x)))) + (n - 1)
factor = (10 ** power)
return round(x * factor) / factor
def vprint(verbose, printstr):
if verbose:
print(printstr)