Commit 6db85759 authored by Daniel Brown's avatar Daniel Brown

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

parents df890ca8 83b15c39
......@@ -773,23 +773,33 @@ def c2r_errsig(z, phase):
return np.real(z*np.exp(-1j*phase*np.pi/180.0))
def isCavsStable(kat):
isStable = True
is_stable = True
kat1 = kat.deepcopy()
code = ''
cp_names = []
for cav in kat1.getAll(pykat.commands.cavity):
code += 'cp {0} x m\ncp {0} y m\n'.format(cav)
cp_names.append('{0}_x_m'.format(cav))
cp_names.append('{0}_y_m'.format(cav))
code += 'noxaxis'
kat1.parse(code)
out = kat1.run()
ymax = np.abs(out.y).max()
cp_outs = []
for name in cp_names:
cp_outs.append(out[name])
ymax = np.abs(cp_outs).max()
if ymax >= 1:
isStable = False
return isStable, ymax
is_stable = False
return is_stable, ymax
def opt_demod_phase(cdata, x, xbounds=None, err_tol=1e-5, xatol=1e-9, isplot=False):
def opt_demod_phase(cdata, x, xbounds=None, err_tol=1e-5, xatol=1e-9, isplot=False,
verbose=False):
'''
Optimizes the demodulation phase of a complex error signal generated by Finesse.
Demands that the error signal is smaller than err_tol to compute the slope or
......@@ -841,7 +851,7 @@ def opt_demod_phase(cdata, x, xbounds=None, err_tol=1e-5, xatol=1e-9, isplot=Fal
#print('hej!')
return -dy1
sol_2 = minimize_scalar(find_op, bounds=xbounds, method='bounded',
options={'maxiter': 500, 'disp': True, 'xatol': xatol})
options={'maxiter': 500, 'disp': verbose, 'xatol': xatol})
# print(sol_2)
return sol_2.fun
......@@ -852,7 +862,7 @@ def opt_demod_phase(cdata, x, xbounds=None, err_tol=1e-5, xatol=1e-9, isplot=Fal
demod_phase_bounds = slice(min_phase, max_phase, step)
phase_bounds = (demod_phase_bounds,)
# Searching for optimal demodulation phase.
sol = brute(max_og, ranges=phase_bounds, full_output=True, finish=fmin, disp=False)
sol = brute(max_og, ranges=phase_bounds, full_output=True, finish=fmin, disp=verbose)
# Reading out results
demod_phase = sol[0][0]
optical_gain = -sol[1]*yscale/xscale
......
......@@ -1119,6 +1119,24 @@ class ADV_IFO(IFO):
for k,v in DCoffset.items():
self.set_DC_offset(DCoffset=v, offset_type=k, verbose=False)
def _strToDOFs(self, DOFs):
dofs = []
for _ in DOFs:
if isinstance(_, six.string_types):
if _ in self.DOFs:
dofs.append(self.DOFs[_])
else:
raise pkex.BasePyKatException(
"Could not find DOF called `%s`. Possible DOF options: %s" % (
_, str(list(self.DOFs.keys()))))
else:
raise pkex.BasePyKatException(
"'%s' not possible DOF options: %s" % (
_, str(list(self.DOFs.keys()))))
return dofs
def assert_adv_ifo_kat(kat):
if not isinstance(kat.IFO, ADV_IFO):
......
......@@ -125,7 +125,7 @@ def error_signals(_kat, xlimits=[-1,1], DOFs=None, plotDOFs=None, replaceDOFSign
toShow = None
if plotDOFs is not None:
toShow = self._strToDOFs(plotDOFs)
toShow = kat.IFO._strToDOFs(plotDOFs)
# Check if other DOF signals we need to include for plotting
for _ in toShow:
......@@ -171,9 +171,11 @@ def error_signals(_kat, xlimits=[-1,1], DOFs=None, plotDOFs=None, replaceDOFSign
else:
for _ in toShow:
if legend is None:
legend = _.name
_legend = _.name
else:
_legend = legend
ax.plot(out.x, out[_.signal_name()] - DC_Offset, label=legend)
ax.plot(out.x, out[_.signal_name()] - DC_Offset, label=_legend)
ax.set_xlim([np.min(out.x), np.max(out.x)])
ax.set_xlabel("{} [deg]".format(d.name))
......
......@@ -322,6 +322,25 @@ class ALIGO_IFO(IFO):
kat.lp1.L += delta_l
kat.IFO.compute_derived_lengths(kat)
def adjust_SRC_length(self, verbose=False):
"""
Adjust SRC length so that it fulfils the requirement
lSRC = (M/5) * c/(2*f1), see [1] equation C.1
In the current design M=17 and N=3.4.
This function directly alters the lengths of the associated kat object.
"""
kat = self.kat
vprint(verbose, "-- adjusting SRC length")
ltmp = 0.5 * clight / kat.IFO.f1
delta_l = 3.4 * ltmp - kat.IFO.lSRC
vprint(verbose, " adusting kat.ls1.L by {:.4g}m".format(delta_l))
kat.ls1.L += delta_l
kat.IFO.compute_derived_lengths(kat)
def apply_lock_feedback(self, out, idx=None):
"""
......
......@@ -323,6 +323,24 @@ class VOYAGER_IFO(IFO):
kat.lp1.L += delta_l
kat.IFO.compute_derived_lengths(kat)
def adjust_SRC_length(self, verbose=False):
"""
Adjust SRC length so that it fulfils the requirement
lSRC = (M/5) * c/(2*f1), see [1] equation C.1
In the current design M=17 and N=3.4.
This function directly alters the lengths of the associated kat object.
"""
kat = self.kat
vprint(verbose, "-- adjusting SRC length")
ltmp = 0.5 * clight / kat.IFO.f1
delta_l = 3.4 * ltmp - kat.IFO.lSRC
vprint(verbose, " adusting kat.ls1.L by {:.4g}m".format(delta_l))
kat.ls1.L += delta_l
kat.IFO.compute_derived_lengths(kat)
def apply_lock_feedback(self, out, idx=None):
"""
......
import numpy as np
def all_bp_detectors(_base, param='q', direction='x'):
"""Generates a string of kat code consisting of
beam parameter detectors for every node in `_base`.
Parameters
----------
_base : :class:`.kat`
Base object containing the optical configuration.
param : str, optional
Parameter of the `bp` detector to detect (defaults to 'q').
direction : str, optional
Plane of detection ('x' for tangential, 'y' for sagittal,
"both" for generating `bp` detectors for both planes).
Returns
-------
out : str
A string of the kat code for the beam parameter detectors.
"""
code = ""
if direction == "both": direction = ['x', 'y']
else: direction = [direction]
for node in _base.nodes.getNodes().keys():
for d in direction:
code += "bp bp{n} {d} {p} {n}\n".format(
n=node, d=d, p=param
)
return code
def generate_amp_detectors(node, maxtem, f=0, ad_prefix="ad"):
"""Generates a string of kat code consisting of amplitude detectors
for each TEM up to a given maxtem.
Parameters
----------
node : str
Name of the node at which to place the amplitude detectors.
maxtem : int
Maximum mode order up to which to generate amplitude detectors.
f : double
Frequency offset to carrier [Hz].
ad_prefix : str
Prefix of amplitude detector names. These will be named
as ``[<ad_prefix>nm, ...]``.
Returns
-------
out : str
A string of the kat code for the amplitude detectors.
author: George Smetana
date: 04/12/2018
"""
temsize = maxtem + 1
code = ""
for n in range(temsize):
for m in range(temsize - n):
code += "ad {prefix}{n}{m} {n} {m} {f} {node} \n".format(
prefix=ad_prefix, n=n, m=m, f=f, node=node
)
return code
def get_tem_array(out, maxtem, ad_prefix='ad', power=False, order_only=True, table=False, plot=False):
"""
Return array of TEM mode amplitudes/powers obtained from amplitude detectors at a node.
Parameters
----------
out : :class:`.KatRun`
Output from a Finesse simulation.
maxtem : int
Maximum TEM order to which detectors should be generated.
ad_prefix : str, optional
Prefix of amplitude detector names.
power : bool, optional
If `True` returns the TEM powers, else returns the TEM amplitudes.
order_only : bool, optional
If `True`, returns quadrature sum of modes of the same TEM order,
else returns the full 2D array.
table : bool, optional
If `True`, print the output using the ``tabulate`` module.
plot : bool, optional
If `True`, generates a plot of the output.
Returns
-------
If `order_only == True` and noaxis set:
1D array where array[k] is amp/power in kth TEM order
If `order_only == True` and xaxis set:
2D array where array[k, i] is amp/power in kth TEM order of ith xaxis measurement.
If `order_only == False` and noaxis set:
2D array where array[n, m] is amp/power in TEM_nm mode
If `order_only == False` and xaxis set:
3D array where array[n, m, i] is amp/power in TEM_nm mode of ith xaxis measurement.
author: George Smetana
date: 04/12/2018
"""
from tabulate import tabulate
if plot:
import matplotlib as mpl
import matplotlib.pyplot as plt
temsize = maxtem + 1
# Depth is how many measurements each detector makes, this must be set manually to 1 for special
# case of a single measurement with noxaxis set.
depth = 1 if out.x.shape is () else out.x.size
# Initialise the amplitude array with all nans
tem_amps = np.empty((temsize, temsize, depth))
tem_amps[:] = np.nan
for n in range(temsize):
for m in range(temsize - n):
tem_amps[n, m, :] = out["{prefix}{n}{m}".format(
prefix=ad_prefix, n=n, m=m
)]
if depth > 1 and (table or plot):
# If multiple measurements are made, printing a table/plotting doesn't make sense.
raise Exception("Table printing and plotting is not supported "
"for multiple outputs per detector.")
if order_only:
orders = list(range(temsize))
tem_powers = tem_amps*np.conj(tem_amps)
x = np.array(orders)
y = x[::-1]
# Rotate array so that diagonal elements (elements of the same TEM order) are stacked in columns.
rotated_powers = tem_powers[x, np.add.outer(x, y) - (temsize - 1), :]
order_powers = np.nansum(rotated_powers, axis=1)
order_vals = np.sqrt(order_powers) if not power else order_powers
if depth == 1:
# If only one measurement is made per detector, keeping a 2D array is redundant.
order_vals = order_vals[:, 0]
if table and depth == 1:
value_name = ['Power'] if power else ['Amplitude']
headers = ['Order'] + orders
row = value_name + list(order_vals)
print(tabulate([row], headers=headers, tablefmt='grid'))
if plot and depth == 1:
plt.figure()
ax = plt.axes()
ax.scatter(orders, order_vals)
ax.set_xticks(orders)
ax.grid(True, 'both', 'both')
ax.set_yscale('log')
ax.set_ylim((min(order_vals[order_vals > 0])*0.8, max(order_vals)*1.2))
ax.set_xlabel('TEM modes order')
ax.set_ylabel('Power' if power else 'Amplitude')
ax.set_title('{} of Each TEM Order'.format('Power' if power else 'Amplitude'))
plt.show()
return order_vals
else:
if depth == 1:
# If only one measurement is made per detector, keeping a 3D array is redundant.
tem_amps = tem_amps[:, :, 0]
tem_vals = tem_amps*np.conj(tem_amps) if power else tem_amps
orders = list(range(temsize))
if table and depth == 1:
value_name = 'Power' if power else 'Amplitude'
print('{} of All Modes - Location (n, m) for TEM_nm'.format(value_name))
print(tabulate(tem_vals, tablefmt='grid'))
if plot and depth == 1:
m_list = orders * temsize
n_list = sorted(m_list)
plt.figure()
ax = plt.axes()
grid = ax.scatter(n_list, m_list, c=tem_vals.flatten(), norm=mpl.colors.LogNorm())
plt.colorbar(grid)
ax.set_xticks(orders)
ax.set_yticks(orders)
ax.grid(True, 'major', 'both')
ax.set_xlabel('n index in $TEM_{nm}$')
ax.set_ylabel('m index in $TEM_{nm}$')
ax.set_title('{} of Each TEM Mode'.format('Power' if power else 'Amplitude'))
plt.show()
return tem_vals
\ No newline at end of file
Markdown is supported
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