#!/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 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), (200, 1650, 3500, np.deg2rad(88)), ] } else: raise ValueError('unknown type') 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)) return geo def geometry2raw(self, geo): r1,r2,alpha,beta,gamma=geo # returns raw motor positions mt=gtz=gty1=gty2=grx=gtx=dtz=dty1=dty2=drx=None aa=np.deg2rad(90-alpha) bb=np.deg2rad(90-beta) gg=np.deg2rad(90-beta) gtz=r1 grx=90-alpha dtz=np.cos(aa+bb)*r2 dty1=dty2=np.sin(aa+bb)*r2 drx=gamma raw=(mt,gtz,gty1,gty2,grx,gtx,dtz,dty1,dty2,drx) return raw # ---------- 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))