add geometry VLS calculation

This commit is contained in:
2022-11-07 12:10:47 +01:00
parent 28822472eb
commit 6b9c6e5611
2 changed files with 364 additions and 39 deletions

280
geometry.py Normal file
View File

@@ -0,0 +1,280 @@
#!/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:
#self._param={'R':65.144m , # VLS curcature radius
# 'a0':1160.759, # gr/mm (gratings per mm)
# 'a1':0.3409, # gr/mm
# 'a2':-3.74E-5., # gr/mm^2
# 'a3':1.98E-7 # gr/mm^3
# } #spec
# gratings/mm = N(y)=a0+a1*y+a2*y^2+a3*y^3
#
#self._lut_energy2geometry
#lookuptable for starting values. The lut is a numpy array with collums:
#energy (sort asc), R1, R2, alpha, beta, gamma'
pass
def setup(R,a0,a1,a2,a3,lut):
self._param={'R':R,'a0':a0,'a1':a1,'a2':a2,'a3':a3}
self._lut_energy2geometry=np.array(lut)
pass
def setup_npz(self, fn):
#setup using a numpy npz file
pass
def save_npz(self, fn):
#save current parameters as an numpy npz file
pass
def update_lut(self, energyLst):
# recalculates the energy2geometry lut
# with the given energies
pass
def energy2geometry(self, energy,mode,prec):
# mode raytrace or interpolate lut
# prec precision requitements (precision iteration)
# returns R1,R2,aa,gg
pass
def geometry2raw(self, meas, debug=False):
# returns raw motor positions
pass
# ---------- eugenio calculation functions ----------
#def bring_to_focus(self, eV):
def solve_focus_equ(self, eV,p0=None):
es=self._equ
L, r1, alpha=es['L'], es['r1'], es['alpha'] # Take the current values from a dictionary
r2=L-r1 # starting R2 is taken as the difference between R1 and the maximum allowed length L
# p0=(np.deg2rad(88), 1650, 3500)
if not p0:
p0=[alpha, r1, r2]
p0=(np.deg2rad(88), 1650, 3500)
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['alpha'], es['r1'], es['r2']=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, alpha, r1, r2):
p=self._param;equ=self._equ
a1, R, a0, k, Lam_nm=p['a1'], p['R'], p['a0'], equ['k'], equ['Lam_nm']
return np.cos(alpha)**2/r1+np.cos(self.beta(alpha, a0))**2/r2-\
(np.cos(alpha)+np.cos(self.beta(alpha, a0)))/R-a1*k*Lam_nm/1e6
def Eq2(self, alpha, r1, r2):
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']
return np.sin(alpha)/r1*(np.cos(alpha)**2/r1-np.cos(alpha)/R)-\
np.sin(self.beta(alpha, a0))/r2*(
np.cos(self.beta(alpha, a0))**2/r2-np.cos(self.beta(alpha, a0))/R)+2*a2*k*Lam_nm/1e6/3
def Eq4(self, alpha, r1, r2):
p=self._param;equ=self._equ
a3, R, k, Lam_nm, a0=p['a3'], p['R'], equ['k'], equ['Lam_nm'], p['a0']
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(self.beta(alpha, a0))**2/r2**2*(
np.cos(self.beta(alpha, a0))**2/r2-np.cos(self.beta(alpha, a0))/R)\
-1/r2*(np.cos(self.beta(alpha, a0))**2/r2-np.cos(self.beta(alpha, a0))/R)**2\
+1/R**2*(1/r2-np.cos(self.beta(alpha, a0))/R)-2*a3*k*Lam_nm/1e6
def beta(self, alpha, a0): # calculate beta as a function of alpha, energy, and diffraction order (k)
p=self._param;equ=self._equ
k, Lam_nm=equ['k'], equ['Lam_nm']
return np.arcsin(np.sin(alpha)-a0*k*Lam_nm/1e6)
def chi(self, alpha, r1, r2): # compute detector lean angle
p=self._param;equ=self._equ
alpha, a0, r2, R, a1=equ['alpha'], p['a0'], equ['r2'], p['R'], p['a1']
return np.arctan(np.cos(self.beta(alpha, a0))/(2*np.sin(self.beta(alpha, a0))-r2*(np.tan(self.beta(alpha, a0))/R+a1/a0)))
def testFunc():
def bring_to_focus(e, pars, std_guess):
L, r1, alpha=pars['L'], pars['r1'], pars['alpha'] # Take the current values from a dictionary
r2=L-r1 # starting R2 is taken as the difference between R1 and the maximum allowed length L
if (abs(pars['En']-e)>200): # If the new energy is far from the current energy, use the user-defined starting guess
p0=std_guess
else: # otherwise, take the current values as starting guess
p0=[alpha, r1, r2]
pars['En']=e # Update the dictionary with the energy at which I want to calculate the positions
pars['Lam_nm']=1239.842/e # Convert energy to wavelength
alpha, r1, r2=fsolve(equations, p0, maxfev=100000000, args=pars)
pars['alpha'], pars['r1'], pars['r2']=alpha, r1, r2
return pars
def equations(p, pars): # system of 3 equations for 3 parameters
alpha, r1, r2=p
return (Eq1(alpha, r1, r2, pars), Eq2(alpha, r1, r2, pars), Eq4(alpha, r1, r2, pars))
def Eq1(alpha, r1, r2, pars):
a1, R, a0, k, Lam_nm=pars['a1'], pars['R'], pars['a0'], pars['k'], pars['Lam_nm']
return np.cos(alpha)**2/r1+np.cos(beta(alpha, a0, pars))**2/r2-\
(np.cos(alpha)+np.cos(beta(alpha, a0, pars)))/R-a1*k*Lam_nm/1e6
def Eq2(alpha, r1, r2, pars):
a0, a1, a2, R, k, Lam_nm=pars['a0'], pars['a1'], pars['a2'], pars['R'], pars['k'], pars['Lam_nm']
return np.sin(alpha)/r1*(np.cos(alpha)**2/r1-np.cos(alpha)/R)-\
np.sin(beta(alpha, a0, pars))/r2*(
np.cos(beta(alpha, a0, pars))**2/r2-np.cos(beta(alpha, a0, pars))/R)+2*a2*k*Lam_nm/1e6/3
def Eq4(alpha, r1, r2, pars):
a3, R, k, Lam_nm, a0=pars['a3'], pars['R'], pars['k'], pars['Lam_nm'], pars['a0']
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(alpha, a0, pars))**2/r2**2*(np.cos(beta(alpha, a0, pars))**2/r2-np.cos(beta(alpha, a0, pars))/R)\
-1/r2*(np.cos(beta(alpha, a0, pars))**2/r2-np.cos(beta(alpha, a0, pars))/R)**2\
+1/R**2*(1/r2-np.cos(beta(alpha, a0, pars))/R)-2*a3*k*Lam_nm/1e6
def beta(alpha, a0, pars): # calculate beta as a function of alpha, energy, and diffraction order (k)
k, Lam_nm=pars['k'], pars['Lam_nm']
return np.arcsin(np.sin(alpha)-a0*k*Lam_nm/1e6)
def get_chi(pars): # compute detector lean angle
alpha, a0, r2, R, a1=pars['alpha'], pars['a0'], pars['r2'], pars['R'], pars['a1']
return np.arctan(
np.cos(beta(alpha, a0, pars))/(2*np.sin(beta(alpha, a0, pars))-r2*(np.tan(beta(alpha, a0, pars))/R+a1/a0)))
gratings=np.load('../doc/ForThierry/Standard.obj', allow_pickle='TRUE').item()
pars=gratings['Typ_80meV']['grating']
#Es = np.arange(250,1550, 50)
Es=[300,]
mypars = np.zeros(5)
for i,e in enumerate(Es):
pars = bring_to_focus(e,pars,[np.deg2rad(88),1650,3500])
#pars = bring_to_focus(e,pars)
alpha, r1, r2, chi = pars['alpha'],pars['r1'],pars['r2'], get_chi(pars)
mypars = np.vstack((mypars, [e,np.rad2deg(alpha),r1,r2,np.rad2deg(chi)]))
print(mypars)
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:
testFunc()
gratings=np.load('../doc/ForThierry/Standard.obj', allow_pickle='TRUE').item()
pars=gratings['Typ_80meV']['grating']
###
#self._param={'R':65.144,'a0':1160.759,'a1':0.3409,'a2':-3.74E-5.,'a3':1.98E-7} #spec
#self._param={'R':65.144,'a0':1160.77 ,'a1':0.3409,'a2':-3.73E-5.,'a3':1.98E-7} #production param
param={'R':65.144,'a0':1160.759,'a1':0.3409,'a2':-3.74E-5,'a3':1.98E-7,
'lut': [#energy (eV), R1(mm), R2(mm), alpha(deg), beta(deg), gamma(deg)
[200 , 100., 300., .10, 20.0, 5.0],
[300 , 101., 301., .11, 20.1, 5.1],
[400 , 102., 302., .12, 20.2, 5.2],
[500 , 103., 303., .13, 20.3, 5.3],
],
}
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._param=param
vlsg._equ=equ
Es=np.arange(250, 1550, 50)
Es=[300, ]
mypars=[]
for i, eV in enumerate(Es):
pn=vlsg.solve_focus_equ(eV)
chi=vlsg.chi(*pn)
alpha, r1, r2=pn
mypars.append((eV, np.rad2deg(alpha), r1, r2, np.rad2deg(chi)))
mypars=np.array(mypars)
print(mypars)
for eV in range(300,1500,50):
alpha, r1, r2 = vlsg.solve_focus_equ(eV)
#vlsg.setup(param)
print(json.dumps(param, cls=MyJsonEncoder))

