diff --git a/geometry.py b/geometry.py new file mode 100644 index 0000000..b1435f4 --- /dev/null +++ b/geometry.py @@ -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)) diff --git a/graphExample.py b/graphExample.py index 4acd561..5b2089e 100755 --- a/graphExample.py +++ b/graphExample.py @@ -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()