reflectometry.reduction.qxqz

1 # This code is public domain 2 # Author: Paul Kienzle 3 """ 4 Translate coordinates between real and reciprical space. 5 """ 6 import numpy 7 from numpy import sin, cos, pi, arcsin, sqrt, arctan2, degrees, radians 8 9

10 -def ABL_to_QxQz(sample_angle, detector_angle, wavelength):

11 """ 12 Compute Qx,Qz given incident and reflected angles and wavelength 13 14 Returns (Qx,Qz) 15 16 If all inputs are scalars the result will be a scalar. 17 """ 18 A,B = sample_angle, detector_angle 19 Qz = 2*pi/wavelength * ( sin(radians(B - A)) + sin(radians(A))) 20 Qx = 2*pi/wavelength * ( cos(radians(B - A)) - cos(radians(A))) 21 return Qx,Qz

22

23 -def QxQzL_to_AB(Qx, Qz, wavelength):

24 """ 25 Guess incident and reflected angles given Qx, Qz and wavelength. 26 27 This transform is not invertible: for any Qx, Qz there are two 28 separate choices for sample angle and detector angle which 29 give the same Q. We choose the one for which detector angle 30 matches the sign of Qz. 31 32 Returns sample angle, detector angle 33 34 Sample angle is also called the incident angle, theta_i. 35 Dectector angle is incident angle plus reflected angle, 36 or theta_i + theta_f. 37 38 If all inputs are scalars the result will be a scalar. 39 40 Error is in the order of 1e-12 for off-specular coordinates. 41 """ 42 # Algorithm for converting Qx-Qz-lambda to alpha-beta: 43 # beta = 2 asin(L/(2 pi) sqrt(Qx^2+Qz^2)/2) * 180/pi 44 # = asin(L/(4 pi) sqrt(Qx^2+Qz^2)) * 360/pi 45 # if Qz < 0, negate beta 46 # theta = atan2(Qx,Qz) * 180/pi 47 # if theta > 90, theta -= 360 48 # alpha = theta + beta/2 49 # if Qz < 0, alpha += 180 50 beta = 2*degrees(arcsin(wavelength/(4*pi) * sqrt(Qx**2+Qz**2))) 51 beta *= -2*(Qz<0)+1 # -2+1=-1, 0+1=1, so -2*(cond)+1 = +/-1 52 theta = degrees(arctan2(Qx,Qz)) 53 theta -= 360*(theta>90) 54 alpha = theta + beta/2 55 alpha += 180*(Qz<0) 56 return alpha,beta

57

58 -def QxQzA_to_BL(Qx,Qz,alpha):

59 """ 60 Guess detector angle and wavelength given Qx, Qz and sample angle. 61 62 The detector angle is branch cut between -180 and 180 degrees. 63 64 Returns detector angle, wavelength 65 66 Sample angle is also called the incident angle, theta_i. 67 Dectector angle is incident angle plus reflected angle, 68 or theta_i + theta_f. 69 70 If all inputs are scalars the result will be a scalar. 71 72 Error is in the order of 1e-12 for off-specular. 73 """ 74 theta = degrees(arctan2(Qx,Qz)) 75 theta -= 360*(theta>90) 76 beta = 2*(alpha-theta) 77 beta += 360*(beta<-180) 78 beta -= 360*(beta>=180) 79 beta -= 360*(beta>=180) 80 wavelength = 4*pi*sin(radians(abs(beta))/2) / sqrt(Qx**2 +Qz**2) 81 return beta,wavelength

82

83 -def _errchk(err,tol=1e-15):

84 chk = (numpy.abs(err) < tol).all() 85 if not chk: 86 print err 87 print "norm",numpy.linalg.norm(err) 88 return chk

89

90 -def _test1(A,B,L,X,Z):

