234 lines
9.0 KiB
Python
234 lines
9.0 KiB
Python
#!/usr/bin/env python
|
|
# *-----------------------------------------------------------------------*
|
|
# | |
|
|
# | Copyright (c) 2022 by Paul Scherrer Institute (http://www.psi.ch) |
|
|
# | |
|
|
# | Author Thierry Zamofing (thierry.zamofing@psi.ch) |
|
|
# *-----------------------------------------------------------------------*
|
|
'''
|
|
coordinate systems, optical center, xray axis, pixel sizes etc.
|
|
|
|
ppm= pixel per mm
|
|
update_pix2pos
|
|
|
|
modes:
|
|
0x01: update_pix2pos
|
|
0x02: update_optical_center
|
|
'''
|
|
import logging
|
|
_log=logging.getLogger(__name__)
|
|
|
|
import numpy as np
|
|
from scipy.optimize import fsolve
|
|
|
|
class VLSgrating:
|
|
|
|
def __init__(self):
|
|
# NEED TO CALL SETUP for parameters: _param and _equ
|
|
pass
|
|
|
|
vlsgTypes=('80meV','80meV_k=2','80meV_k=3','gratingTest1')
|
|
|
|
def setup(self,vlsType='80meV'):
|
|
if vlsType=='80meV':
|
|
self._param={'R':65144.054, # radius of the grating mirror
|
|
'a0':1160.759,'a1':0.3409,'a2':-3.74E-5,'a3':1.98E-7} # polynome for line density of the mirror
|
|
self._equ={'k': 1, # deflection order (default=1)
|
|
'lut': # lookup table for starting guess values at a given energy for equation solver
|
|
[ # energy (eV), r1(mm), r2(mm), alpha(rad),
|
|
(300, 1193, 4127, np.deg2rad(88)), ] }
|
|
elif vlsType=='80meV_k=2':
|
|
self._param={'R':65144.054, # radius of the grating mirror
|
|
'a0':1160.759,'a1':0.3409,'a2':-3.74E-5,'a3':1.98E-7} # polynome for line density of the mirror
|
|
self._equ={'k': 2, # deflection order (default=1)
|
|
'lut': # lookup table for starting guess values at a given energy for equation solver
|
|
[ # energy (eV), r1(mm), r2(mm), alpha(rad),
|
|
(350, 1045, 4557, np.deg2rad(88)),
|
|
(500, 1398, 3763, np.deg2rad(88)), ] }
|
|
elif vlsType=='80meV_k=3':
|
|
self._param={'R':65144.054, # radius of the grating mirror
|
|
'a0':1160.759,'a1':0.3409,'a2':-3.74E-5,'a3':1.98E-7} # polynome for line density of the mirror
|
|
self._equ={'k': 3, # deflection order (default=1)
|
|
'lut': # lookup table for starting guess values at a given energy for equation solver
|
|
[ # energy (eV), r1(mm), r2(mm), alpha(rad),
|
|
(400, 988, 4778, np.deg2rad(88)),
|
|
(800, 1156, 4218, np.deg2rad(88)), ] }
|
|
elif vlsType=='gratingTest1':
|
|
self._param={'R':65144.054, # radius of the grating mirror
|
|
'a0':1360.759,'a1':0.3409,'a2':-3.74E-5,'a3':1.98E-7} # polynome for line density of the mirror
|
|
self._equ={'k': 1, # deflection order (default=1)
|
|
'lut': # lookup table for starting guess values at a given energy for equation solver
|
|
[ # energy (eV), r1(mm), r2(mm), alpha(rad),
|
|
(300, 1133, 4468, np.deg2rad(88)), ] }
|
|
else:
|
|
raise ValueError(f'unknown type:{vlsType}')
|
|
|
|
def energy2geometry(self, energy):
|
|
# mode raytrace or interpolate lut
|
|
# prec precision requitements (precision iteration)
|
|
# returns R1,R2,aa,bb,cc
|
|
pn=self.solve_focus_equ(energy)
|
|
r1, r2, alpha=pn
|
|
beta=self.beta(*pn)
|
|
gamma=self.gamma(*pn)
|
|
#geo=(r1,r2,np.rad2deg(alpha), np.rad2deg(beta), np.rad2deg(gamma))
|
|
geo={'r1':r1,'r2':r2,'aa':np.rad2deg(alpha), 'bb':np.rad2deg(beta), 'cc':np.rad2deg(gamma)}
|
|
|
|
return geo
|
|
|
|
def geometry2raw(self,r1,r2,aa,bb,cc):
|
|
|
|
# returns raw motor positions
|
|
# offset detector plane to deflected beam: 34deg
|
|
|
|
mt=gtz=gty1=gty2=grx=gtx=dtz=dty1=dty2=drx=None
|
|
degArm=90-aa+90-bb
|
|
radArm=np.deg2rad(degArm)
|
|
gtz=r1
|
|
grx=90-aa
|
|
dtz=np.cos(radArm)*r2
|
|
dty1=dty2=np.sin(radArm)*r2
|
|
drx=90-aa+90-bb+cc-34
|
|
|
|
dd=cc-34 # angle of bellow to detector
|
|
|
|
e=None
|
|
if degArm>10:
|
|
e=ValueError('angle arm > 10°')
|
|
elif degArm<1:
|
|
e=ValueError('angle arm < 1°')
|
|
elif abs(dd)>15:
|
|
e=ValueError('angle bellow to detector > 15°')
|
|
#elif drx<-1:
|
|
# e=ValueError('angle drx < -1°')
|
|
#elif drx>20:
|
|
# e=ValueError('angle drx > 20°')
|
|
|
|
#raw=(mt,gtz,gty1,gty2,grx,gtx,dtz,dty1,dty2,drx)
|
|
raw={'mt':mt,'gtz':gtz,'gty1':gty1,'gty2': gty2,'grx':grx,'gtx':gtx,'dtz':dtz,'dty1':dty1,'dty2':dty2,'drx':drx}
|
|
return raw,e
|
|
|
|
# ---------- eugenio calculation functions ----------
|
|
#def bring_to_focus(self, eV):
|
|
def solve_focus_equ(self, eV):
|
|
es=self._equ
|
|
p0=es['lut'][0][1:] # fsolve starting values for [r1, r2, alpha]
|
|
es['En']=eV # Update the dictionary with the energy at which I want to calculate the positions
|
|
es['Lam_nm']=1239.842/eV # Convert energy to wavelength
|
|
pn=fsolve(self.equations, p0, maxfev=100000000)
|
|
#es['r1'], es['r2'], es['alpha']=pn
|
|
return pn
|
|
|
|
def equations(self, pn): # system of 3 equations for 3 parameters
|
|
return (self.Eq1(*pn), self.Eq2(*pn), self.Eq4(*pn))
|
|
|
|
def Eq1(self, r1, r2, alpha):
|
|
p=self._param;equ=self._equ
|
|
a1, R, a0, k, Lam_nm=p['a1'], p['R'], p['a0'], equ['k'], equ['Lam_nm']
|
|
beta=self.beta(r1, r2, alpha)
|
|
return np.cos(alpha)**2/r1+np.cos(beta)**2/r2-\
|
|
(np.cos(alpha)+np.cos(beta))/R-a1*k*Lam_nm/1e6
|
|
|
|
def Eq2(self, r1, r2, alpha):
|
|
p=self._param;equ=self._equ
|
|
a0, a1, a2, R, k, Lam_nm=p['a0'], p['a1'], p['a2'], p['R'], equ['k'], equ['Lam_nm']
|
|
beta=self.beta(r1, r2, alpha)
|
|
return np.sin(alpha)/r1*(np.cos(alpha)**2/r1-np.cos(alpha)/R)-\
|
|
np.sin(beta)/r2*(
|
|
np.cos(beta)**2/r2-np.cos(beta)/R)+2*a2*k*Lam_nm/1e6/3
|
|
|
|
def Eq4(self, r1, r2, alpha):
|
|
p=self._param;equ=self._equ
|
|
a3, R, k, Lam_nm, a0=p['a3'], p['R'], equ['k'], equ['Lam_nm'], p['a0']
|
|
beta=self.beta(r1, r2, alpha)
|
|
return 4*np.sin(alpha)**2/r1**2*(np.cos(alpha)**2/r1-np.cos(alpha)/R)\
|
|
-1/r1*(np.cos(alpha)**2/r1-np.cos(alpha)/R)**2+1/R**2*(1/r1-np.cos(alpha)/R)\
|
|
+4*np.sin(beta)**2/r2**2*(
|
|
np.cos(beta)**2/r2-np.cos(beta)/R)\
|
|
-1/r2*(np.cos(beta)**2/r2-np.cos(beta)/R)**2\
|
|
+1/R**2*(1/r2-np.cos(beta)/R)-2*a3*k*Lam_nm/1e6
|
|
|
|
def beta(self, r1, r2, alpha): # calculate beta as a function of alpha, energy, and diffraction order (k)
|
|
p=self._param;equ=self._equ
|
|
a0, k, Lam_nm=p['a0'], equ['k'], equ['Lam_nm']
|
|
return np.arcsin(np.sin(alpha)-a0*k*Lam_nm/1e6)
|
|
|
|
def gamma(self, r1, r2, alpha): # compute detector lean angle in respect of the r2 arm
|
|
p=self._param;equ=self._equ
|
|
R, a0, a1=p['R'], p['a0'], p['a1']
|
|
beta=self.beta(r1, r2, alpha)
|
|
return np.arctan(np.cos(beta)/(2*np.sin(beta)-r2*(np.tan(beta)/R+a1/a0)))
|
|
|
|
if __name__=="__main__":
|
|
import argparse
|
|
|
|
logging.basicConfig(level=logging.DEBUG, format='%(levelname)s:%(module)s:%(lineno)d:%(funcName)s:%(message)s ')
|
|
|
|
#(h, t)=os.path.split(sys.argv[0]);cmd='\n '+(t if len(h)>20 else sys.argv[0])+' '
|
|
#exampleCmd=('', '-m0xf -v0')
|
|
epilog=__doc__ #+'\nExamples:'+''.join(map(lambda s:cmd+s, exampleCmd))+'\n'
|
|
parser=argparse.ArgumentParser(epilog=epilog, formatter_class=argparse.RawDescriptionHelpFormatter)
|
|
parser.add_argument("-m", "--mode", type=lambda x: int(x,0), help="mode (see bitmasks) default=0x%(default)x", default=0xff)
|
|
|
|
args=parser.parse_args()
|
|
_log.info('Arguments:{}'.format(args.__dict__))
|
|
|
|
logging.getLogger('matplotlib').setLevel(logging.INFO)
|
|
import numpy as np
|
|
import matplotlib as mpl
|
|
import matplotlib.pyplot as plt
|
|
import json
|
|
|
|
class MyJsonEncoder(json.JSONEncoder):
|
|
""" Special json encoder for numpy types """
|
|
|
|
def default(self, obj):
|
|
if isinstance(obj, np.integer):
|
|
return int(obj)
|
|
elif isinstance(obj, np.floating):
|
|
return float(obj)
|
|
elif isinstance(obj, np.ndarray):
|
|
return obj.tolist()
|
|
elif type(obj) not in (dict, list, str, int):
|
|
_log.error('dont know how to json')
|
|
return repr(obj)
|
|
return json.JSONEncoder.default(self, obj)
|
|
|
|
|
|
np.set_printoptions(suppress=True)
|
|
np.set_printoptions(linewidth=196)
|
|
|
|
if args.mode&0x01:
|
|
gratings=np.load('../doc/ForThierry/Standard.obj', allow_pickle='TRUE').item()
|
|
pars=gratings['Typ_80meV']['grating']
|
|
|
|
###
|
|
param=dict()
|
|
for k in ('R','a0','a1','a2','a3',): #fixed VLS parameters
|
|
param[k]=pars.pop(k) #=pars[k]
|
|
equ=dict() #equation solve parameters
|
|
for k in ('L', 'r1', 'r2', 'alpha', 'k',): # starting guess values
|
|
equ[k]=pars.pop(k) #=pars[k]
|
|
|
|
vlsg=VLSgrating()
|
|
vlsg.setup('80meV')
|
|
#vlsg._param=param
|
|
#vlsg._equ=equ
|
|
print(vlsg.energy2geometry(300))
|
|
|
|
Es=np.arange(250, 1550, 50)
|
|
#Es=[300, ]
|
|
mypars=[]
|
|
for i, eV in enumerate(Es):
|
|
pn=vlsg.solve_focus_equ(eV)
|
|
r1, r2, alpha=pn
|
|
beta=vlsg.beta(*pn)
|
|
gamma=vlsg.gamma(*pn)
|
|
mypars.append((eV, r1, r2, np.rad2deg(alpha), np.rad2deg(beta), np.rad2deg(gamma)))
|
|
|
|
mypars=np.array(mypars)
|
|
print(mypars)
|
|
|
|
#vlsg.setup(param)
|
|
print(json.dumps(param, cls=MyJsonEncoder))
|