Files
PBSwissMX/python/helicalscan.py
2017-11-07 14:26:00 +01:00

206 lines
6.7 KiB
Python
Executable File

#!/usr/bin/env python
# *-----------------------------------------------------------------------*
# | |
# | Copyright (c) 2017 by Paul Scherrer Institute (http://www.psi.ch) |
# | |
# | Author Thierry Zamofing (thierry.zamofing@psi.ch) |
# *-----------------------------------------------------------------------*
'''
tools to setup and execute a helical scan of a cristal
#THIS IS JUST TESTING CODE TO SOLVE FINDING THE ROTATION CENTER
'''
import os, sys, json
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from utilities import *
class HelicalScan:
def __init__(self,args):
if args.cfg:
fh=open(args.cfg,'r')
s=fh.read()
cfg=json.loads(s, object_hook=ConvUtf8)
s=json.dumps(cfg, indent=2, separators=(',', ': '));print(s)
else:
fn='/tmp/shapepath4'
#fn='/home/zamofing_t/Documents/prj/SwissFEL/epics_ioc_modules/ESB_MX/python/data/'+time.strftime('%y-%m-%d-%H_%M_%S')
#cfg = {"sequencer": ['gen_grid_points(w=5,h=5,pitch=100,rnd=0.4)', 'sort_points()','gen_prog(file="'+fn+'.prg",host="SAR-CPPM-EXPMX1",mode=1,pt2pt_time=10,cnt=1)', 'plot_gather("'+fn+'.npz")']}
#cfg = {"sequencer": ['test_find_rot_ctr()']}
#cfg = {"sequencer": ['test_find_rot_ctr(n=5. ,per=1.,bias=2.31,ampl=4.12,phi=24.6)']}
cfg = {"sequencer": ['test_coord_trf()']}
self.cfg=dotdict(cfg)
self.args=args
def run(self):
print('args='+str(self.args))
print('cfg='+str(self.cfg))
#try:
# self.points=np.array(self.cfg.points)
#except AttributeError:
# pass
try:
sequencer= self.cfg.pop('sequencer')
except KeyError:
print('no command sequence to execute')
else:
dryrun=self.args.dryrun
for cmd in sequencer:
print('>'*5+' '+cmd+' '+'<'*5)
if not dryrun:
eval('self.' + cmd)
def test_coord_trf(self):
n = 3.; per = 1.; t = np.arange(n)
p=((2.3,2.31,4.12,24.6),(6.2,2.74,32.1,3.28)) #(y, bias, ampl, phi)
self.param=param=np.ndarray((len(p),5))
z=4.5 # fix z position
for i in range(2):
(y, bias, ampl, phi) =p[i]
x= ampl * np.cos(2 * np.pi * (per / n * t + phi / 360.)) + bias
print('yMeas_%d='%i+str(y)+' xMeas_%d='%i+str(x))
#param[i]=(z_i, y_i, x_i, r_i,phi_i)
param[i][0] =z
param[i][1] =y
param[i][2:]=HelicalScan.meas_rot_ctr(x) #(bias,ampl,phase)
pass
print param
#self.fwd_transform(param[0][1],0.,param[0][2],param[0][1])
# y_0 ,120deg ,x_0 ,z_0)
self.fwd_transform(param[0][1],2*np.pi/3.,param[0][2],param[0][0])
#self.fwd_transform(param[1][1],0.,param[1][2],param[1][3])
def fwd_transform(self,y,w,cx,cz):
#cx,cy: coarse stage
#TODO: NOT WORKING AT ALL NOW...
param=self.param
# param[i]=(z_i, y_i, x_i, r_i,phi_i)
p=np.ndarray((param.shape[0], 3))
for i in range(2):
#p[i][0]=param[i][2]+param[i][3]*np.cos(param[i][4]+w) # x= x_i+r_i*cos(phi_i*w)+cx
p[i][0]=cx+param[i][3]*np.cos(param[i][4]+w) # x= x_i+r_i*cos(phi_i*w)+cx
p[i][1]=param[i][1] # y= y_i
#p[i][2]=param[i][2]+param[i][3]*np.sin(param[i][4]+w) # z= z_i+r_i*sin(phi_i*w)
p[i][2] =cz + param[i][3] * np.sin(param[i][4] + w) # z= z_i+r_i*sin(phi_i*w)
print p
v=p[1]-p[0]
v=v/np.sqrt(v.dot(v)) # v/|v|
v=v*(y-param[0][1])/(param[1][1]-param[0][1]) # v(y)=v*(v-y_0)/(y_1-y_0)
v=p[0]+v
#v=v/abs(v)
print v
#
#
#
#x,y,z
#returns y,w,dx,dz
pass
def inv_transform(y,phi,dx=0,dz=0):
#dx,dy: deviation from cristal center line
#ps= #x,y,z
#returns y,phi,cx,cz
pass
@staticmethod
def meas_rot_ctr(y,per=1):
# find the amplitude bias and phase of an equidistant sampled sinus
# it needs at least 3 measurements e.g. at 0,120 240 deg or 0 90 180 270 deg
# per is the number of persiods, default is 1 period =360 deg
n=len(y)
f = np.fft.fft(y)
idx=int(per)
bias=np.absolute(f[0]/n)
phase=np.angle(f[idx])
ampl=np.absolute(f[idx])*2/n
return (bias,ampl,phase)
@staticmethod
def test_find_rot_ctr(n=3.,per=1.,bias=4.1,ampl=2.4,phi=37):
# find the rotation center, amplitude out of n (niminum 3) measurements
# n number of equidistant measurements
# per number of periods (full rotation of all measurements nut be a interger value for precise measurements)
# phi phase
# bias bias value
# ampl amplitude
t = np.arange(n)
y=ampl*np.cos(2*np.pi*(per/n*t+phi/360.))+bias
plt.figure(1)
plt.subplot(311)
plt.plot(t,y,'b.-')
plt.subplot(312)
f = np.fft.fft(y)
plt.step(t, f.real,'b.-', t, f.imag,'r.-', where='mid')
(bias,ampl,phase)=HelicalScan.meas_rot_ctr(y, per)
print('bias: '+str(bias))
print('amplitude: '+str(ampl))
print('phase: '+str(phase*360./2/np.pi))
plt.subplot(313)
t2 = np.linspace(0,2*np.pi,64)
y2=ampl*np.cos(t2+phase)+bias
plt.plot(t2,y2,'g-')
plt.stem(t/n*2*np.pi*per,y,'b-')
plt.show()
pass
if __name__=='__main__':
from optparse import OptionParser, IndentedHelpFormatter
class MyFormatter(IndentedHelpFormatter):
'helper class for formating the OptionParser'
def __init__(self):
IndentedHelpFormatter.__init__(self)
def format_epilog(self, epilog):
if epilog:
return epilog
else:
return ""
def parse_args():
'main command line interpreter function'
#usage: gpasciiCommunicator.py --host=PPMACZT84 myPowerBRICK.cfg
(h, t)=os.path.split(sys.argv[0]);cmd='\n '+(t if len(h)>3 else sys.argv[0])+' '
exampleCmd=('-n',
'-v15'
)
epilog=__doc__+'''
Examples:'''+''.join(map(lambda s:cmd+s, exampleCmd))+'\n '
fmt=MyFormatter()
parser=OptionParser(epilog=epilog, formatter=fmt)
parser.add_option('-v', '--verbose', type="int", dest='verbose', help='verbosity bits (see below)', default=0)
parser.add_option('-n', '--dryrun', action='store_true', help='dryrun to stdout')
parser.add_option('--xy', action='store_true', help='sort x,y instead y,x')
parser.add_option('--cfg', help='config file containing json configuration structure')
(args, other)=parser.parse_args()
args.other=other
sp=HelicalScan(args)
sp.run()
#------------------ Main Code ----------------------------------
#ssh_test()
ret=parse_args()
exit(ret)