master4.py 11.5 KB
Newer Older
1
2
from pykat import finesse
from pykat.commands import *
3
from pykat.optics.gaussian_beams import gauss_param
4
5
6
7

import pylab as pl
import scipy
from scipy.optimize import minimize_scalar
8
9
10
11
12
import numpy as np
import shelve
import copy
import sys

13
14
15
16
17
18
def set_thermal_lens(kat, f):
    kat.ITM_TL.f=f
    # if a bs-based cavity is used, we need to set second lens
    if "ITM_TL_r" in kat._kat__components:
        kat.ITM_TL_r.f=f
    return (kat)
19
20

def main():
21

22
23
24
25
    print """
    --------------------------------------------------------------
    Example file for using PyKat to automate Finesse simulations
    Finesse: http://www.gwoptics.org/finesse
26
27
28
    PyKat:   https://pypi.python.org/pypi/PyKat/

    The file runs through the various Finesse simulations
29
30
    to generate the Finesse results reported in the document:
    `Comparing Finesse simulations, analytical solutions and OSCAR 
31
32
    simulations of Fabry-Perot alignment signals', LIGO-T1300345,
    freely available online: http://arxiv.org/abs/1401.5727
33

34
35
36
37
38
    Run this file after master2.py to create data which can be
    plotted using master4_plot.py. Results are saved after 
    each step and plots can be created at any time.
        
    Andreas Freise 16.01.2014
39
40
41
42
43
44
45
46
47
    --------------------------------------------------------------
    """
    
    # shall we clear the workspace?
    # %reset -f

    # making these global during testing and debugging
    #global kat
    #global out
48

49
50
    kat = finesse.kat(tempdir=".",tempname="test")
    kat.verbose = False
51

52
53
54
55
56
57
58
59
60
    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))
Andreas Freise's avatar
Andreas Freise committed
61
    
62
63
    # this does not work yet due to the scale command
    kat.PDrefl_p.enabled = False
Andreas Freise's avatar
Andreas Freise committed
64
    kat.PDrefl_q.enabled = False
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
    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']

    print "--------------------------------------------------------"
    print " 11. Do beam tracing to measure beam parameters"
    # get beam parameters at nodes: "npsl", "nITM1", "nWFS1", "nWFS2", "npo2"
    beam1 = get_qs(kat,1e13)
    beam2 = get_qs(kat,50e3)
    beam3 = get_qs(kat,5e3)

    # starting with cold beam at npsl
    kat.psl.npsl.node.setGauss(kat.psl, beam1[0])
    kat.parseKatCode("startnode npsl")
Andreas Freise's avatar
Andreas Freise committed
82

83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
    global beam1, beam2, beam3

    # if we use bs-based cavity we try to set good beam
    # parameter for reflected beam, first by
    # computing 'average' beam from cold and hot beams
    x1=0.70
    x2=0.30
    if "ITM_TL_r" in kat._kat__components:
        beam50 = gauss_param(z=(x1*beam1[1].z+x2*beam2[1].z), w0=(x1*beam1[1].w0+x2*beam2[1].w0))
        beam5  = gauss_param(z=(x1*beam1[1].z+x2*beam3[1].z), w0=(x1*beam1[1].w0+x2*beam3[1].w0))
        node_text = "at ITM->nITM1r"
        t_comp=kat.ITM
        t_node=kat.ITM.nITM1r
    else:
        beam50 = gauss_param(z=(x1*beam1[4].z+x2*beam2[4].z), w0=(x1*beam1[4].w0+x2*beam2[4].w0))
        beam5  = gauss_param(z=(x1*beam1[4].z+x2*beam3[4].z), w0=(x1*beam1[4].w0+x2*beam3[4].w0))
        node_text = "at s2->npo2"
        t_comp=kat.s2
        t_node=kat.s2.npo2
    
    kat = set_thermal_lens(kat,50e3)
