diff --git a/examples/asc_test/master5.py b/examples/asc_test/master5.py new file mode 100644 index 0000000000000000000000000000000000000000..e2a73a428886d9bcc2c76179f616c6fd30b52541 --- /dev/null +++ b/examples/asc_test/master5.py @@ -0,0 +1,391 @@ +from pykat import finesse +from pykat.commands import * +from pykat.utilities.optics.gaussian_beams import gauss_param + +import pylab as pl +import scipy +from scipy.optimize import minimize_scalar +import numpy as np +import shelve +import copy +import sys + + +def main(): + + print """ + -------------------------------------------------------------- + Example file for using PyKat to automate Finesse simulations + Finesse: http://www.gwoptics.org/finesse + PyKat: https://pypi.python.org/pypi/PyKat/ + + The file runs through the various pykat files which are used + to generate the Finesse results reported in the document: + `Comparing Finesse simulations, analytical solutions and OSCAR + simulations of Fabry-Perot alignment signals', LIGO-T1300345 + + This file is part of a collection. Run this after master2.py + + Andreas Freise 06.12.2013 + -------------------------------------------------------------- + """ + + # shall we clear the workspace? + # %reset -f + + # making these global during testing and debugging + global kat + global out + + kat = finesse.kat(tempdir=".",tempname="test") + kat.verbose = False + + tmpresultfile = 'myshelf2.dat' + + # loading data saved by master.py + kat.loadKatFile('asc_base3.kat') + try: + tmpfile = shelve.open(tmpresultfile) + result=tmpfile['result'] + tmpfile.close() + except: raise Exception("Could not open temprary results file {0}".format(tmpresultfile)) + + + # overwriting some variables + kat.maxtem=3 + Lambda=1064.0e-9 + + # this does not work yet due to the scale command + kat.PDrefl_p.enabled = False + kat.PDrefl_q.enabled = False + kat.WFS1_I.enabled = False + kat.WFS1_Q.enabled = False + kat.WFS2_I.enabled = False + kat.WFS2_Q.enabled = False + + kat.ETM.phi=result['phi_tuned'] + + + (beam1, beam2, beam3) = get_qs(kat) + print " Measured beam parameter:" + print " - At front of ITM (no thermal lens) q={0}".format(beam1.q) + print " (eqals w0={0}, z={1})".format(beam1.w0, beam1.z) + print " - At pick of mirror 'po' (50k lens) q={0}".format(beam2.q) + print " (eqals w0={0}, z={1})".format(beam2.w0, beam2.z) + print " - At pick of mirror 'po' (5k lens) q={0}".format(beam3.q) + print " (eqals w0={0}, z={1})".format(beam3.w0, beam3.z) + #print " Setting these now view Gauss command and adding thermal lens" + + kat.ITM.nITM1.node.setGauss(kat.ITM,beam1) + + print "--------------------------------------------------------" + print " 11. computing beam sizes with thermal lens" + beam_size(kat, beam2, beam3) + + kat.ITM_TL.f=50e3 + kat.maxtem = 8 + print "--------------------------------------------------------" + print " 11. computing beam tilt with thermal lens (f={0}, maxtem={1})".format(kat.ITM_TL.f, kat.maxtem) + #gravity_tilt(kat) + + kat.ITM_TL.f=5e3 + kat.maxtem = 23 + + print "--------------------------------------------------------" + print " 12. computing beam tilt with thermal lens (f={0}, maxtem={1})".format(kat.ITM_TL.f, kat.maxtem) + #gravity_tilt(kat) + + + print "--------------------------------------------------------" + print " 12. compute beam center with thermal lens" + + + + + print "--------------------------------------------------------" + print " Saving results in temp. files to be read by master6.py" + tmpkatfile = "asc_base4.kat" + tmpresultfile = "myshelf3.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) + # now the result variables: + tmpfile = shelve.open(tmpresultfile) + tmpfile['result']=result + tmpfile.close() + + +#----------------------------------------------------------------------------------- + +def get_qs(tmpkat): + kat = copy.deepcopy(tmpkat) + nodename0="npsl" + nodename1="nITM1" + nodename2="nWFS1" + nodename3="nWFS2" + # measure beam parameter for the 'cold beam' i.e. the laser beam + # matched to the cavity without any thermal lens + code_bp = "bp w0 y q {0}\nbp w1 y q {1}\nbp w2 y q {2}\nbp w3 y q {3}".format(nodename0,nodename1,nodename2,nodename3) + + kat.parseKatCode(code_bp) + kat.parseKatCode('yaxis re:im') + kat.noxaxis = True + kat.maxtem=0 + + def carrier_size(tmpkat, f): + kat = copy.deepcopy(tmpkat) + kat.ITM_TL.f=f + out = kat.run(printout=0,printerr=0) + raw_input("Press enter to continue") + # computing beam size at WFS1 and WFS2 + q2 = complex(out['w2'][0],out['w2'][1]) + beam1 = gauss_param(q=q2) + q3 = complex(out['w3'][0],out['w3'][1]) + beam2 = gauss_param(q=q3) + print " Carrier (cavity eigenmode) beam size with thermal lens f={0}".format(f) + print " - WFS1 w={0}".format(beam1.w) + print " (q={0}, w0={1}, z={2})".format(beam1.q, beam1.w0, beam1.z) + print " - WFS2 w={0}".format(beam2.w) + print " (q={0}, w0={1}, z={2})".format(beam2.q, beam2.w0, beam2.z) + + return(beam1, beam2) + + + def sideband_size(tmpkat, f): + kat = copy.deepcopy(tmpkat) + # 1. run finesse with input laser mode matched to cavity (no thermal lens) + out = kat.run(printout=0,printerr=0) + raw_input("Press enter to continue") + # beam at laser when matched to cold cavity + # (note the sign flip of the real part to change direction of gauss param) + q0 = complex(-1.0*out['w0'][0],out['w0'][1]) + beam0 = gauss_param(q=q0) + kat.psl.npsl.node.setGauss(kat.psl, beam0) + kat.parseKatCode("startnode npsl") + + # add thermal lens and propagate input beam to ITM + kat.ITM_TL.f=f + out = kat.run(printout=0,printerr=0) + raw_input("Press enter to continue") + + # computing beam size at ITM an set it as new startnode + q1 = complex(out['w1'][0],out['w1'][1]) + beam1 = gauss_param(q=q1) + kat.ITM.nITM1.node.setGauss(kat.ITM, beam1) + kat.removeLine("startnode") + kat.parseKatCode("startnode nITM1") + out = kat.run(printout=0,printerr=0) + + # computing beam size at WFS1 and WFS2 + q2 = complex(out['w2'][0],out['w2'][1]) + beam2 = gauss_param(q=q2) + q3 = complex(out['w3'][0],out['w3'][1]) + beam3 = gauss_param(q=q3) + print " Sideband (input mode) beam size with thermal lens f={0}".format(f) + print " - WFS1 w={0}".format(beam2.w) + print " (q={0}, w0={1}, z={2})".format(beam2.q, beam2.w0, beam2.z) + print " - WFS2 w={0}".format(beam2.w) + print " (q={0}, w0={1}, z={2})".format(beam3.q, beam3.w0, beam3.z) + + return(beam1, beam2, beam3) + + f=50e3 + carrier_size(kat,f) + sideband_size(kat,f) + + f=5e3 + carrier_size(kat,f) + (beam1,beam2,beam3)=sideband_size(kat,f) + + return (beam1, beam2,beam3) + +def asc_signal(tmpkat): + kat = copy.deepcopy(tmpkat) + + code_lock = """ + set err PDrefl_p re + lock z $err 900 1p + put* ETM phi $z + noplot z + """ + + kat.parseKatCode(code_lock) + kat.parseKatCode('yaxis abs') + kat.noxaxis = True + kat.maxtem=1 + + signal=np.zeros((2, 2)) + kat.ITM.ybeta=1e-10 + kat.ETM.ybeta=0.0 + out = kat.run(printout=0,printerr=0) + WFS1_idx=out.ylabels.index("WFS1_I") + WFS2_idx=out.ylabels.index("WFS2_I") + signal[0,0] = out.y[WFS1_idx] + signal[1,0] = out.y[WFS2_idx] + + kat.ITM.ybeta=0.0 + kat.ETM.ybeta=-1e-10 + out = kat.run(printout=0,printerr=0) + signal[0,1] = out.y[WFS1_idx] + signal[1,1] = out.y[WFS2_idx] + signal = signal *1e10 + sensors=('WFS1', 'WFS2') + mirrors=('ITM', 'ETM') + print " ASC Matrix:" + for i in range(2): + print " ", sensors[i], " ", + for j in range(2): + print "%12.10g" % signal[i,j], + print mirrors[i] + return signal + + +def asc_phases(tmpkat): + kat = copy.deepcopy(tmpkat) + + kat.parseKatCode('yaxis abs') + kat.noxaxis = True + kat.maxtem=1 + + def demod_phase1(x): + kat.WFS1_I.phi[0]=x + out = kat.run(printout=0,printerr=0) + WFS1_idx=out.ylabels.index("WFS1_I") + signal = out.y[WFS1_idx] + print '\r minimising: function value %g ' % signal , + sys.stdout.flush() + return -1*abs(signal) + + def demod_phase2(x): + kat.WFS2_I.phi[0]=x + out = kat.run(printout=0,printerr=0) + WFS2_idx=out.ylabels.index("WFS2_I") + signal = out.y[WFS2_idx] + print '\r minimising: function value %g ' % signal , + sys.stdout.flush() + return -1*abs(signal) + + kat.ITM.ybeta=1e-10 + kat.ETM.ybeta=0.0 + res = minimize_scalar(demod_phase1, method='brent') + WFS1_phase = res.x + print "" + print " WFS1 demod phase : %.10g deg" % WFS1_phase + + kat.ITM.ybeta=0.0 + kat.ETM.ybeta=-1e-10 + res = minimize_scalar(demod_phase2, method='brent') + WFS2_phase = res.x + print "" + print " WFS2 demod phase : %.10g deg" % WFS2_phase + return(WFS1_phase, WFS2_phase) + +def gravity_tilt(tmpkat): + kat = copy.deepcopy(tmpkat) + + def compute_gravity_tilt(tmpkat): + kat = copy.deepcopy(tmpkat) + out = kat.run(printout=0,printerr=0) + + y1 = out["b1"] + y2 = out["b1_1k"] + # shift of beam center on detector 1 (as m/w0y) + x1 = np.sum(out.x*y1)/np.sum(y1) + # shift of beam center on detector 2 (as m/w0y) + x2 = np.sum(out.x*y2)/np.sum(y2) + # calibrate this in meter by mutliplying with w0y + # and compute the angle geometrically + w0=out["w0y"][0] + detector_distance = 1000.0 + tilt=w0*(x2-x1)/detector_distance + print " Wavefront tilt : %g nrad" % tilt + + code_WFS1 = """ + beam b1 nWFS1 + beam b1_1k nL1_in + bp w0y y w0 nWFS1 + """ + + code_WFS2 = """ + m md 0 1 0 nWFS2 nWFS2b + s sd 1k nWFS2b nWFS2c + beam b1 nWFS2* + beam b1_1k nWFS2c + bp w0y y w0 nWFS2 + """ + + code_xaxis= """ + xaxis b1 y lin -40 40 800 + put b1_1k y $x1 + yaxis abs + """ + print " WFS1:" + print " ITM ybeta 0.1nm" + kat.parseKatCode(code_WFS1) + kat.parseKatCode(code_xaxis) + kat.spo1.L=1000.0 + kat.ITM.ybeta=1e-10 + kat.ETM.ybeta=0.0 + compute_gravity_tilt(kat) + print " ETM ybeta -0.1nm" + kat.ITM.ybeta=0.0 + kat.ETM.ybeta=-1e-10 + compute_gravity_tilt(kat) + + print " WFS1:" + print " ITM ybeta 0.1nm" + kat = copy.deepcopy(tmpkat) + kat.parseKatCode(code_WFS2) + kat.parseKatCode(code_xaxis) + kat.spo1.L=1.0e-9 + kat.ITM.ybeta=1e-10 + kat.ETM.ybeta=0.0 + compute_gravity_tilt(kat) + print " ETM ybeta -0.1nm" + kat.ITM.ybeta=0.0 + kat.ETM.ybeta=-1e-10 + compute_gravity_tilt(kat) + +def beam_size(tmpkat, beam2, beam3): + kat = copy.deepcopy(tmpkat) + + global out + code_bps = """ + bp wWFS1 y w nWFS1 + bp wWFS2 y w nWFS2 + """ + kat.parseKatCode(code_bps) + kat.maxtem = 0 + kat.ITM.R=1.0 + kat.ITM.T=0.0 + + kat.noxaxis = True + kat.ITM_TL.f=50e3 + kat.po.nWFS1.node.setGauss(kat.po,beam2) + + out = kat.run(printout=0,printerr=0) + + WFS1_idx=out.ylabels.index("wWFS1") + WFS2_idx=out.ylabels.index("wWFS2") + + + y1 = out.y[WFS1_idx] + y2 = out.y[WFS2_idx] + print " Beam size with thermal lens f={0}".format(kat.ITM_TL.f) + print " WFS1: {0}cm".format(y1*100.0) + print " WFS2: {0}cm".format(y2*100.0) + + kat.ITM_TL.f=5e3 + kat.po.nWFS1.node.setGauss(kat.po,beam3) + out = kat.run(printout=0,printerr=0) + y1 = out.y[WFS1_idx] + y2 = out.y[WFS2_idx] + print " Beam size with thermal lens f={0}".format(kat.ITM_TL.f) + print " WFS1: {0}cm".format(y1*100.0) + print " WFS2: {0}cm".format(y2*100.0) + + +if __name__ == '__main__': + main() + diff --git a/pykat/node_network.py b/pykat/node_network.py index 9d48e3cf6ac242dc072ff0bf4ae73d91f592f40c..d5e006f1ce7c022499efd978c722c57e7e52869b 100644 --- a/pykat/node_network.py +++ b/pykat/node_network.py @@ -266,9 +266,11 @@ class Node(object): rtn = [] if self.__q_x == self.__q_y: - rtn.append("gauss* g_{node} {comp} {node} {z} {zr}".format(node=self.name, comp=self.__q_comp.name, z=self.__q_x.real, zr=self.__q_x.imag)) + rtn.append("gauss g_{node} {comp} {node} {w0:.15g} {z:.15g}".format(node=self.name, comp=self.__q_comp.name, w0=self.__q_x.w0, z=self.__q_x.z)) + #rtn.append("gauss* g_{node} {comp} {node} {z} {zr}".format(node=self.name, comp=self.__q_comp.name, z=self.__q_x.real, zr=self.__q_x.imag)) else: - rtn.append("gauss* g_{node} {comp} {node} {zx} {zrx} {zy} {zry}".format(node=self.name, comp=self.__q_comp.name, zx=self.__q_x.real, zrx=self.__q_x.imag, zy=self.__q_y.real, zry=self.__q_y.imag)) + rtn.append("gauss g_{node} {comp} {node} {w0x:.15g} {zx:.15g} {w0y:.15g} {zy:.15g}".format(node=self.name, comp=self.__q_comp.name, w0x=self.__q_x.w0, zx=self.__q_x.z, w0y=self.__q_y.w0, zy=self.__q_y.z)) + #rtn.append("gauss* g_{node} {comp} {node} {zx} {zrx} {zy} {zry}".format(node=self.name, comp=self.__q_comp.name, zx=self.__q_x.real, zrx=self.__q_x.imag, zy=self.__q_y.real, zry=self.__q_y.imag)) return rtn