master2.py 8.93 KB
Newer Older
1
2
3
4
5
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

6
7
8
from pykat import finesse
from pykat.commands import *
import pylab as pl
9
import scipy
Andreas Freise's avatar
Andreas Freise committed
10
from scipy.optimize import fmin
11
import numpy as np
12
13
import shelve
import copy
14
15
import sys

16
17

def main():
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
	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 5 and 6 and saves temporary
	data and a new Finesse input file to be read by master3.py,
	and master4.py.
	
	Andreas Freise 16.01.2014
	--------------------------------------------------------------
	""")
	
	# 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
	
49
	tmpresultfile = "myshelf1.dat"
50
51
52
53
	
	# loading data saved by master.py
	kat.loadKatFile('asc_base2.kat')
	try:
54
55
		tmpfile = shelve.open(tmpresultfile, flag="c")
		result=tmpfile[str('result')]
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
		tmpfile.close()
	except: raise Exception("Could not open temprary results file {0}".format(tmpresultfile))
		
	# overwriting some variables
	kat.maxtem=3
	Lambda=1064.0e-9

	# disable PDH photo diode as we won't need it for most of this
	kat.PDrefl_p.enabled = False
	kat.PDrefl_q.enabled = False

	# simulating a tuned cavity
	kat.ETM.phi=result['phi_tuned']
	
	print("--------------------------------------------------------")
	print(" 5. checking wavefronts for ITM/ETM tilt of 0.1nrad")
	tilt(kat)

	print("--------------------------------------------------------")
	print(" 6. compute beam tilt from beam propogation")
	gravity_tilt(kat)

	print("--------------------------------------------------------")
	print(" 7. compute optimal demodulation phase of WFS1 and WFS2")

	# adding wave front sensors to global kat object, will need them later
	# on as well.
	
	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.phi1=WFS1_phase
	kat.WFS1_Q.phi1=WFS1_phase+90.0
	kat.WFS2_I.phi1=WFS2_phase
	kat.WFS2_Q.phi1=WFS2_phase+90.0
	result['WFS1_phase']=WFS1_phase
	result['WFS2_phase']=WFS2_phase

	print("--------------------------------------------------------")
	print(" 8. compute ASC signal matrix at WFS1 and WFS2")
	signal = asc_signal(kat)
	
	print("--------------------------------------------------------")
	print(" Saving results in temp. files to be read by master3.py")
	# re-enable PDH photo diode for savinf
	kat.PDrefl_p.enabled = True
	kat.PDrefl_q.enabled = True

	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)
129
	tmpfile[str('result')]=result
130
	tmpfile.close()
131
132


133
#-----------------------------------------------------------------------------------
134

135
def asc_signal(tmpkat):
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
	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)
	# need to re-enable the photo diode for lock
	kat.PDrefl_p.enabled = True

	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()
	signal[0,0] = out["WFS1_I"]
	signal[1,0] = out["WFS2_I"]

	
	kat.ITM.ybeta=0.0
	kat.ETM.ybeta=-1e-10
	out = kat.run()
	signal[0,1] = out["WFS1_I"]
	signal[1,1] = out["WFS2_I"]
	signal = signal *1e10
	sensors=('WFS1', 'WFS2')
	mirrors=('ITM', 'ETM')
	print("	 ASC Matrix:")
	for i in range(2):
		print("	 ", sensors[i], " ", end=' ')
		for j in range(2):
			print("%12.10g" % signal[i,j], end=' ')
		print(mirrors[i])
	return signal
	
	
178
def asc_phases(tmpkat):
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
210
211
212
213
214
215
	kat = copy.deepcopy(tmpkat)
	
	kat.parseKatCode('yaxis abs')
	kat.noxaxis = True
	kat.maxtem=1

	def demod_phase1(x):
		kat.WFS1_I.phi1=x[0]
		out = kat.run()
		signal = out["WFS1_I"]
		print('\r minimising: function value %g					   ' % signal, end=' ')
		sys.stdout.flush()
		return -1*abs(signal)

	def demod_phase2(x):
		kat.WFS2_I.phi1=x[0]
		out = kat.run()
		signal = out["WFS2_I"]
		print('\r minimising: function value %g					   ' % signal, end=' ')
		sys.stdout.flush()
		return -1*abs(signal)

	kat.ITM.ybeta=1e-10
	kat.ETM.ybeta=0.0
	res = fmin(demod_phase1, [0.0], xtol=1e-8, disp=False)
	WFS1_phase = res[0]
	print("")
	print(" WFS1 demod phase : %.10g deg" % WFS1_phase)
	 
	kat.ITM.ybeta=0.0
	kat.ETM.ybeta=-1e-10
	res = fmin(demod_phase2, [0.0], xtol=1e-8, disp=False)
	WFS2_phase = res[0]
	print("")
	print(" WFS2 demod phase : %.10g deg" % WFS2_phase)
	return(WFS1_phase, WFS2_phase)	  
	
216
def gravity_tilt(tmpkat):
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
	kat = copy.deepcopy(tmpkat)

	def compute_gravity_tilt(tmpkat):
		kat = copy.deepcopy(tmpkat)
		out = kat.run()

		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(" WFS2:")
	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)

	
283
def tilt(tmpkat):
284
285
286
287
288
289
290
291
292
	kat = copy.deepcopy(tmpkat)
	
	def compute_tilt(tmpkat):
		global kat
		kat = copy.deepcopy(tmpkat)
		global out
		out = kat.run()

		# compute data x range in meters
293
		beamsize = out["w0y"][0].real
294
295
296
297
298
		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
299
300
		diff_l = (out["PDrefl_low"].real-out["PDrefl_car"].real)/stepsize
		diff_u = (out["PDrefl_up"].real-out["PDrefl_car"].real)/stepsize
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
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
		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 = """
	beam PDrefl_car 0 nWFS2
	beam PDrefl_up 9M nWFS2
	beam PDrefl_low -9M nWFS2
	bp w0y y w0 nWFS2
	"""
	code_comm = """
	xaxis PDrefl_car y lin -1 1 100
	put PDrefl_up y $x1
	put PDrefl_low y $x1
	yaxis abs:deg
	"""

	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)

	print(" ETM ybeta -0.1nm")
	kat.ITM.ybeta=0.0
	kat.ETM.ybeta=-1e-10
	(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)

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

	return 
	
356
if __name__ == '__main__':
357
	main()
358