104
    print "--------------------------------------------------------"
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
    print " 12. computing beam tilt with thermal lens (f={0})".format(kat.ITM_TL.f)

    print " Setting compromise beam parameter {0}: w0={1}, z={2}".format(node_text, beam50.w0, beam50.z)
    t_node.node.setGauss(t_comp, beam50)
    kat.maxtem=8
    print " Calculating maxtem = %d " % kat.maxtem
    tmp = gravity_tilt(kat)

    kat = set_thermal_lens(kat,5e3)
    print "--------------------------------------------------------"
    print " 13. computing beam tilt with thermal lens (f={0})".format(kat.ITM_TL.f)
    print " Setting compromise beam parameter {0}: w0={1}, z={2}".format(node_text, beam5.w0, beam5.z)
    t_node.node.setGauss(t_comp, beam5)
    #maxtems = [1, 3, 5, 9, 11, 13, 15, 19, 23, 25, 27, 29]
    maxtems = [1, 3, 5, 7]
    converge_tilt(kat, maxtems)
   
    kat.PDrefl_p.enabled = True
    kat.WFS1_I.enabled = True
    kat.WFS2_I.enabled = True
    
    kat = set_thermal_lens(kat,50e3)
    print "--------------------------------------------------------"
    print " 12. computing alignment signal with thermal lens (f={0})".format(kat.ITM_TL.f)
    t_node.node.setGauss(t_comp, beam50)
    #maxtems = [1, 3, 5, 9, 11, 13, 15, 17, 19]
    maxtems = [1, 3, 5, 7]
    converge_asc(kat, maxtems, 'asc_signals_50.txt')

    kat = set_thermal_lens(kat,5e3)
    print "--------------------------------------------------------"
    print " 13. computing alignment signal with thermal lens (f={0})".format(kat.ITM_TL.f)
    t_node.node.setGauss(t_comp, beam5)
    #maxtems = [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31]
    maxtems = [1, 3, 5, 7]
    converge_asc(kat, maxtems, 'asc_signals_5.txt')

#-----------------------------------------------------------------------------------

def converge_asc(tmpkat, maxtems, filename):
    kat = copy.deepcopy(tmpkat)
    savedata = np.zeros([len(maxtems),5])
    for idx, tem in enumerate(maxtems):
        savedata[idx,0]=tem
        print " Calculating maxtem = %d " % tem
        kat.maxtem = tem
        tmp = asc_signal(kat)
        savedata[idx,1:5]=tmp.reshape([1,4])
        print " Saving results in file: {0}".format(filename)
        np.savetxt(filename, savedata[0:idx+1,:], fmt='%.18e', delimiter=' ')    

def converge_tilt(tmpkat, maxtems):
    kat = copy.deepcopy(tmpkat)
    savedata = np.zeros([len(maxtems),5])
    filename = "thermal_gravity.txt"
    for idx, tem in enumerate(maxtems):
        savedata[idx,0]=tem
        print " Calculating maxtem = %d " % tem
        kat.maxtem = tem
        tmp = gravity_tilt(kat)
        savedata[idx,1:5]=tmp
        print " Saving results in file: {0}".format(filename)
        np.savetxt(filename, savedata[0:idx+1,:], fmt='%.18e', delimiter=' ')    

def get_qs(tmpkat,f):
    kat = copy.deepcopy(tmpkat)

    # measure beam parameter for the 'cold beam' i.e. the laser beam
    # matched to the cavity without any thermal lens
    nodenames=["npsl", "nITM1", "nWFS1", "nWFS2", "npo2"]
    for idx, nname in enumerate(nodenames):
        kat.parseKatCode('bp w{0} y q {1}'.format(idx, nname))

    kat.parseKatCode('yaxis re:im')
    kat.noxaxis = True
    kat.maxtem=0

    def beam_size(tmpkat, f, beam0):
        kat = copy.deepcopy(tmpkat)

        kat.psl.npsl.node.setGauss(kat.psl, beam0)
        kat.parseKatCode("startnode npsl")

        # add thermal lens and propagate input beam to ITM
        kat = set_thermal_lens(kat, f)
        out = kat.run(printout=0,printerr=0)
        
        # computing beam size at ITM 
        # and then we reflect of ITM, an set it as new startnode
        q_in = complex(out['w1'][0],out['w1'][1])