91 msg=", ".join(["%g"%numpy.asarray(v).flatten()[0] for v in A,B,L,X,Z]) 92 x,z = ABL_to_QxQz(A,B,L) 93 mx,mz=["%g"%numpy.asarray(v).flatten()[0] for v in x,z] 94 assert _errchk(z-Z),"%s -> %s"%(msg,mz) 95 assert _errchk(x-X),"%s -> %s"%(msg,mx) 96 97 a,b2 = QxQzL_to_AB(X,Z,L) 98 ma,mb2=["%g"%numpy.asarray(v).flatten()[0] for v in a,b2] 99 100 if numpy.any(Z<0): 101 # Full test; make sure the forward transform from the 102 # inverted transform is correct even if it doesn't 103 # happen to match the choice of detector angle. 104 x,z=ABL_to_QxQz(a,b2,L) 105 assert _errchk(x-X,1e-12),"incorrect inverse Qx" 106 assert _errchk(z-Z,1e-12),"incorrect inverse Qz" 107 assert _errchk(numpy.sign(Z)-numpy.sign(b2),1.5),"incorrect branch" 108 idx = (numpy.abs(a-A)<1e-12)|(numpy.abs(b2-B)<1e-12) 109 if False: 110 import pylab 111 #print 1*idx 112 Bs = B*(1+a-a) 113 As = A*(1+a-a) 114 pylab.quiver(Bs[~idx],As[~idx],(b2-B)[~idx],(a-A)[~idx]) 115 pylab.plot(Bs[idx],As[idx],'o') 116 pylab.show() 117 assert _errchk((A-a)[idx],1e-12),"incorrect sample angle" 118 assert _errchk((B-b2)[idx],1e-12),"incorrect detector angle" 119 else: 120 assert _errchk(b2-B,1e-12),"%s -> %s"%(msg,mb2) 121 assert _errchk(a-A,1e-12),"%s -> %s"%(msg,ma) 122 123 b,l = QxQzA_to_BL(X,Z,A) 124 mb,ml=["%g"%numpy.asarray(v).flatten()[0] for v in b,l] 125 if False and not _errchk(b-B, 1e-12): 126 # Debug problems with BL; should be less since results are 127 # not ambiguous. 128 idx = abs(b-B)>1e-12 129 Bs = B*(1+b-b) 130 Ls = L*(1+b-b) 131 import pylab 132 pylab.quiver(Bs[idx],Ls[idx],(b-B)[idx],(l-L)[idx]) 133 pylab.plot(Bs[~idx],Ls[~idx],'o') 134 pylab.plot(Bs[idx],Ls[idx],'o') 135 pylab.plot(b[idx],l[idx],'x') 136 pylab.show() 137 assert _errchk(b-B, 1e-12),"%s -> %s"%(msg,mb) 138 assert _errchk(l-L, 1e-12),"%s -> %s"%(msg,ml)

139 140 141

142 -def test():

143 A,B,L = 3,6,4.5 144 X,Z = 0,4*pi/L*sin(radians(A)) 145 vec = numpy.ones((2,3,4)) 146 147 # Check combos of scalar/vector for Qx=0 148 _test1(A,B,L,X,Z) 149 _test1(vec*A,B,L,X,Z) 150 _test1(A,vec*B,L,X,Z) 151 _test1(A,B,vec*L,X,Z) 152 _test1(A,B,L,vec*X,Z) 153 _test1(A,B,L,X,vec*Z) 154 _test1(vec*A,vec*B,vec*L,vec*X,vec*Z) 155 _test1(vec[:,0:1,0:1]*A,vec[0:1,:,0:1]*B,vec[0:1,0:1,:]*L,X,Z) 156 157 # Check combos of scalar/vector for Qx!=0 158 A,B,L = 3,2,4.5 159 ti,tf = radians(A),radians(B-A) 160 X,Z = 2*pi/L*(cos(tf)-cos(ti)),2*pi/L*(sin(tf)+sin(ti)) 161 _test1(A,B,L,X,Z) 162 _test1(A,B,vec*L,X,Z) 163 _test1(vec*A,B,L,X,Z) 164 _test1(A,vec*B,L,X,Z) 165 _test1(A,B,vec*L,X,Z) 166 _test1(A,B,L,vec*X,Z) 167 _test1(A,B,L,X,vec*Z) 168 _test1(vec*A,vec*B,vec*L,vec*X,vec*Z) 169 _test1(vec[:,0:1,0:1]*A,vec[0:1,:,0:1]*B,vec[0:1,0:1,:]*L,X,Z) 170 171 # Check the whole coordinate space, avoiding Qz=0. 172 A = numpy.linspace(-170,170,5).reshape((1,1,5)) 173 B = numpy.linspace(-170,170,6).reshape((1,6,1)) 174 L = numpy.linspace(0.1,7.1,4).reshape((4,1,1)) 175 #L=4 176 X,Z = ABL_to_QxQz(A,B,L) 177 _test1(A,B,L,X,Z)

178 179 if __name__ == "__main__": 180 test() 181

Source Code for Module reflectometry.reduction.qxqz