master2.py 9.84 KB
Newer Older
1
2
3
from pykat import finesse
from pykat.commands import *
import pylab as pl
4
5
import scipy
from scipy.optimize import minimize_scalar
6
import numpy as np
7
8
import shelve
import copy
9
10
import sys

11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

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.
    
    Andreas Freise 06.12.2013
    --------------------------------------------------------------
    """
    
    # shall we clear the workspace?
    # %reset -f
Andreas Freise's avatar
Andreas Freise committed
33
34

    # making these global during testing and debugging
35
36
    #global kat
    #global out
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
    
    kat = finesse.kat(tempdir=".",tempname="test")
    kat.verbose = False
    
    tmpresultfile = 'myshelf1.dat'
    
    # loading data saved by master.py
    kat.loadKatFile('asc_base2.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
54

Andreas Freise's avatar
Andreas Freise committed
55
56
57
    # disable PDH photo diode as we won't need it for most of this
    kat.PDrefl_p.enabled = False
    kat.PDrefl_q.enabled = False
58

Andreas Freise's avatar
Andreas Freise committed
59
    # simulating a tuned cavity
Andreas Freise's avatar
Andreas Freise committed
60
    kat.ETM.phi=result['phi_tuned']
Andreas Freise's avatar
Andreas Freise committed
61
62
63
    
    print "--------------------------------------------------------"
    print " 5. checking wavefronts for ITM/ETM tilt of 0.1nrad"
64
65
66
    tilt(kat)

    print "--------------------------------------------------------"
67
    print " 6. compute beam tilt from center of gravity calculation"
68
69
    gravity_tilt(kat)

70
71
72
    print "--------------------------------------------------------"
    print " 7. compute optimal demodulation phase of WFS1 and WFS2"

Andreas Freise's avatar
Andreas Freise committed
73
74
75
    # adding wave front sensors to global kat object, will need them later
    # on as well.
    
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
    code_WFS1 = """
    pd1 WFS1_I 9M 0 nWFS1
    pdtype WFS1_I y-split
    pd1 WFS1_Q 9M 90 nWFS1
    pdtype WFS1_Q y-split
    scale 2 WFS1_I % compensate the 0.5 gain of the demodulation
    scale 2 WFS1_Q % compensate the 0.5 gain of the demodulation
    """
    code_WFS2 = """
    pd1 WFS2_I 9M 0 nWFS2
    pdtype WFS2_I y-split
    pd1 WFS2_Q 9M 90 nWFS2
    pdtype WFS2_Q y-split
    scale 2 WFS2_I % compensate the 0.5 gain of the demodulation
    scale 2 WFS2_Q % compensate the 0.5 gain of the demodulation
    """
    kat.parseKatCode(code_WFS1)
    kat.parseKatCode(code_WFS2)
    
    (WFS1_phase, WFS2_phase) = asc_phases(kat)
    kat.WFS1_I.phi[0]=WFS1_phase
    kat.WFS1_Q.phi[0]=WFS1_phase+90.0
    kat.WFS2_I.phi[0]=WFS2_phase
    kat.WFS2_Q.phi[0]=WFS2_phase+90.0
    result['WFS1_phase']=WFS1_phase
    result['WFS2_phase']=WFS2_phase
102

103
104
105
    print "--------------------------------------------------------"
    print " 8. compute ASC signal matrix at WFS1 and WFS2"
    signal = asc_signal(kat)
106
    
107
    print "--------------------------------------------------------"
108
    print " Saving results in temp. files to be read by master3.py"
Andreas Freise's avatar
Andreas Freise committed
109
110
111
112
    # re-enable PDH photo diode for savinf
    kat.PDrefl_p.enabled = True
    kat.PDrefl_q.enabled = True

113
114
115
116
117
118
119
120
121
122
123
124
    tmpkatfile = "asc_base3.kat"
    tmpresultfile = "myshelf2.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()


125
#-----------------------------------------------------------------------------------
126

127
128
129
130
131
132
133
134
135
136
137
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)
Andreas Freise's avatar
Andreas Freise committed
138
139
140
    # need to re-enable the photo diode for lock
    kat.PDrefl_p.enabled = True

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
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
    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)    
210
    
211
212
def gravity_tilt(tmpkat):
    kat = copy.deepcopy(tmpkat)
213

214
215
216
217
218
219
    def compute_gravity_tilt(tmpkat):
        kat = copy.deepcopy(tmpkat)
        out = kat.run(printout=0,printerr=0)

        y1 = out["b1"]
        y2 = out["b1_1k"]
220
        # shift of beam center  on detector 1 (as m/w0y)
221
        x1 = np.sum(out.x*y1)/np.sum(y1) 
222
        # shift of beam center  on detector 2 (as m/w0y)
223
        x2 = np.sum(out.x*y2)/np.sum(y2)
224
        # calibrate this in meter by mutliplying with w0y
225
226
227
228
229
230
231
232
233
234
235
        # 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
    """