195
        from pykat.optics.ABCD import apply, mirror_refl
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
        abcd = mirror_refl(1,-2500)
        q_out = apply(abcd,q_in,1,1)
        beam1 = gauss_param(q=q_out)    
        kat.removeLine("startnode")
        kat.psl.npsl.node.removeGauss()
        if "ITM_TL_r" in kat._kat__components:
            kat.ITM.nITM1r.node.setGauss(kat.ITM, beam1)
            kat.parseKatCode("startnode nITM1r")
        else:
            kat.ITM.nITM1.node.setGauss(kat.ITM, beam1)
            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)    

        # computing beam size at pick off
        q4 = complex(out['w4'][0],out['w4'][1])    
        beam4 = gauss_param(q=q4)    
        print " Input mode beam size with thermal lens f={0}".format(f)
        print " - ITM  w={0:.6}cm  (w0={1}, z={2})".format(100.0*beam1.w,beam1.w0, beam1.z)
        print " - WFS1 w={0:.6}cm  (w0={1}, z={2})".format(100.0*beam2.w,beam2.w0, beam2.z)
        print " - WFS2 w={0:.6}cm  (w0={1}, z={2})".format(100.0*beam3.w,beam3.w0, beam3.z)
        print " - npo2 w={0:.6}cm  (w0={1}, z={2})".format(100.0*beam4.w,beam4.w0, beam4.z)
        #raw_input("Press enter to continue")
        
        return [beam1, beam2, beam3, beam4]

    # run finesse with input laser mode matched to cavity (no thermal lens)
    out = kat.run(printout=0,printerr=0)

    # 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)   
    # compute beam sizes when tracing this beam back through the system
    (beam1,beam2,beam3, beam4)=beam_size(kat,f,beam0)
236
    
237
    return (beam0, beam1, beam2,beam3, beam4)
238

239
def asc_signal(tmpkat):
240
241
242
243
244
245
246
247
    kat = copy.deepcopy(tmpkat)

    code_lock = """
    set err PDrefl_p re
    lock z $err 900 1p
    put* ETM phi $z
    noplot z
    """
248
    
249
250
    kat.parseKatCode(code_lock)
    kat.parseKatCode('yaxis abs')
251
252
253
254
255
256
    kat.noxaxis = True

    signal=np.zeros((2, 2))
    kat.ITM.ybeta=1e-10
    kat.ETM.ybeta=0.0
    out = kat.run(printout=0,printerr=0)
257
258
    signal[0,0] = out["WFS1_I"]
    signal[1,0] = out["WFS2_I"]
259
260
261
262

    kat.ITM.ybeta=0.0
    kat.ETM.ybeta=-1e-10
    out = kat.run(printout=0,printerr=0)
263
264
    signal[0,1] = out["WFS1_I"]
    signal[1,1] = out["WFS2_I"]
265
266
267
268
269
270
271
272
273
274
    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
275
    
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
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
        return 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
    """
    kat.parseKatCode(code_WFS1)
    kat.parseKatCode(code_xaxis)
    kat.spo1.L=1000.0
    kat.ITM.ybeta=1e-10
    kat.ETM.ybeta=0.0
    t1=compute_gravity_tilt(kat)
    kat.ITM.ybeta=0.0
    kat.ETM.ybeta=-1e-10
    t2=compute_gravity_tilt(kat)

    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
    t3=compute_gravity_tilt(kat)
    kat.ITM.ybeta=0.0
    kat.ETM.ybeta=-1e-10
    t4=compute_gravity_tilt(kat)
    print "  WFS1 ITM {0:.4} nrad".format(t1*1e9)    
    print "  WFS2 ITM {0:.4} nrad".format(t3*1e9)    
    print "  WFS1 ETM {0:.4} nrad".format(t2*1e9)    
    print "  WFS2 ETM {0:.4} nrad".format(t4*1e9)    
    return [t1,t2,t3,t4]
341
342
343
    
if __name__ == '__main__':
    main()
344