swithcing from gauss* to gauss (more intuitive), also

added proper formatting to the gauss parameters i.e.
{.15g}, probably need to do that for many other params as well.

examples/asc_test/master5.py

+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()
+

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