Commit 6816af56 authored by Michele Valentini's avatar Michele Valentini
Browse files

Merge branch 'master' of git.ligo.org:finesse/pykat

parents 003845bb 3b286d5c
v0.0, 23/07/2013 -- Initial release.
\ No newline at end of file
PyKat-1.1.297
Adding in A+ model and files
pykat.ifo plot_strain_sensitivity now returns the kat object that it runs to get the plot
Plotting is now done using Matplotlib styles now
Pykat can be run without having Finesse installed now
\ No newline at end of file
This diff is collapsed.
......@@ -22,6 +22,7 @@ import six
########################
# Global helper functions
isContainer = lambda c: (not isinstance(c, six.string_types)) and hasattr(c, "__iter__")
is_container = lambda c: (not isinstance(c, six.string_types)) and hasattr(c, "__iter__")
########################
import imp
......@@ -49,17 +50,19 @@ from pykat.style import use as set_plot_style
from .SIfloat import SIfloat
msg = "Could not find the finesse executable 'kat'" \
"or you do not have the permissions to run it."
class nokat(object):
def __getattribute__(self, attr):
warn(msg)
try:
kat = finesse.kat()
v = kat.finesse_version()
except pkex.MissingFinesse:
from warnings import warn
msg = "Could not find the finesse executable 'kat'" \
"or you do not have the permissions to run it."
warn(msg)
class nokat(object):
def __getattribute__(self, attr):
warn(msg)
kat = nokat()
v = str(__min_req_finesse__)
if float(v.split('-')[0]) < __min_req_finesse__:
......@@ -67,6 +70,16 @@ if float(v.split('-')[0]) < __min_req_finesse__:
str(__min_req_finesse__),
v))
def info():
print("Pykat version: " + __version__)
print("Pykat loaded from: " + __file__)
if kat != nokat():
print("Finesse version: " + str(v))
print("Finesse loaded from: " + str(kat._finesse_exec()))
else:
print("Finesse could not be initialised")
SI = {'yocto': 1E-24, # yocto
'zepto': 1E-21, # zepto
'atto': 1E-18, # atto
......
......@@ -343,7 +343,8 @@ class cavity(Command):
base: The kat object to do the tracing with
"""
kat = self._kat.deepcopy()
kat.maxtem = 0 # always need a maxtem
for cav in kat.getAll(pykat.commands.cavity):
if self.name != cav.name:
cav.remove()
......@@ -357,7 +358,7 @@ class cavity(Command):
kat.parse("pd p %s" % self.__n1)
_, T = kat.run(getTraceData=True)
qx, qy, _ = T[0][node]
qx, qy, _ = T[0][str(node)]
return qx, qy
......
......@@ -479,17 +479,6 @@ class KatRun(object):
if yaxis is not None:
kat.yaxis = yaxis
if "log" in kat.yaxis:
if kat.xaxis.scale == "log":
plot_cmd = pyplot.loglog
else:
plot_cmd = pyplot.semilogy
else:
if kat.xaxis.scale == "log":
plot_cmd = pyplot.semilogx
else:
plot_cmd = pyplot.plot
dual_plot = False
_func1 = np.abs
_func2 = None
......@@ -497,20 +486,43 @@ class KatRun(object):
plot_cmd1 = None
plot_cmd2 = None
if "re:im" in kat.yaxis or "abs:deg" in kat.yaxis or "db:deg" in kat.yaxis:
dual_plot = True
if dual_plot:
fig = plt.figure(width="full", height=1)
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
else:
fig = plt.figure(width="full")
ax1 = fig.add_subplot(1, 1, 1)
ax2 = ax1
if "log" in kat.yaxis:
if kat.xaxis.scale == "log":
plot_cmd1 = ax1.loglog
plot_cmd2 = ax2.loglog
else:
plot_cmd1 = ax1.semilogy
plot_cmd2 = ax2.semilogy
else:
if kat.xaxis.scale == "log":
plot_cmd1 = ax1.semilogx
plot_cmd2 = ax2.semilogx
else:
plot_cmd1 = ax1.plot
plot_cmd2 = ax2.plot
if "re:im" in kat.yaxis:
_func1 = np.real
_func2 = np.imag
plot_cmd1 = plot_cmd2 = plot_cmd
dual_plot = True
elif "abs:deg" in kat.yaxis:
_func1 = np.abs
_func2 = lambda x: np.rad2deg(np.angle(x))
plot_cmd1 = plot_cmd
plot_cmd2 = pyplot.plot if kat.xaxis.scale == "lin" else pyplot.semilogx
plot_cmd2 = ax2.plot if kat.xaxis.scale == "lin" else ax2.semilogx
dual_plot = True
elif "db:deg" in kat.yaxis:
if "db" not in original_yaxis:
_func1 = lambda x: 10*np.log10(x)
......@@ -519,33 +531,24 @@ class KatRun(object):
_func2 = lambda x: np.rad2deg(np.angle(x))
plot_cmd1 = plot_cmd
plot_cmd2 = pyplot.plot if kat.xaxis.scale == "lin" else pyplot.semilogx
plot_cmd2 = ax2.plot if kat.xaxis.scale == "lin" else ax2.semilogx
dual_plot = True
elif "abs" in kat.yaxis:
# _func1 = np.abs
_func1 = np.real
plot_cmd1 = plot_cmd
elif "db" in kat.yaxis:
if "db" not in original_yaxis:
_func1 = lambda x: 10*np.log10(x)
else:
_func1 = lambda x: x
plot_cmd1 = plot_cmd
elif "deg" in kat.yaxis:
_func1 = lambda x: np.rad2deg(np.angle(x))
plot_cmd1 = plot_cmd
if dual_plot:
fig = plt.figure(width="full", height=1)
else:
fig = plt.figure(width="full")
if detectors is None:
detectors = [lbl.split()[0] for lbl in self.ylabels]
detectors = list(set(detectors))
detectors.sort()
......@@ -553,17 +556,12 @@ class KatRun(object):
for det in detectors:
if not hasattr(kat, det) or (hasattr(kat, det) and not getattr(kat, det).noplot):
if dual_plot:
ax = pyplot.subplot(2,1,1)
if styles is not None and det in styles:
l, = plot_cmd1(_x, _func1(self[det])/pykat.SI[y1scale], styles[det], label=det)
else:
l, = plot_cmd1(_x, _func1(self[det])/pykat.SI[y1scale], label=det)
if dual_plot:
pyplot.subplot(2,1,2)
plot_cmd2(_x, _func2(self[det])/pykat.SI[y2scale], color=l.get_color(), ls=l.get_linestyle(), label=det)
if dual_plot:
......@@ -598,33 +596,33 @@ class KatRun(object):
font_label_size = pyplot.rcParams["font.size"]-1
if dual_plot:
ax = pyplot.subplot(2,1,1)
pyplot.xlabel(xlabel, fontsize=font_label_size)
pyplot.ylabel(ylabel, fontsize=font_label_size)
pyplot.xlim(xlim[0], xlim[1])
#ax = plt.subplot(2,1,1)
ax1.set_xlabel(xlabel, fontsize=font_label_size)
ax1.set_ylabel(ylabel, fontsize=font_label_size)
ax1.set_xlim(xlim[0], xlim[1])
if ylim is not None:
pyplot.ylim(ylim[0],ylim[1])
ax1.set_ylim(ylim[0],ylim[1])
if title is not None:
pyplot.title(title, fontsize=font_label_size)
ax1.set_title(title, fontsize=font_label_size)
pyplot.subplot(2,1,2)
pyplot.xlabel(x2label, fontsize=font_label_size)
pyplot.ylabel(y2label, fontsize=font_label_size)
#plt.subplot(2,1,2)
ax2.set_xlabel(x2label, fontsize=font_label_size)
ax2.set_ylabel(y2label, fontsize=font_label_size)
pyplot.xlim(x2lim[0], x2lim[1])
ax2.set_xlim(x2lim[0], x2lim[1])
if y2lim is not None:
pyplot.ylim(y2lim[0],y2lim[1])
ax2.set_ylim(y2lim[0],y2lim[1])
else:
pyplot.xlabel(xlabel, fontsize=font_label_size)
pyplot.ylabel(ylabel)
pyplot.xlim(xlim[0], xlim[1])
ax1.set_xlabel(xlabel, fontsize=font_label_size)
ax1.set_ylabel(ylabel)
ax1.set_xlim(xlim[0], xlim[1])
if title is not None:
pyplot.title(title, fontsize=font_label_size)
ax1.set_title(title, fontsize=font_label_size)
if ylim is not None:
pyplot.ylim(ylim[0],ylim[1])
ax1.set_ylim(ylim[0],ylim[1])
pyplot.margins(0, 0.05)
pyplot.tight_layout()
......
......@@ -1291,7 +1291,7 @@ class IFO(object):
import pykat
from pykat.ifo import aligo
base = aligo.make_kat("design_with_IMC_HAM2_FI_OMC")
base = aligo.make_kat("design")
base.maxtem = 2
base = aligo.setup(base)
......@@ -1486,21 +1486,65 @@ class Output(object):
You can add many detectors at a given port, which readout different quadratures or types of transfer functions
or signals.
"""
def __init__(self, IFO, name, nodeNames, f_property_name=None, phase=0, block=None):
def __init__(self, IFO, name, nodeNames, f_property_name=None, phase=0, block=None,
alternate_node_name=None):
"""
An Output object is a generic description of an output port of an interferometer.
At an output port, such as the AS port, we can do DC or RF demodulations. Each
output can then compute the transfer functions or or photodiode error signals,
or optical amplitudes.
Parameters
----------
IFO : pykat.ifo.IFO
Interferometer class to associate with
name : str
Name of the output
nodeNames : list(str)
List of nodes that this output uses, one for normal detectors, multiple for
balanced homodyne like outputs.
f_property_name : str
Name of the frequency attribute in IFO. This should be a name of a frequency
attribute in the IFO class, this is used when computing what frequency to
demodulate at. None means a DC output
phase : foat
Demodulation phase
block : str
Which block in the kat object to add the commands to
"""
if f_property_name is not None and not isinstance(f_property_name, six.string_types):
raise Exception("f_property_name should be a property name of a frequency in class {}".format(IFO.__class__))
self.__IFO = IFO
self.name = name
self.nodeNames = make_list_copy(nodeNames)
self.__nodeNames = make_list_copy(nodeNames)
if alternate_node_name is None:
self.__alternate_node_name = None
else:
self.__alternate_node_name = make_list_copy(alternate_node_name)
self.__f_property_name = f_property_name
self.phase = phase # demodulation phase for I quadrature, float
self._block = block
if pykat.is_container(nodeNames):
self.num_nodes = len(nodeNames)
self.__nodeNames = [self.__nodeNames]
else:
self.num_nodes = 1
self.__nodeNames = [self.__nodeNames]
if self.num_nodes > 1 and self.__f_property_name is not None:
raise Exception("Can't use multiple node outputs and RF calculations")
@property
def kat(self):
"""For referencing the kat object this DOF is associated with"""
"""For referencing the actual kat object this output is associated with"""
return self.__IFO.kat
@property
def f(self):
......@@ -1509,20 +1553,23 @@ class Output(object):
return getattr(self.__IFO, self.__f_property_name)
def check_nodeName(self):
self.nodeName = None
for node in self.nodeNames:
_node = node
if _node[-1] == "*":
_node = _node[:-1]
if _node in self.__IFO.kat.nodes:
self.nodeName=node
if self.__alternate_node_name is not None:
_nodes = zip(self.__nodeNames, self.__alternate_node_name)
else:
_nodes = self.__nodeNames
for nodes in _nodes:
if all([_.strip('* ') in self.__IFO.kat.nodes for _ in nodes]):
self.nodeName = nodes
break
if self.nodeName==None:
raise pkex.BasePyKatException("port {}: cannot find any of these nodes: '{}'".format(self.name,self.nodeNames))
if self.nodeName is None:
raise pkex.BasePyKatException("port {}: cannot find any of these nodes: '{}'".format(self.name, self.__nodeNames))
def get_amplitude_name(self, f, n=None, m=None):
''''
......@@ -1547,7 +1594,11 @@ class Output(object):
rtn = "{}_{}_{}{}_ad".format(self.name,fstr,n,m)
return rtn
def get_amplitude_cmds(self, f, n=None, m=None):
if self.num_nodes != 1:
raise Exception("Amplitude output doesn't work for multiple node outputs")
rtn = []
self.check_nodeName()
......@@ -1555,11 +1606,12 @@ class Output(object):
name = self.get_amplitude_name(f, n=n, m=m)
if n==None and m==None:
rtn.append("ad {} {} {}".format(name, f, self.nodeName))
rtn.append("ad {} {} {}".format(name, f, self.nodeName[0]))
else:
rtn.append("ad {} {} {} {} {}".format(name, f, n, m, self.nodeName))
rtn.append("ad {} {} {} {} {}".format(name, f, n, m, self.nodeName[0]))
return rtn
def _pdtype(self, name, sigtype):
rtn = []
......@@ -1571,6 +1623,7 @@ class Output(object):
return rtn
def add_signal(self, quad=None, sigtype=None):
"""
Adds a photodiode detector to the kat object at this port. Must
......@@ -1595,6 +1648,7 @@ class Output(object):
self.__IFO.kat.parse(cmds, addToBlock=self._block)
return self.get_signal_name(quad, sigtype)
def get_signal_name(self, quad="I", sigtype='z'):
name = self.name
......@@ -1609,6 +1663,7 @@ class Output(object):
return name
def get_signal_cmds(self, dof=None, **kwargs):
"""
Returns the Finesse commands for a detector added to this port's location.
......@@ -1654,19 +1709,20 @@ class Output(object):
name = self.get_signal_name(quad=quad, sigtype=sigtype)
if self.f==None:
rtn.append("pd {} {}".format(name, self.nodeName))
rtn.append("pd {} {}".format(name, self.nodeName[0]))
else:
if quad !="I" and quad != "Q":
raise pkex.BasePyKatException("quadrature must be 'I' or 'Q'")
phase = self.IQ_phase(quad, self.phase)
rtn.append("pd1 {} {} {} {}".format(name, self.f, phase, self.nodeName))
rtn.append("pd1 {} {} {} {}".format(name, self.f, phase, self.nodeName[0]))
rtn.extend(self._pdtype(name, sigtype))
return rtn
def IQ_phase(self, quad, phase):
if phase is None:
raise pkex.BasePyKatException("Phase cannot be None")
......@@ -1679,11 +1735,13 @@ class Output(object):
return phase
def add_transfer(self, quad=None, sigtype="z"):
cmds = self.get_transfer_cmds(quad=quad, sigtype=sigtype)
self.__IFO.kat.parse(cmds, addToBlock=self._block)
return self.get_transfer_name(quad, sigtype)
def get_transfer_name(self, quad="I", sigtype="z"):
name = self.name
......@@ -1697,32 +1755,49 @@ class Output(object):
return name + "_TF"
def get_transfer_cmds(self, quad="I", sigtype="z", phase2=None):
self.check_nodeName()
name = self.get_transfer_name(quad=quad, sigtype=sigtype)
rtn = []
if self.num_nodes == 2:
cmd = 'hd'
elif self.num_nodes == 1:
cmd = 'pd1'
else:
raise Exception("Unexpected number of nodes")
if self.f is None:
if phase2 is None:
rtn.append("pd1 {} {} {}".format(name, "$fs", self.nodeName))
if cmd == 'pd1':
rtn.append("{} {} {} {}".format(cmd, name, "$fs", *self.nodeName))
elif cmd == 'hd':
rtn.append("{} {} 180 {} {}".format(cmd, name, *self.nodeName))
else:
raise Expection("Unexpected command")
else:
rtn.append("pd1 {} {} {} {}".format(name, "$fs", phase2, self.nodeName))
if self.num_nodes != 1: raise Exception("Can't do RF homodyne output")
rtn.append("pd1 {} {} {} {}".format(name, "$fs", phase2, self.nodeName[0]))
else:
if self.num_nodes != 1: raise Exception("Can't do RF homodyne output")
if quad not in ("I", "Q"):
raise pkex.BasePyKatException("quadrature must be 'I' or 'Q'")
phase = self.IQ_phase(quad, self.phase)
if phase2 is None:
rtn.append("pd2 {} {} {} {} {}".format(name, self.f , phase, "$fs", self.nodeName))
rtn.append("pd2 {} {} {} {} {}".format(name, self.f , phase, "$fs", self.nodeName[0]))
else:
rtn.append("pd2 {} {} {} {} {} {}".format(name, self.f, phase, "$fs", phase2, self.nodeName))
rtn.append("pd2 {} {} {} {} {} {}".format(name, self.f, phase, "$fs", phase2, self.nodeName[0]))
rtn.extend(self._pdtype(name, sigtype))
return rtn
def get_diff_cmds(self, quad='I', sigtype='z'):
'''
......@@ -1734,14 +1809,14 @@ class Output(object):
name = self.get_transfer_name(quad=quad, sigtype=sigtype)
rtn = []
if self.f is None:
rtn.append("pd {} {}".format(name, self.nodeName))
rtn.append("pd {} {}".format(name, self.nodeName[0]))
else:
if quad not in ("I", "Q"):
raise pkex.BasePyKatException("quadrature must be 'I' or 'Q'")
phase = self.IQ_phase(quad, self.phase)
rtn.append("pd1 {} {} {} {}".format(name, self.f, phase, self.nodeName))
rtn.append("pd1 {} {} {} {}".format(name, self.f, phase, self.nodeName[0]))
rtn.extend(self._pdtype(name, sigtype))
......
......@@ -749,22 +749,11 @@ def make_kat(name="design", katfile=None, verbose = False, debug=False, use_RF_D
The `name` argument selects from default aLIGO files included in Pykat:
- design: A file based on the design parameters for the final aLIGO setup.
125W input, T_SRM = 20%.
- design_low_power: A file based on the design parameters for the final aLIGO setup.
20W input, T_SRM = 35%. The higher SRM transmission mirror is used for low power
operation. 20W input power from O1 observation.
- design_with_IMC_HAM2: A file based on `design` but has the IMC and HAM2 blocks
which contain design parameter input optics
- design_with_IMC_HAM2_FI_OMC: A file with the OMC and IMC, most complete file
keepComments: If true it will keep the original comments from the file
preserveComments: If true it will keep the const commands in the kat
"""
names = ['design', 'design_low_power', 'design_with_IMC_HAM2', 'design_with_IMC_HAM2_FI_OMC']
names = ['design']
if debug:
kat = finesse.kat(tempdir=".",tempname="test")
......@@ -834,13 +823,15 @@ def make_kat(name="design", katfile=None, verbose = False, debug=False, use_RF_D
# If we don't have an OMC then we need to attach
# directly to the AS node. Otherwise use OMC refl
if "OMC" in kat.getBlocks():
nAS_RF = ["nOMC_ICb","nAS"] # Output() class accepts list of node names and match to the first one it finds
alt_nAS = "nAS"
nAS_RF = "nOMC_ICb"
else:
alt_nAS = None
nAS_RF = "nAS"
kat.IFO.AS_f1 = Output(kat.IFO, "AS_f1", nAS_RF, "f1", phase=101)
kat.IFO.AS_f2 = Output(kat.IFO, "AS_f2", nAS_RF, "f2", phase=14)
kat.IFO.AS_f36 = Output(kat.IFO, "AS_f36", nAS_RF, "f36M", phase=14)
kat.IFO.AS_f1 = Output(kat.IFO, "AS_f1", nAS_RF, "f1", phase=101, alternate_node_name=alt_nAS)
kat.IFO.AS_f2 = Output(kat.IFO, "AS_f2", nAS_RF, "f2", phase=14, alternate_node_name=alt_nAS)
kat.IFO.AS_f36 = Output(kat.IFO, "AS_f36", nAS_RF, "f36M", phase=14, alternate_node_name=alt_nAS)
kat.IFO.AS_DC = Output(kat.IFO, "AS_DC", "nAS")
kat.IFO.POW_BS = Output(kat.IFO, "PowBS", "nPRBS*")
......@@ -1126,62 +1117,6 @@ def pretune_SRCL(_kat, verbose=False, debug=False):
deg = (out.x[SRCL_metric.argmin()] % 360) - 180
IFO.preSRCL.apply_tuning(-deg)
def setup(base, old=True, DC_offset_pm=20, verbose=False):
"""
Runs a preparation routine to produce a LIGO model at a resonable operating point.
This uses the pretune2 and pretune_SRCL methods which allow you
"""
assert_aligo_ifo_kat(base)
base = base.deepcopy()
base.verbose = False
base.removeBlock('locks', False)
base.removeBlock('errsigs', False)
base.removeBlock('powers', False)
base.phase = 2
base.IFO.fix_mirrors()
kat = base.deepcopy()
kat.IFO.remove_modulators()
if old:
pretune(kat,pretune_precision=1e-3, verbose=verbose)
else:
pretune_2(kat, pretune_precision=1e-3, verbose=verbose)
# Apply the tunings to our base kat file
base.IFO.apply_tunings(kat.IFO.get_tunings())
base.IFO.adjust_PRC_length(verbose=verbose)
if verbose:
pretune_status(base)
base.IFO.lengths_status()
# Set DC offset and plot
DCoffset = DC_offset_pm*1e-12 / base.lambda0 * 180.0
base.IFO.set_DC_offset(DCoffset=DCoffset, verbose=verbose)