Commit bca3b43a authored by Daniel Toyra's avatar Daniel Toyra

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

parents e9ce2b6f 1227471a
......@@ -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.
......@@ -59,10 +59,7 @@ class aLIGO(object):
self.kat.loadKatFile(katfile)
self.rawBlocks.read(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)
"""
# TODO different files not yet implemented
if _name != "default":
printf("aLIGO name `{}' not recognised, using 'default'",_name)
self.kat.loadKatFile(self._data_path+"design_with_IMC_HAM2_FI_OMC.kat")#"aLIGO.kat"
......@@ -265,7 +262,7 @@ put f1m f $mx1
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, precision = 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))
......@@ -276,7 +273,7 @@ put f1m f $mx1
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, precision = 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))
......@@ -285,7 +282,7 @@ put f1m f $mx1
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, precision = 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))
......@@ -294,7 +291,7 @@ put f1m f $mx1
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, precision = 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))
......@@ -302,20 +299,13 @@ put f1m f $mx1
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, precision = 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
......@@ -931,6 +921,17 @@ def set_tunings(kat, tunings):
for comp in keys:
kat.components[comp].phi = tunings[comp]
def scan_to_precision(kat, DOF, pretune_precision, minmax="max", phi=0.0, precision=60.0):
"""
find a maximum or minimum in a DOF within a given range
"""
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 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)
......
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