236

237
238
239
240
241
242
243
    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
    """
244

245
246
247
248
249
250
251
252
253
254
255
256
257
    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)
258
    print " ETM ybeta -0.1nm"
259
260
261
262
263
264
265
266
267
268
269
270
271
    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)
272
    print " ETM ybeta -0.1nm"
273
274
275
276
277
    kat.ITM.ybeta=0.0
    kat.ETM.ybeta=-1e-10
    compute_gravity_tilt(kat)

    
278
def tilt(tmpkat):
279
    kat = copy.deepcopy(tmpkat)
Andreas Freise's avatar
Andreas Freise committed
280
    
281
282
283
284
    def compute_tilt(tmpkat):
        kat = copy.deepcopy(tmpkat)
        out = kat.run(printout=0,printerr=0)

285
        # compute data x range in meters
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
        beamsize = out["w0y"][0,0] 
        xrange = beamsize*(out.x.max()-out.x.min())
        stepsize=xrange/(len(out.x)-1)
        print " Beamsize %e m" % beamsize
        print " Measurement range: %e m, stepszie: %e m" % (xrange, stepsize)
        # compute difference in angle between wavefront of carrier and sidebands
        diff_l = (out["PDrefl_low"][:,1]-out["PDrefl_car"][:,1])/stepsize
        diff_u = (out["PDrefl_up"][:,1]-out["PDrefl_car"][:,1])/stepsize
        tilt_l = diff_l[1:-1]-diff_l[0:-2]
        tilt_u = diff_u[1:-1]-diff_u[0:-2]
        print " Tilt (upper  - car), mean: %e m/deg, stddev %e m/deg" % (np.mean(tilt_u), np.std(tilt_u))
        print " Tilt (lower  - car), mean: %e m/deg, stddev %e m/deg" % (np.mean(tilt_l), np.std(tilt_l))
        return (np.mean(tilt_l), np.mean(tilt_u))

    code_WFS1 = """
    beam PDrefl_car 0 nWFS1
    beam PDrefl_up 9M nWFS1
    beam PDrefl_low -9M nWFS1
    bp w0y y w0 nWFS1
    """

    code_WFS2 = """
308
309
310
311
    beam PDrefl_car 0 nWFS2
    beam PDrefl_up 9M nWFS2
    beam PDrefl_low -9M nWFS2
    bp w0y y w0 nWFS2
312
313
314
315
316
    """
    code_comm = """
    xaxis PDrefl_car y lin -1 1 100
    put PDrefl_up y $x1
    put PDrefl_low y $x1
317
318
    yaxis abs:deg
    """
Andreas Freise's avatar
Andreas Freise committed
319

320
321
322
323
324
325
326
327
    print " WFS1:"
    print " ITM ybeta 0.1nm"
    kat.parseKatCode(code_comm)
    kat.parseKatCode(code_WFS1)
    kat.ITM.ybeta=1e-10
    kat.ETM.ybeta=0.0
    (a1, a2) = compute_tilt(kat)

328
    print " ETM ybeta -0.1nm"
329
    kat.ITM.ybeta=0.0
330
    kat.ETM.ybeta=-1e-10
331
332
333
334
335
336
337
338
339
340
341
    (a3, a4) = compute_tilt(kat)
    
    print " WFS2:"
    print " ITM ybeta 0.1nm"
    kat = copy.deepcopy(tmpkat)
    kat.parseKatCode(code_comm)
    kat.parseKatCode(code_WFS2)
    kat.ITM.ybeta=1e-10
    kat.ETM.ybeta=0.0
    (a5, a6) = compute_tilt(kat)

342
    print " ETM ybeta -0.1nm"
343
    kat.ITM.ybeta=0.0
344
    kat.ETM.ybeta=-1e-10
345
346
    (a6, a7) = compute_tilt(kat)

347
    return 
348
349
350
351
    
if __name__ == '__main__':
    main()