View File

@@ -53,15 +53,14 @@ def obj_info(obj,p=''):
except AttributeError:
pass
def drawArrow():
pass
class RIXSgirder(QWidget):
def __init__(self):
super().__init__()
self._param={'ctr':(200, 300), # location of big chamber
'aa':45, # RIXS arm angle
'cc':45, # UNSUSED
}
self.initUI()
@@ -70,25 +69,50 @@ class RIXSgirder(QWidget):
self.setWindowTitle('RIXS girder')
#label = QLabel('Python', self)
#label.move(50,50)
b1 = QPushButton("Button1",self)
b1.move(400,150)
#b1 = QPushButton("Button1",self)
#b1.move(400,150)
w=QtGui.QGroupBox("Drawing",self)
w.move(10,10)
l=QtGui.QVBoxLayout(w)
sld={}
for key,rng,tk,pos in (('aa',(0,180),5,20),('cc',(0,180),5,40),):
sl=QtGui.QSlider(QtCore.Qt.Horizontal)
sl.setFixedWidth(200);sl.setMinimum(rng[0]);sl.setMaximum(rng[1])
sl.setValue(self._param[key])
sl.setTickPosition(QtGui.QSlider.TicksBelow);sl.setTickInterval(tk)
l.addWidget(sl)
sl.valueChanged.connect(lambda val,key=key: self.sldChanged(key,val))
sld[key]=sl
self.show()
def sldChanged(self,key,val,*args,**kwargs):
print(key,val)
self._param[key]=val
self.update()
def paintEvent(self, e):
p=self._param
ctr=p['ctr']
aa=p['aa']
qp = QPainter()
qp.begin(self)
qp.translate(200,200)
qp.translate(ctr[0],ctr[1])
qp.setPen(QtGui.QPen(QtCore.Qt.black, 2, QtCore.Qt.SolidLine))
qp.setBrush(QColor(255, 80, 0, 128))
qp.drawEllipse(-100, -100, 200, 200) # big chamber
qp.drawEllipse(-10, -10, 20, 20) # target
qp.drawLine(0,-150,0,-10) #beam
qp.rotate(-45)
rCmb=100 # radius of chamber
rTrg=10 # radius of target
qp.drawEllipse(-rCmb, -rCmb, 2*rCmb, 2*rCmb) # big chamber
qp.drawEllipse(-rTrg, -rTrg, 2*rTrg, 2*rTrg) # target
qp.drawLine(0,-rCmb-50,0,-rTrg) #beam
qp.rotate(-aa)
qp.drawRect(-50, 100, 100, 300) # girder
qp.drawEllipse(-50, 100+150, 20, 20) #pusher left
qp.drawEllipse(-50, 100+150, 20, 20) #pusher left
qp.drawEllipse( 50-20, 100+150, 20, 20) #pusher right
qp.drawEllipse(-50, 100+250, 20, 20) #air pad left
qp.drawEllipse( 50-20, 100+250, 20, 20) #air pad right
qp.end()
@@ -106,10 +130,10 @@ class RIXSgrating(QWidget):
super().__init__()
self._param={'ctr':(400,400), # location of vlsg
'aa':20, #grating angle
'gg':90, #detector angle
'r1':350,#distance probe grating
'r2':450,#distance grating detector
'aa':22, #grating angle
'cc':58, #detector angle
'r1':314,#distance probe grating
'r2':447,#distance grating detector
'szG':(200,5), #size VLS grating
'szD':(150,5), #size detector
}
@@ -159,7 +183,7 @@ class RIXSgrating(QWidget):
l=QtGui.QVBoxLayout(w)
sld={}
for key,rng,tk,pos in (('r1',(50,800),50,60),('r2',(50,800),50,80),('aa',(0,90),5,20),('gg',(0,180),5,40),):
for key,rng,tk,pos in (('r1',(50,800),50,60),('r2',(50,800),50,80),('aa',(0,90),5,20),('cc',(0,180),5,40),):
sl=QtGui.QSlider(QtCore.Qt.Horizontal)
sl.setFixedWidth(200);sl.setMinimum(rng[0]);sl.setMaximum(rng[1])
sl.setValue(self._param[key])
@@ -182,7 +206,7 @@ class RIXSgrating(QWidget):
r2=p['r2']
aa=p['aa']
bb=aa*1.5
gg=p['gg']
cc=p['cc']
szG=p['szG']
szD=p['szD']
ctr=p['ctr']
@@ -221,7 +245,7 @@ class RIXSgrating(QWidget):
tf2=qp.transform()
qp.drawLine(0,0,r2,0) #beam
qp.translate(r2,0)
qp.rotate(-gg)
qp.rotate(-cc)
tf3=qp.transform()
qp.drawRect(int(-szD[0]/2),0 , szD[0], szD[1]) # detector
@@ -235,18 +259,43 @@ class RIXSgrating(QWidget):
qp.drawLine(p0,p1)
# beam vlsg to detector
n=16
pen=QtGui.QPen(QtCore.Qt.black, 3, QtCore.Qt.SolidLine)
col=pen.color()
#n=4;mode=2;alpha=190
#n=4;mode=1;alpha=200
#n=32;mode=1;alpha=120
#n=8;mode=1;alpha=196
n=32;mode=0;alpha=120;width=3
#n=32;mode=2;alpha=255
scl=256/n
p0=QtCore.QPointF(*tf1.map(-szG[0]/2, 0))
p1=QtCore.QPointF(*tf1.map(szG[0]/2, 0))
for i in range(n):
col.setHsv(i*16,255,255,120)
pen.setColor(col)
qp.setPen(pen)
p2=QtCore.QPointF(*tf3.map(-szD[0]/2+szD[0]*i/(n-1),0))
qp.drawLine(p0,p2)
qp.drawLine(p1,p2)
if mode==0:
pen=QtGui.QPen(QtCore.Qt.black, width, QtCore.Qt.SolidLine)
col=QtGui.QColor()
for i in range(n):
p2=QtCore.QPointF(*tf3.map(szD[0]/2-szD[0]*i/(n-1),0))
col.setHsv(i*scl,255,255,alpha)
pen.setColor(col)
qp.setPen(pen)
qp.drawLine(p0,p2)
qp.drawLine(p1,p2)
if mode==1:
col=QtGui.QColor()
qp.setPen(QtGui.QPen(QtCore.Qt.black, 0, QtCore.Qt.SolidLine))
for i in range(n):
p2=QtCore.QPointF(*tf3.map(szD[0]/2-szD[0]*i/(n-1), 0))
col.setHsv(i*scl, 255, 255, alpha)
qp.setBrush(col)
qp.drawPolygon(p0, p1, p2)
if mode==2:
col=QtGui.QColor()
qp.setPen(QtGui.QPen(QtCore.Qt.black, 0, QtCore.Qt.SolidLine))
for i in range(n):
p2=QtCore.QPointF(*tf3.map(szD[0]/2-szD[0]*(i+1)/n, 0))
p3=QtCore.QPointF(*tf3.map(szD[0]/2-szD[0]*i/n, 0))
col.setHsv(i*scl, 255, 255, alpha)
qp.setBrush(col)
qp.drawPolygon(p0, p1, p2, p3)
qp.setCompositionMode(QtGui.QPainter.CompositionMode_SourceOver)
pen=QtGui.QPen(QtCore.Qt.red, 3, QtCore.Qt.SolidLine)
@@ -271,22 +320,18 @@ class RIXSgrating(QWidget):
qp.drawLine(p2,p9)
s=50
qp.drawArc(int(p2.x())-s, int(p2.y())-s, 2*s, 2*s, 180*16, (aa)*16)
qp.drawText(p2+QtCore.QPointF(-s,s),f'aa={aa:.5g}°')
s1=s*np.cos(aa*np.pi/180); s2=s*np.sin(aa*np.pi/180)
qp.drawText(p2+QtCore.QPointF(-s1,s2+20),f'aa={aa:.5g}°')
qp.drawArc(int(p2.x())-s, int(p2.y())-s, 2*s, 2*s, aa*16, bb*16)
qp.drawText(p2+QtCore.QPointF(s,s),f'bb={bb:.5g}°')
qp.drawText(p2+QtCore.QPointF(s1,-s2+20),f'bb={bb:.5g}°')
qp.drawArc(int(p6.x())-s, int(p6.y())-s, 2*s, 2*s, (180+aa+bb)*16, (gg)*16)
qp.drawArc(int(p6.x())-s, int(p6.y())-s, 2*s, 2*s, (180+aa+bb)*16, (cc)*16)
p10=QtCore.QPointF(*tf2.map(r2+20, 0))
qp.drawText(p10,f'gg={gg:.5g}°')
qp.drawText(p10,f'cc={cc:.5g}°')
#for i,p in enumerate((p0,p1,p2,p3,p4,p5,p6,p7,p8,p9)):
# qp.drawText(p,f'P{i}')
#p5=QtCore.QPointF(*tf1.map(-r2-10, 50+w))
qp.end()