1 """
2 Vector and coordinate identities and transformations
3
4 General conventions:
5
6 - theta is the colatitude, i.e. the angle from the Z-axis
7 - phi is the longitude
8
9
10 Transform cartesian coordinates to spherical coordinates and back
11
12 >>> x,y,z = np.random.uniform(low=-1,high=1,size=(3,10))
13 >>> r,phi,theta = cart2spher_coord(x,y,z)
14 >>> x_,y_,z_ = spher2cart_coord(r,phi,theta)
15 >>> np.allclose(x,x_),np.allclose(y,y_),np.allclose(z,z_)
16 (True, True, True)
17
18 Transform cartesian vectors to spherical vectors and back
19
20 >>> x0,y0,z0 = np.random.uniform(low=-1,high=1,size=(3,10))
21 >>> x1,y1,z1 = np.random.uniform(low=-1,high=1,size=(3,10))
22 >>> r0,phi0,theta0 = cart2spher_coord(x0,y0,z0)
23 >>> r1,phi1,theta1 = cart2spher((x0,y0,z0),(x1,y1,z1))
24 >>> x1_,y1_,z1_ = spher2cart((r0,phi0,theta0),(r1,phi1,theta1))
25 >>> np.allclose(x1,x1_),np.allclose(y1,y1_),np.allclose(z1,z1_)
26 (True, True, True)
27
28 """
29
30 import numpy as np
31 from numpy import cos,sin,pi,sqrt
32
33
34
35
37 """
38 Spherical to Cartesian coordinates
39 """
40 x = r*cos(phi)*sin(theta)
41 y = r*sin(phi)*sin(theta)
42 z = r*cos(theta)
43 return x,y,z
44
46 """
47 theta colatitude
48 phi longitude
49 """
50 rho = np.sqrt(x**2+y**2+z**2)
51 phi = np.arctan2(y,x)
52 theta = np.arctan2(np.sqrt(x**2+y**2),z)
53 return rho,phi,theta
54
55 -def rotate(x,y,theta,x0=0.,y0=0.):
56 """
57 Rotate counter clockwise.
58 """
59 xnew = cos(theta)*x - sin(theta)*y - x0
60 ynew = sin(theta)*x + cos(theta)*y - y0
61 return xnew,ynew
62
63
64
65
67 """
68 theta is angle from z-axis (colatitude)
69 phi is longitude
70
71 E.g. http://www.engin.brown.edu/courses/en3/Notes/Vector_Web2/Vectors6a/Vectors6a.htm
72
73 >>> np.random.seed(111)
74 >>> r,phi,theta = np.random.uniform(low=-1,high=1,size=(3,2))
75 >>> a_r,a_phi,a_theta = np.random.uniform(low=-1,high=1,size=(3,2))
76 >>> a_x,a_y,a_z = spher2cart((r,phi,theta),(a_r,a_phi,a_theta))
77 """
78 ax = sin(theta)*cos(phi)*a_r + cos(theta)*cos(phi)*a_theta - sin(phi)*a_phi
79 ay = sin(theta)*sin(phi)*a_r + cos(theta)*sin(phi)*a_theta + cos(phi)*a_phi
80 az = cos(theta) *a_r - sin(theta) *a_theta
81 return ax,ay,az
82
83
85 """
86 theta is angle from z-axis (colatitude)
87 phi is longitude
88
89 return r,phi,theta
90 """
91 r,phi,theta = cart2spher_coord(x0,y0,z0)
92 ar = sin(theta)*cos(phi)*x1 + sin(theta)*sin(phi)*y1 + cos(theta)*z1
93 atheta = cos(theta)*cos(phi)*x1 + cos(theta)*sin(phi)*y1 - sin(theta) *z1
94 aphi = -sin(phi) *x1 + cos(phi) *y1
95 return ar,aphi,atheta
96
97
98
99
100
102 """
103 Euclidic norm of a vector (or of a grid of vectors)
104
105 Input vectors should be numpy arrays.
106 """
107 return sqrt((vec**2).sum(axis=0))
108
109
111 """
112 Compute angle between two vectors (or between two grids of vectors).
113
114 Input vectors should be numpy arrays.
115 """
116 return np.arccos( (vec1*vec2).sum(axis=0) / (norm(vec1)*norm(vec2)))
117
118
120 """
121 Compute cosine of angle between two vectors (or between two grids of vectors).
122
123 Input vectors should be numpy arrays.
124 """
125 return (vec1*vec2).sum(axis=0) / (norm(vec1)*norm(vec2))
126
127
128
129
130 if __name__=="__main__":
131 import doctest
132 import pylab as pl
133 doctest.testmod()
134