Commit d639a43b authored by Daniel Brown's avatar Daniel Brown

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

parents 4a3b014d bca3b43a
......@@ -18,7 +18,7 @@ def run(tmpkat):
# function for root finding
def PD_q_test(x):
kat.PDrefl_q.phi[0]=x
kat.PDrefl_q.phase[0]=x
out = kat.run(printout=0,printerr=0)
print '\r root finding: function value %g ' % out.y,
sys.stdout.flush()
......
......@@ -4,7 +4,7 @@ from __future__ import print_function
from __future__ import unicode_literals
from pykat import finesse
from pykat.commands import *
#from pykat.commands import *
import copy
import pickle
import sys
......@@ -12,83 +12,83 @@ import scipy.optimize
def main():
print("""
--------------------------------------------------------------
Example file for using PyKat to automate Finesse simulations
Finesse: http://www.gwoptics.org/finesse
PyKat: http://www.gwoptics.org/pykat
The file runs through the various Finesse simulations
to generate the Finesse results reported in the document:
`Comparing Finesse simulations, analytical solutions and OSCAR
simulations of Fabry-Perot alignment signals', LIGO-T1300345,
freely available online: http://arxiv.org/abs/1401.5727
This file is part of a collection; it outputs the results
shown the document's sections 3 and 4 and saves temporary
data and a new Finesse input file to be read by master2.py.
Andreas Freise 16.01.2014
--------------------------------------------------------------
""")
# for debugging we might need to see the temporay file:
global kat
kat = finesse.kat(tempdir=".",tempname="test")
kat.verbose = False
kat.loadKatFile('asc_base.kat')
kat.maxtem=3
Lambda=1064.0e-9
result = {}
# defining variables as global for debugging
#global kat, out, result
print("--------------------------------------------------------")
print(" 1. tunes ETM position to find resonance")
kat.ETM.phi=resonance(kat)
print("--------------------------------------------------------")
print(" 2. print sideband and carrier powers/amplitudes")
powers(kat)
print("--------------------------------------------------------")
print(" 3. determine the optimal phase for the PDH signal")
(result['p_phase'], result['q_phase']) = pd_phase(kat)
# setting demodulation phase
code_det = """
pd1 PDrefl_p 9M 0 nWFS1
scale 2 PDrefl_p
pd1 PDrefl_q 9M 90 nWFS1
scale 2 PDrefl_q
"""
kat.parseKatCode(code_det)
kat.PDrefl_p.phase1=result['p_phase']
kat.PDrefl_q.phase1=result['q_phase']
print("--------------------------------------------------------")
print(" 4. adding a 0.1nm offset to ETM and compute PDH signal")
result['phi_tuned']=float(kat.ETM.phi)
result['phi_detuned'] = result['phi_tuned'] + 0.1*360.0/1064.0
kat.ETM.phi=result['phi_detuned']
print(" new ETM phi tuning = %g " % kat.ETM.phi)
(result['pd_p'], result['pd_q']) = pd_signal(kat)
print(" PDH inphase = %e " % result['pd_p'])
print(" PDH quadrtature = %e " % result['pd_q'])
print("--------------------------------------------------------")
print(" Saving results in temp. files to be read by master2.py")
tmpkatfile = "asc_base2.kat"
tmpresultfile = "myshelf1.dat"
print(" kat object saved in: {0}".format(tmpkatfile))
print(" current results saved in: {0}".format(tmpresultfile))
# first the current kat file
kat.saveScript(tmpkatfile)
with open(tmpresultfile, 'wb') as handle:
pickle.dump(result, handle)
print("""
--------------------------------------------------------------
Example file for using PyKat to automate Finesse simulations
Finesse: http://www.gwoptics.org/finesse
PyKat: http://www.gwoptics.org/pykat
The file runs through the various Finesse simulations
to generate the Finesse results reported in the document:
`Comparing Finesse simulations, analytical solutions and OSCAR
simulations of Fabry-Perot alignment signals', LIGO-T1300345,
freely available online: http://arxiv.org/abs/1401.5727
This file is part of a collection; it outputs the results
shown the document's sections 3 and 4 and saves temporary
data and a new Finesse input file to be read by master2.py.
Andreas Freise 16.01.2014
--------------------------------------------------------------
""")
# for debugging we might need to see the temporay file:
global kat
kat = finesse.kat(tempdir=".",tempname="test")
kat.verbose = False
kat.loadKatFile('asc_base.kat')
kat.maxtem=3
Lambda=1064.0e-9
result = {}
# defining variables as global for debugging
#global kat, out, result
print("--------------------------------------------------------")
print(" 1. tunes ETM position to find resonance")
kat.ETM.phi=resonance(kat)
print("--------------------------------------------------------")
print(" 2. print sideband and carrier powers/amplitudes")
powers(kat)
print("--------------------------------------------------------")
print(" 3. determine the optimal phase for the PDH signal")
(result['p_phase'], result['q_phase']) = pd_phase(kat)
# setting demodulation phase
code_det = """
pd1 PDrefl_p 9M 0 nWFS1
scale 2 PDrefl_p
pd1 PDrefl_q 9M 90 nWFS1
scale 2 PDrefl_q
"""
kat.parseKatCode(code_det)
kat.PDrefl_p.phase1=result['p_phase']
kat.PDrefl_q.phase1=result['q_phase']
print("--------------------------------------------------------")
print(" 4. adding a 0.1nm offset to ETM and compute PDH signal")
result['phi_tuned']=float(kat.ETM.phi)
result['phi_detuned'] = result['phi_tuned'] + 0.1*360.0/1064.0
kat.ETM.phi=result['phi_detuned']
print(" new ETM phi tuning = %g " % kat.ETM.phi)
(result['pd_p'], result['pd_q']) = pd_signal(kat)
print(" PDH inphase = %e " % result['pd_p'])
print(" PDH quadrtature = %e " % result['pd_q'])
print("--------------------------------------------------------")
print(" Saving results in temp. files to be read by master2.py")
tmpkatfile = "asc_base2.kat"
tmpresultfile = "myshelf1.dat"
print(" kat object saved in: {0}".format(tmpkatfile))
print(" current results saved in: {0}".format(tmpresultfile))
# first the current kat file
kat.saveScript(tmpkatfile)
with open(tmpresultfile, 'wb') as handle:
pickle.dump(result, handle)
#---------------------------------------------------------------------------
def pd_signal(tmpkat):
......@@ -107,100 +107,100 @@ def pd_signal(tmpkat):
def pd_phase(tmpkat):
kat = copy.deepcopy(tmpkat)
code_det = """
pd1 PDrefl_q 9M 90 nWFS1
"""
kat.parseKatCode(code_det)
kat.noxaxis= True
# function for root finding
def PD_q_test(x):
kat.PDrefl_q.phase1=x
out = kat.run()
print('\r root finding: function value {0:<16g}'.format(float(out.y)), end='')
sys.stdout.flush()
return float(out.y)
# do root finding
xtol=1e-8
(result, info)=scipy.optimize.bisect(PD_q_test,80.0,100.0, xtol=xtol, maxiter=500, full_output=True)
print("")
if info.converged:
p_phase=result-90.0
q_phase=result
print(" Root has been found:")
print(" p_phase %8f" % (p_phase))
print(" q_phase %8f" % (q_phase))
print(" (%d iterations, %g tolerance)" % (info.iterations, xtol))
return (p_phase, q_phase)
else:
raise Exception("Root has not been found")
kat = copy.deepcopy(tmpkat)
code_det = """
pd1 PDrefl_q 9M 90 nWFS1
"""
kat.parseKatCode(code_det)
kat.noxaxis= True
# function for root finding
def PD_q_test(x):
kat.PDrefl_q.phase1=x
out = kat.run()
print('\r root finding: function value {0:<16g}'.format(float(out.y)), end='')
sys.stdout.flush()
return float(out.y)
# do root finding
xtol=1e-8
(result, info)=scipy.optimize.bisect(PD_q_test,80.0,100.0, xtol=xtol, maxiter=500, full_output=True)
print("")
if info.converged:
p_phase=result-90.0
q_phase=result
print(" Root has been found:")
print(" p_phase %8f" % (p_phase))
print(" q_phase %8f" % (q_phase))
print(" (%d iterations, %g tolerance)" % (info.iterations, xtol))
return (p_phase, q_phase)
else:
raise Exception("Root has not been found")
def powers(tmpkat):
kat = copy.deepcopy(tmpkat)
code1 = """
ad EOM_up 9M nEOM1
ad EOM_low -9M nEOM1
pd cav_pow nITM2
ad cav_c 0 nITM2
ad WFS1_u 9M nWFS1
ad WFS1_l -9M nWFS1
ad WFS1_c 0 nWFS1
ad WFS2_u 9M nWFS2
ad WFS2_l -9M nWFS2
ad WFS2_c 0 nWFS2
noxaxis
"""
kat.parseKatCode(code1)
global out
out = kat.run()
for i in range(len(out.y[0])):
print(" %8s: %.4e" % (out.ylabels[i], out.y[0,i]))
kat = copy.deepcopy(tmpkat)
code1 = """
ad EOM_up 9M nEOM1
ad EOM_low -9M nEOM1
pd cav_pow nITM2
ad cav_c 0 nITM2
ad WFS1_u 9M nWFS1
ad WFS1_l -9M nWFS1
ad WFS1_c 0 nWFS1
ad WFS2_u 9M nWFS2
ad WFS2_l -9M nWFS2
ad WFS2_c 0 nWFS2
noxaxis
"""
kat.parseKatCode(code1)
global out
out = kat.run()
for i in range(len(out.y[0])):
print(" %8s: %.4e" % (out.ylabels[i], out.y[0,i]))
def resonance(tmpkat):
kat = copy.deepcopy(tmpkat)
code1 = """
ad carr2 0 nITM1*
ad carr3 0 nITM2
yaxis deg
"""
kat.parseKatCode(code1)
kat.noxaxis = True
# function for root finding
def carrier_resonance(x):
kat.ETM.phi=x
out = kat.run()
phase = (out.y[0,0]-out.y[0,1]-90)%360-180.0
print('\r root finding: function value {0:<16g}'.format(float(phase)), end='')
sys.stdout.flush()
return phase
# do root finding
xtol=1e-8
(result, info)=scipy.optimize.bisect(carrier_resonance,0.0,40.0, xtol=xtol, maxiter=500, full_output=True)
print("")
if info.converged:
print(" Root has been found:")
print(" ETM phi %8f" % (result))
print(" (%d iterations, %g tolerance)" % (info.iterations, xtol))
return result
else:
raise Exception(" Root has not been found")
kat = copy.deepcopy(tmpkat)
code1 = """
ad carr2 0 nITM1*
ad carr3 0 nITM2
yaxis deg
"""
kat.parseKatCode(code1)
kat.noxaxis = True
# function for root finding
def carrier_resonance(x):
kat.ETM.phi=x
out = kat.run()
phase = (out.y[0,0]-out.y[0,1]-90)%360-180.0
print('\r root finding: function value {0:<16g}'.format(float(phase)), end='')
sys.stdout.flush()
return phase
# do root finding
xtol=1e-8
(result, info)=scipy.optimize.bisect(carrier_resonance,0.0,40.0, xtol=xtol, maxiter=500, full_output=True)
print("")
if info.converged:
print(" Root has been found:")
print(" ETM phi %8f" % (result))
print(" (%d iterations, %g tolerance)" % (info.iterations, xtol))
return result
else:
raise Exception(" Root has not been found")
if __name__ == '__main__':
main()
main()
......@@ -3,11 +3,16 @@
# Firstly you need to start an ipython cluster on your computer. To do this open
# a new terminal and type the command:
#
# ipcluster start -n 4
# ipcluster start --n=4
#
# (or something similar)
# This will start a cluster with 4 workers. You should set this number to how many
# cores you have.
#
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
import pykat
from pykat.parallel import parakat
......
This diff is collapsed.
......@@ -11,6 +11,8 @@ from scipy.optimize import brute
from scipy.optimize import fmin
from scipy.interpolate import interp1d
from scipy.optimize import minimize_scalar
from scipy.optimize import minimize
from scipy.misc import comb
# THE PLOTTING OPTION WILL BE REMOVED
import matplotlib.pyplot as plt
......@@ -101,9 +103,9 @@ def round_to_n(x, n):
factor = (10 ** power)
return round(x * factor) / factor
def vprint(verbose, printstr):
def vprint(verbose, printstr, end='\n'):
if verbose:
print(printstr)
print(printstr,end=end)
def BS_optical_path(thickness, n=1.44963098985906, angle=45.0):
"""
......@@ -599,11 +601,117 @@ def mismatch_scan_L(base, node, length, lower, upper, steps):
return out['qx'], out['qy']
def modematch(kat, components, cavs, node, verbose = False):
'''
Mode matches the cavity eigenmmodes for the cavities in cavs by varying the
components in components. Computes mode overlaps between the cavity eigenmodes.
Minimises the maximum cavity mismatch.
Inputs
--------
kat - kat-object to run
components - list of names of components to vary. Each entry must be a lens (varies f),
mirror (Rc), BS (Rc) or space (L).
cavs - list with names of cavities to match
node - name of node where to compute mismatches. Must be first or
last node in IFO for reliable result.
verbose - If true, prints the new optimised paramaters
Returns
--------
kat1 - kat-object with the optimised component parameters. Deepcopy of kat.
'''
Nc = len(cavs)
Np = len(components)
kat1 = kat.deepcopy()
kat2 = kat.deepcopy()
# Storing parameters to tune, and their initial values, in lists
attrs = []
attrs2 = []
p0 = []
for c in components:
if isinstance(kat1.components[c], pykat.components.lens):
p0.append(kat1.components[c].f.value)
elif (isinstance(kat1.components[c], pykat.components.mirror) or
isinstance(kat1.components[c], pykat.components.beamSplitter)):
p0.append(kat1.components[c].Rc.value)
elif isinstance(kat1.components[c], pykat.components.space):
p0.append(kat1.components[c].L.value)
attrs.append(kat1.components[c])
attrs2.append(kat2.components[c])
# Switching off cavity commands for cavities not in cavs
for cav in kat1.getAll(pykat.commands.cavity):
if not cav.name in cavs:
cav.remove()
# Cost function
def func(p):
for k in range(Np):
if isinstance(attrs[k], pykat.components.lens):
attrs[k].f = p[k]
elif isinstance(attrs[k], pykat.components.space):
attrs[k].L = p[k]
elif (isinstance(attrs[k], pykat.components.mirror) or
isinstance(attrs[k], pykat.components.beamSplitter)):
attrs[k].Rc = p[k]
mmx, mmy, qs = pykat.ifo.mismatch_cavities(kat1, node)
mm = np.zeros(comb(Nc,2,exact=True), dtype=float)
cs = deepcopy(cavs)
k = 0
for c1 in cavs:
cs.pop(0)
for c2 in cs:
mm[k] = np.sqrt(mmx[c1][c2])*np.sqrt(mmy[c1][c2])
k += 1
# print(kat1.CPN_TL.f.value, kat1.CPW_TL.f.value, mm.mean())
return mm.max()
if verbose:
# Computing initial mismatch. Only for display
mmx, mmy, qs = pykat.ifo.mismatch_cavities(kat1, node)
mm0 = np.zeros(comb(Nc,2,exact=True), dtype=float)
cs = deepcopy(cavs)
k = 0
for c1 in cavs:
cs.pop(0)
for c2 in cs:
mm0[k] = np.sqrt(mmx[c1][c2])*np.sqrt(mmy[c1][c2])
k += 1
# Optimising
opts = {'xtol': 1.0, 'ftol': 1.0e-7, 'disp': False}
out = minimize(func, p0, method='Nelder-Mead', options=opts)
if not out['success']:
pkex.printWarning(out.message)
# Setting new parameters to kat-object
for k in range(Np):
if isinstance(attrs2[k], pykat.components.lens):
attrs2[k].f = out.x[k]
elif isinstance(attrs2[k], pykat.components.space):
attrs2[k].L = out.x[k]
elif (isinstance(attrs2[k], pykat.components.mirror) or
isinstance(attrs2[k], pykat.components.beamSplitter)):
attrs2[k].Rc = out.x[k]
if verbose:
print('Maximum mismatch: {:.2e} --> {:.2e}'.format(mm0.max(), out.fun))
for c in components:
if isinstance(kat.components[c], pykat.components.lens):
print(' {}.f: {:.5e} m --> {:.5e} m'.format(c, kat.components[c].f.value, kat2.components[c].f.value ))
elif isinstance(kat.components[c], pykat.components.space):
print(' {}.L: {:.5e} m --> {:.5e} m'.format(c, kat.components[c].L.value, kat2.components[c].L.value ))
elif (isinstance(kat.components[c], pykat.components.mirror) or
isinstance(kat.components[c], pykat.components.beamSplitter)):
print(' {}.Rc = {:.5e} m --> {:.5e} m'.format(c, kat.components[c].Rc.value, kat2.components[c].Rc.value ))
return kat2, out.x
......
......@@ -493,7 +493,7 @@ class ADV_IFO(IFO):
else:
raise pkex.BasePyKatException("\033[91m offset_type must be DARM or MICH. \033[0m")
print("-- applying user-defined DC offset to {}:".format(offset_type))
vprint(verbose, "-- applying user-defined DC offset to {}:".format(offset_type))
_kat = self.kat
m = self.mirrors
......@@ -533,11 +533,11 @@ class ADV_IFO(IFO):
out = kat.run()
print("-- adjusting {} DCoffset based on light in dark port:".format(offset_type))
vprint(verbose, "-- adjusting {} DCoffset based on light in dark port:".format(offset_type))
waste_light = round(float(out[signame]),1)
print(" waste light in AS port of {:2} W".format(waste_light))
vprint(verbose, " waste light in AS port of {:2} W".format(waste_light))
#kat_lock = _kat.deepcopy()
DCoffset = self.find_DC_offset(5*waste_light, offset_type, verbose=verbose)
......@@ -551,6 +551,7 @@ class ADV_IFO(IFO):
This function directly alters the tunings of the associated kat object.
"""
m = self.mirrors
if offset_type == 'DARM' or offset_type == 'darm':
isDARM = True
......@@ -603,7 +604,7 @@ class ADV_IFO(IFO):
#kat.NI.phi = IXphi - phi/2.0
out = kat.run()
print(" ! ", out[self.B1.get_signal_name()], phi)
# print(verbose, " ! ", out[self.B1.get_signal_name()], phi)
return np.abs(out[self.B1.get_signal_name()] - AS_power)
......@@ -863,7 +864,131 @@ class ADV_IFO(IFO):
for _ in inspect.getmembers(self, lambda x: isinstance(x, Output)):
self.Outputs[_[0]] = _[1]
def thermal_lensing(self, thermal_mirror_list):
out = compute_thermal_effect(self.kat, thermal_mirror_list, nScale=True)
mirrors = self.kat.IFO.mirrors
# Setting values to kat-object
for k,v in out.items():
# Setting new RoC (No RoC calculations yet)
# self.kat.components[k].Rc = v[1][0]
# Setting new lens
if k == mirrors['IX']:
self.kat.CPN_TL.f = float(v['f_CP_new'])
elif k == mirrors['IY']:
self.kat.CPW_TL.f = float(v['f_CP_new'])
def find_maxtem(self, tol=5e-3, start = 0, stop = 10, verbose=False):
'''
Finding the minimum required maxtem for the power to converge to within the relative tolerance tol.
'''
kat1 = self.kat.deepcopy()
sigs = []
sigs.append(kat1.IFO.POW_BS.add_signal())
sigs.append(kat1.IFO.POW_X.add_signal())
sigs.append(kat1.IFO.POW_Y.add_signal())
kat1.parse('noxaxis\nyaxis abs')
run = True
P_old = np.zeros([2,3], dtype=float) + 1e-9
P = np.zeros(3, dtype=float) + 1e-9
rdiff = np.ones([2,3],dtype=float)
mxtm = 0
while run and mxtm <= stop:
vprint(verbose, mxtm, ": ")
P_old[0,:] = P_old[1,:]
P_old[1,:] = P
rdiff_old = rdiff[1,:]
kat1.maxtem = mxtm
out = kat1.run()
# print(out.stdout)
for k,s in enumerate(sigs):
P[k] = out[s]
rdiff = np.abs((P-P_old)/P_old)
if rdiff.max()<tol:
run = False
vprint(verbose, "{0[0]:.2e} {0[1]:.2e} {0[2]:.2e}".format(rdiff[1,:]))
# print(kat1.maxtem, rdiff, rdiff.max())