Commit 1a53f2d4 authored by Daniel Brown's avatar Daniel Brown
Browse files

Merge branch 'master' of gitmaster.atlas.aei.uni-hannover.de:pykat/pykat

parents f37ade6a 1c4fc7f7
...@@ -33,6 +33,7 @@ def main(): ...@@ -33,6 +33,7 @@ def main():
""") """)
# for debugging we might need to see the temporay file: # for debugging we might need to see the temporay file:
global kat
kat = finesse.kat(tempdir=".",tempname="test") kat = finesse.kat(tempdir=".",tempname="test")
kat.verbose = False kat.verbose = False
kat.loadKatFile('asc_base.kat') kat.loadKatFile('asc_base.kat')
...@@ -105,7 +106,7 @@ def pd_signal(tmpkat): ...@@ -105,7 +106,7 @@ def pd_signal(tmpkat):
""" """
kat.parseKatCode(code1) kat.parseKatCode(code1)
kat.noxaxis = True kat.noxaxis = True
#global out global out
out = kat.run() out = kat.run()
print(" Cavity power: {0:.6f}W".format(out.y[0,2])) print(" Cavity power: {0:.6f}W".format(out.y[0,2]))
return (out.y[0,0], out.y[0,1]) return (out.y[0,0], out.y[0,1])
...@@ -123,8 +124,7 @@ def pd_phase(tmpkat): ...@@ -123,8 +124,7 @@ def pd_phase(tmpkat):
# function for root finding # function for root finding
def PD_q_test(x): def PD_q_test(x):
kat.PDrefl_q.phi1=x kat.PDrefl_q.phase1=x
out = kat.run() out = kat.run()
print('\r root finding: function value %g ' % out.y, end=' ') print('\r root finding: function value %g ' % out.y, end=' ')
sys.stdout.flush() sys.stdout.flush()
...@@ -132,6 +132,8 @@ def pd_phase(tmpkat): ...@@ -132,6 +132,8 @@ def pd_phase(tmpkat):
# do root finding # do root finding
xtol=1e-8 xtol=1e-8
#print("Starting values for bisect are: %e and %e \n" % (PD_q_test(60.0),PD_q_test(100.0)))
(result, info)=scipy.optimize.bisect(PD_q_test,80.0,100.0, xtol=xtol, maxiter=500, full_output=True) (result, info)=scipy.optimize.bisect(PD_q_test,80.0,100.0, xtol=xtol, maxiter=500, full_output=True)
......
...@@ -101,10 +101,10 @@ def main(): ...@@ -101,10 +101,10 @@ def main():
kat.parseKatCode(code_WFS2) kat.parseKatCode(code_WFS2)
(WFS1_phase, WFS2_phase) = asc_phases(kat) (WFS1_phase, WFS2_phase) = asc_phases(kat)
kat.WFS1_I.phi1=WFS1_phase kat.WFS1_I.phase1=WFS1_phase
kat.WFS1_Q.phi1=WFS1_phase+90.0 kat.WFS1_Q.phase1=WFS1_phase+90.0
kat.WFS2_I.phi1=WFS2_phase kat.WFS2_I.phase1=WFS2_phase
kat.WFS2_Q.phi1=WFS2_phase+90.0 kat.WFS2_Q.phase1=WFS2_phase+90.0
result['WFS1_phase']=WFS1_phase result['WFS1_phase']=WFS1_phase
result['WFS2_phase']=WFS2_phase result['WFS2_phase']=WFS2_phase
...@@ -183,7 +183,7 @@ def asc_phases(tmpkat): ...@@ -183,7 +183,7 @@ def asc_phases(tmpkat):
kat.maxtem=1 kat.maxtem=1
def demod_phase1(x): def demod_phase1(x):
kat.WFS1_I.phi1=x[0] kat.WFS1_I.phase1=x[0]
out = kat.run() out = kat.run()
signal = out["WFS1_I"] signal = out["WFS1_I"]
print('\r minimising: function value %g ' % signal, end=' ') print('\r minimising: function value %g ' % signal, end=' ')
...@@ -191,7 +191,7 @@ def asc_phases(tmpkat): ...@@ -191,7 +191,7 @@ def asc_phases(tmpkat):
return -1*abs(signal) return -1*abs(signal)
def demod_phase2(x): def demod_phase2(x):
kat.WFS2_I.phi1=x[0] kat.WFS2_I.phase1=x[0]
out = kat.run() out = kat.run()
signal = out["WFS2_I"] signal = out["WFS2_I"]
print('\r minimising: function value %g ' % signal, end=' ') print('\r minimising: function value %g ' % signal, end=' ')
......
...@@ -97,7 +97,7 @@ def asc_large(tmpkat, mir_name): ...@@ -97,7 +97,7 @@ def asc_large(tmpkat, mir_name):
done_maxtems.append(tem) done_maxtems.append(tem)
print(" Calculating maxtem = %d " % tem) print(" Calculating maxtem = %d " % tem)
kat.maxtem = tem kat.maxtem = tem
out[str(tem)] = kat.run(printout=0,printerr=1) out[str(tem)] = kat.run()
import os.path import os.path
if os.path.isfile(tmpfilename): if os.path.isfile(tmpfilename):
shutil.copyfile(tmpfilename, backupname) shutil.copyfile(tmpfilename, backupname)
......
...@@ -190,7 +190,7 @@ def get_qs(tmpkat,f): ...@@ -190,7 +190,7 @@ def get_qs(tmpkat,f):
# add thermal lens and propagate input beam to ITM # add thermal lens and propagate input beam to ITM
kat = set_thermal_lens(kat, f) kat = set_thermal_lens(kat, f)
out = kat.run(printout=0,printerr=0) out = kat.run()
# computing beam size at ITM # computing beam size at ITM
# and then we reflect of ITM, an set it as new startnode # and then we reflect of ITM, an set it as new startnode
...@@ -207,7 +207,7 @@ def get_qs(tmpkat,f): ...@@ -207,7 +207,7 @@ def get_qs(tmpkat,f):
else: else:
kat.ITM.nITM1.node.setGauss(kat.ITM, beam1) kat.ITM.nITM1.node.setGauss(kat.ITM, beam1)
kat.parseKatCode("startnode nITM1") kat.parseKatCode("startnode nITM1")
out = kat.run(printout=0,printerr=0) out = kat.run()
# computing beam size at WFS1 and WFS2 # computing beam size at WFS1 and WFS2
q2 = out['w2'] q2 = out['w2']
...@@ -228,7 +228,7 @@ def get_qs(tmpkat,f): ...@@ -228,7 +228,7 @@ def get_qs(tmpkat,f):
return [beam1, beam2, beam3, beam4] return [beam1, beam2, beam3, beam4]
global out global out
# run finesse with input laser mode matched to cavity (no thermal lens) # run finesse with input laser mode matched to cavity (no thermal lens)
out = kat.run(printout=0,printerr=0) out = kat.run()
# beam at laser when matched to cold cavity # beam at laser when matched to cold cavity
# (note the sign flip of the real part to change direction of gauss param) # (note the sign flip of the real part to change direction of gauss param)
...@@ -256,13 +256,13 @@ def asc_signal(tmpkat): ...@@ -256,13 +256,13 @@ def asc_signal(tmpkat):
signal=np.zeros((2, 2)) signal=np.zeros((2, 2))
kat.ITM.ybeta=1e-10 kat.ITM.ybeta=1e-10
kat.ETM.ybeta=0.0 kat.ETM.ybeta=0.0
out = kat.run(printout=0,printerr=0) out = kat.run()
signal[0,0] = out["WFS1_I"] signal[0,0] = out["WFS1_I"]
signal[1,0] = out["WFS2_I"] signal[1,0] = out["WFS2_I"]
kat.ITM.ybeta=0.0 kat.ITM.ybeta=0.0
kat.ETM.ybeta=-1e-10 kat.ETM.ybeta=-1e-10
out = kat.run(printout=0,printerr=0) out = kat.run()
signal[0,1] = out["WFS1_I"] signal[0,1] = out["WFS1_I"]
signal[1,1] = out["WFS2_I"] signal[1,1] = out["WFS2_I"]
signal = signal *1e10 signal = signal *1e10
...@@ -281,7 +281,7 @@ def gravity_tilt(tmpkat): ...@@ -281,7 +281,7 @@ def gravity_tilt(tmpkat):
def compute_gravity_tilt(tmpkat): def compute_gravity_tilt(tmpkat):
kat = copy.deepcopy(tmpkat) kat = copy.deepcopy(tmpkat)
out = kat.run(printout=0,printerr=0) out = kat.run()
y1 = out["b1"] y1 = out["b1"]
y2 = out["b1_1k"] y2 = out["b1_1k"]
......
...@@ -146,7 +146,7 @@ def get_qs(tmpkat): ...@@ -146,7 +146,7 @@ def get_qs(tmpkat):
def beam_size(tmpkat, f): def beam_size(tmpkat, f):
kat = copy.deepcopy(tmpkat) kat = copy.deepcopy(tmpkat)
# 1. run finesse with input laser mode matched to cavity (no thermal lens) # 1. run finesse with input laser mode matched to cavity (no thermal lens)
out = kat.run(printout=0,printerr=0) out = kat.run()
# beam at laser when matched to cold cavity # beam at laser when matched to cold cavity
# (note the sign flip of the real part to change direction of gauss param) # (note the sign flip of the real part to change direction of gauss param)
...@@ -159,7 +159,7 @@ def get_qs(tmpkat): ...@@ -159,7 +159,7 @@ def get_qs(tmpkat):
kat.ITM_TL.f=f kat.ITM_TL.f=f
if "ITM_TL_r" in kat._kat__components: if "ITM_TL_r" in kat._kat__components:
kat.ITM_TL_r.f=f kat.ITM_TL_r.f=f
out = kat.run(printout=0,printerr=0) out = kat.run()
# computing beam size at ITM # computing beam size at ITM
# and then we reflect of ITM, an set it as new startnode # and then we reflect of ITM, an set it as new startnode
...@@ -176,7 +176,7 @@ def get_qs(tmpkat): ...@@ -176,7 +176,7 @@ def get_qs(tmpkat):
else: else:
kat.ITM.nITM1.node.setGauss(kat.ITM, beam1) kat.ITM.nITM1.node.setGauss(kat.ITM, beam1)
kat.parseKatCode("startnode nITM1") kat.parseKatCode("startnode nITM1")
out = kat.run(printout=0,printerr=0) out = kat.run()
# computing beam size at WFS1 and WFS2 # computing beam size at WFS1 and WFS2
q2 = out['w2'] q2 = out['w2']
...@@ -215,7 +215,7 @@ def asc_signal(tmpkat): ...@@ -215,7 +215,7 @@ def asc_signal(tmpkat):
signal=np.zeros((2, 2)) signal=np.zeros((2, 2))
kat.ITM.ybeta=1e-10 kat.ITM.ybeta=1e-10
kat.ETM.ybeta=0.0 kat.ETM.ybeta=0.0
out = kat.run(printout=0,printerr=0) out = kat.run()
WFS1_idx=out.ylabels.index("WFS1_I") WFS1_idx=out.ylabels.index("WFS1_I")
WFS2_idx=out.ylabels.index("WFS2_I") WFS2_idx=out.ylabels.index("WFS2_I")
signal[0,0] = out.y[WFS1_idx] signal[0,0] = out.y[WFS1_idx]
...@@ -223,7 +223,7 @@ def asc_signal(tmpkat): ...@@ -223,7 +223,7 @@ def asc_signal(tmpkat):
kat.ITM.ybeta=0.0 kat.ITM.ybeta=0.0
kat.ETM.ybeta=-1e-10 kat.ETM.ybeta=-1e-10
out = kat.run(printout=0,printerr=0) out = kat.run()
signal[0,1] = out.y[WFS1_idx] signal[0,1] = out.y[WFS1_idx]
signal[1,1] = out.y[WFS2_idx] signal[1,1] = out.y[WFS2_idx]
signal = signal *1e10 signal = signal *1e10
...@@ -247,7 +247,7 @@ def asc_phases(tmpkat): ...@@ -247,7 +247,7 @@ def asc_phases(tmpkat):
def demod_phase1(x): def demod_phase1(x):
kat.WFS1_I.phi[0]=x kat.WFS1_I.phi[0]=x
out = kat.run(printout=0,printerr=0) out = kat.run()
WFS1_idx=out.ylabels.index("WFS1_I") WFS1_idx=out.ylabels.index("WFS1_I")
signal = out.y[WFS1_idx] signal = out.y[WFS1_idx]
print('\r minimising: function value %g ' % signal, end=' ') print('\r minimising: function value %g ' % signal, end=' ')
...@@ -256,7 +256,7 @@ def asc_phases(tmpkat): ...@@ -256,7 +256,7 @@ def asc_phases(tmpkat):
def demod_phase2(x): def demod_phase2(x):
kat.WFS2_I.phi[0]=x kat.WFS2_I.phi[0]=x
out = kat.run(printout=0,printerr=0) out = kat.run()
WFS2_idx=out.ylabels.index("WFS2_I") WFS2_idx=out.ylabels.index("WFS2_I")
signal = out.y[WFS2_idx] signal = out.y[WFS2_idx]
print('\r minimising: function value %g ' % signal, end=' ') print('\r minimising: function value %g ' % signal, end=' ')
...@@ -283,7 +283,7 @@ def gravity_tilt(tmpkat): ...@@ -283,7 +283,7 @@ def gravity_tilt(tmpkat):
def compute_gravity_tilt(tmpkat): def compute_gravity_tilt(tmpkat):
kat = copy.deepcopy(tmpkat) kat = copy.deepcopy(tmpkat)
out = kat.run(printout=0,printerr=0) out = kat.run()
y1 = out["b1"] y1 = out["b1"]
y2 = out["b1_1k"] y2 = out["b1_1k"]
...@@ -364,7 +364,7 @@ def beam_size(tmpkat, beam2, beam3): ...@@ -364,7 +364,7 @@ def beam_size(tmpkat, beam2, beam3):
kat.po.nWFS1.node.setGauss(kat.po,beam2) kat.po.nWFS1.node.setGauss(kat.po,beam2)
out = kat.run(printout=0,printerr=0) out = kat.run()
WFS1_idx=out.ylabels.index("wWFS1") WFS1_idx=out.ylabels.index("wWFS1")
WFS2_idx=out.ylabels.index("wWFS2") WFS2_idx=out.ylabels.index("wWFS2")
...@@ -380,7 +380,7 @@ def beam_size(tmpkat, beam2, beam3): ...@@ -380,7 +380,7 @@ def beam_size(tmpkat, beam2, beam3):
if "ITM_TL_r" in kat._kat__components: if "ITM_TL_r" in kat._kat__components:
kat.ITM_TL_r.f=5e3 kat.ITM_TL_r.f=5e3
kat.po.nWFS1.node.setGauss(kat.po,beam3) kat.po.nWFS1.node.setGauss(kat.po,beam3)
out = kat.run(printout=0,printerr=0) out = kat.run()
y1 = out.y[WFS1_idx] y1 = out.y[WFS1_idx]
y2 = out.y[WFS2_idx] y2 = out.y[WFS2_idx]
print(" Beam size with thermal lens f={0}".format(kat.ITM_TL.f)) print(" Beam size with thermal lens f={0}".format(kat.ITM_TL.f))
......
...@@ -80,7 +80,11 @@ def manualProcess(): ...@@ -80,7 +80,11 @@ def manualProcess():
# Doing this here just because the figure gains resolution. To # Doing this here just because the figure gains resolution. To
# replicate the mirror maps in the official figure measurement # replicate the mirror maps in the official figure measurement
# reports, use: smap.crop(0.0802) # reports, use: smap.crop(0.0802)
smap.crop(0.155)
# smap.crop(0.155)
# smap.crop(0.0802)
smap.crop(0.062)
smap.plot() smap.plot()
# Recentering is useful. To show the effect of xyOffset first: # Recentering is useful. To show the effect of xyOffset first:
smap.xyOffset = (0.02,0.05) smap.xyOffset = (0.02,0.05)
...@@ -88,8 +92,7 @@ def manualProcess(): ...@@ -88,8 +92,7 @@ def manualProcess():
smap.plot() smap.plot()
# And now we recentering. The origo is again at the mirror centre. # And now we recentering. The origo is again at the mirror centre.
smap.recenter() smap.recenter()
smap.plot() smap.plot()# Splitting into two versions: One where we process by convolving the mirror surface
# Splitting into two versions: One where we process by convolving the mirror surface
# with Zernike-polynomials, and one where we process by fitting surfaces to the # with Zernike-polynomials, and one where we process by fitting surfaces to the
# mirror. In the latter case Gaussian weighting is used to make the centre of the # mirror. In the latter case Gaussian weighting is used to make the centre of the
# mirror more important. # mirror more important.
...@@ -206,8 +209,8 @@ def manualProcess(): ...@@ -206,8 +209,8 @@ def manualProcess():
smap.find_radius(method='min',unit='meters'), smap.center) smap.find_radius(method='min',unit='meters'), smap.center)
amap.plot() amap.plot()
print('Aperture map created') print('Aperture map created')
filename = self.name + '_finesse.txt' filename = smap.name + '_finesse.txt'
self.write_map(filename) smap.write_map(filename)
print(' Phase map written to file {:s}'.format(filename)) print(' Phase map written to file {:s}'.format(filename))
filename = amap.name + '_aperture.txt' filename = amap.name + '_aperture.txt'
amap.write_map(filename) amap.write_map(filename)
...@@ -242,4 +245,5 @@ if isAutoProcessing: ...@@ -242,4 +245,5 @@ if isAutoProcessing:
isManualProcessing = True isManualProcessing = True
if isManualProcessing: if isManualProcessing:
manualProcess() manualProcess()
import numpy as np import numpy as np
from scipy.misc import factorial as fac from scipy.misc import factorial as fac
from six.moves import xrange
import math import math
def zernike_R(m, n, rho): def zernike_R(m, n, rho):
......
...@@ -371,7 +371,7 @@ class surfacemap(object): ...@@ -371,7 +371,7 @@ class surfacemap(object):
self.data = data self.data = data
# xlim and ylim given in centimeters # xlim and ylim given in centimeters
def plot(self, show=True, clabel=None, xlim=None, ylim=None): def plot(self, show=True, clabel=None, xlim=None, ylim=None, isBlock=False):
import pylab import pylab
if xlim is not None: if xlim is not None:
...@@ -434,7 +434,7 @@ class surfacemap(object): ...@@ -434,7 +434,7 @@ class surfacemap(object):
cbar.set_label(clabel) cbar.set_label(clabel)
if show: if show:
pylab.show() pylab.show(block=isBlock)
return fig return fig
...@@ -594,7 +594,10 @@ class surfacemap(object): ...@@ -594,7 +594,10 @@ class surfacemap(object):
nVals = range(0,n+1) nVals = range(0,n+1)
for k in nVals: for k in nVals:
mVals.append(range(-k,k+1,2)) mVals.append(range(-k,k+1,2))
elif isinstance(m,list): elif isinstance(m, list):
nVals = range(n,n+1)
mVals.append(m)
elif isinstance(m, range):
nVals = range(n,n+1) nVals = range(n,n+1)
mVals.append(m) mVals.append(m)
elif isinstance(m,int): elif isinstance(m,int):
...@@ -801,7 +804,7 @@ class surfacemap(object): ...@@ -801,7 +804,7 @@ class surfacemap(object):
return A1,xbeta,ybeta,zOff return A1,xbeta,ybeta,zOff
def rmSphericalSurf(self, Rc0, w=None, zOff=None, isCenter=[False,False]): def rmSphericalSurf(self, Rc0, w=None, zOff=None, isCenter=[False,False], maxfev=2000):
''' '''
Fits spherical surface to the mirror map and removes it. Fits spherical surface to the mirror map and removes it.
...@@ -885,7 +888,7 @@ class surfacemap(object): ...@@ -885,7 +888,7 @@ class surfacemap(object):
# Using the simplex Nelder-Mead method. This is the same or very # Using the simplex Nelder-Mead method. This is the same or very
# similar to the method used in 'FT_remove_curvature_from_mirror_map.m', # similar to the method used in 'FT_remove_curvature_from_mirror_map.m',
# but there are probably better methods to use. # but there are probably better methods to use.
opts = {'xtol': 1.0e-5, 'ftol': 1.0e-9, 'maxiter': 1000, 'maxfev': 1000, 'disp': False} opts = {'xtol': 1.0e-5, 'ftol': 1.0e-9, 'maxiter': 10000, 'maxfev': maxfev, 'disp': False}
out = minimize(costFunc, p, method='Nelder-Mead', options=opts) out = minimize(costFunc, p, method='Nelder-Mead', options=opts)
if not out['success']: if not out['success']:
msg = ' Warning: ' + out['message'].split('.')[0] + ' (nfev={:d}).'.format(out['nfev']) msg = ' Warning: ' + out['message'].split('.')[0] + ' (nfev={:d}).'.format(out['nfev'])
...@@ -1368,7 +1371,12 @@ class surfacemap(object): ...@@ -1368,7 +1371,12 @@ class surfacemap(object):
rho, phi= self.createPolarGrid() rho, phi= self.createPolarGrid()
rho = rho/R rho = rho/R
if isinstance(m,list): if isinstance(m, list):
for k in range(len(m)):
Z = zernike(m[k], n, rho, phi)
self.data[self.notNan] = self.data[self.notNan]-A[k]*Z[self.notNan]
self.zernikeRemoved = (m[k], n, A[k])
elif isinstance(m, range):
for k in range(len(m)): for k in range(len(m)):
Z = zernike(m[k], n, rho, phi) Z = zernike(m[k], n, rho, phi)
self.data[self.notNan] = self.data[self.notNan]-A[k]*Z[self.notNan] self.data[self.notNan] = self.data[self.notNan]-A[k]*Z[self.notNan]
...@@ -1555,7 +1563,7 @@ class mergedmap: ...@@ -1555,7 +1563,7 @@ class mergedmap:
for m in self.__maps: for m in self.__maps:
m.interpolate(nx, ny) m.interpolate(nx, ny)
def plot(self, mode="absorption", show=True, clabel=None, xlim=None, ylim=None, wavelength=1064e-9): def plot(self, mode="absorption", show=True, clabel=None, xlim=None, ylim=None, wavelength=1064e-9, isBlock=False):
import pylab import pylab
...@@ -1609,7 +1617,7 @@ class mergedmap: ...@@ -1609,7 +1617,7 @@ class mergedmap:
cbar.set_label(clabel) cbar.set_label(clabel)
if show: if show:
pylab.show() pylab.show(block=isBlock)
return fig return fig
...@@ -1725,7 +1733,7 @@ def read_map(filename, mapFormat='finesse', scaling=1.0e-9, mapType='phase', fie ...@@ -1725,7 +1733,7 @@ def read_map(filename, mapFormat='finesse', scaling=1.0e-9, mapType='phase', fie
Reads surface map files and returns a surfacemap object. Reads surface map files and returns a surfacemap object.
filename - name of surface map file. filename - name of surface map file.
mapFormat - 'finesse', 'ligo', 'zygo'. Currently only for ascii formats. mapFormat - 'finesse', 'ligo', 'zygo', 'metroPro' (binary).
scaling - scaling of surface height of the mirror map [m]. scaling - scaling of surface height of the mirror map [m].
''' '''
...@@ -2088,9 +2096,6 @@ def readHeaderMP(f): ...@@ -2088,9 +2096,6 @@ def readHeaderMP(f):
return hData return hData
class BinaryReaderEOFException(Exception): class BinaryReaderEOFException(Exception):
def __init__(self): def __init__(self):
pass pass
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment