eic-opticks/ana/spherical.py

0001 #!/usr/bin/env python
0002 """
0003 spherical.py
0004 =============
0005
0006 Draws points on a sphere together with a (r,phi,theta) derivative vectors
0007 which are normals and tangents to the sphere and form orthogonal
0008 frames at every point on the sphere (other than perhaps the poles).
0009
0010                          .
0011
0012                             .
0013
0014                               .     pvec (blue)   increasing phi : eg around equator or other lines of latitude, from West to East
0015                                    /
0016                                .  /
0017                                  /
0018       +                         +-------->
0019     center                      |         rvec (red)   radial normal vector
0020                                .|
0021                                 |
0022                              . tvec (green) increasing theta : eg down lines of longitude : from North to South
0023
0024                             .
0025
0026                          ,
0027
0028 See also:
0029
0030 * ana/spherical_0.py ana/tangential.py
0031 * https://mathworld.wolfram.com/SphericalCoordinates.html
0032 * https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates
0033
0034 """
0035
0036 import numpy as np
0037 import pyvista as pv
0038
0039 _white = "ffffff"
0040 _red = "ff0000"
0041 _green = "00ff00"
0042 _blue = "0000ff"
0043
0044 DTYPE = np.float64
0045 SIZE = np.array([1280, 720])
0046
0047 class Spherical(object):
0048     @classmethod
0049     def Grid(cls, radius=20., n_theta=50, n_phi=50):
0050         v_theta = np.linspace(0,1,n_theta, dtype=DTYPE )
0051         v_phi = np.linspace(0,2,n_phi, dtype=DTYPE)
0052         num = len(v_theta)*len(v_phi)
0053         theta, phi = np.meshgrid(v_theta, v_phi)
0054         print("v_theta.shape %s theta.shape %s " % (str(v_theta.shape),str(theta.shape)) )
0055         print("v_phi.shape %s phi.shape %s " % (str(v_phi.shape),str(phi.shape)))
0056         sph = cls((radius, theta, phi))
0057         return sph
0058
0059     @classmethod
0060     def One(cls, radius=20., theta=0.5, phi=0.0):
0061         sph = cls((radius, theta, phi))
0062         return sph
0063
0064     def __init__(self, rtp ):
0065         r,t,p = rtp
0066
0067         num_r = 1 if type(r) is float else len(r)
0068         num_t = 1 if type(t) is float else len(t)
0069         num_p = 1 if type(p) is float else len(p)
0070         num = num_r*num_t*num_p
0071
0072         assert num_r == 1
0073
0074         radius = r
0075         theta = t*np.pi
0076         phi = p*np.pi
0077
0078         # derivative with radius of definition of spherical coordinates : direction of increasing radius : normal vector
0079         x = np.cos(phi)*np.sin(theta)
0080         y = np.sin(phi)*np.sin(theta)
0081         z = np.cos(theta)
0082
0083         one = np.sqrt(x*x + y*y + z*z)
0084
0085         rvec = np.zeros( (num,4),  dtype=DTYPE )
0086         pvec = np.zeros( (num,4),  dtype=DTYPE )
0087         tvec = np.zeros( (num,4),  dtype=DTYPE )
0088
0089         rvec[:,0] = x.ravel()
0090         rvec[:,1] = y.ravel()
0091         rvec[:,2] = z.ravel()
0092         rvec[:,3] = 0.
0093
0094         # derivative with theta  : direction of increasing theta : theta tangent vector
0095         tvec[:,0] =  (np.cos(phi)*np.cos(theta)).ravel()
0096         tvec[:,1] =  (np.sin(phi)*np.cos(theta)).ravel()
0097         tvec[:,2] = -np.sin(theta).ravel()
0098         tvec[:,3] = 0.
0099
0100         # derivative with phi : direction of increasing phi : phi tangent vector
0101         pvec[:,0] = -np.sin(phi).ravel()
0102         pvec[:,1] =  np.cos(phi).ravel()
0103         pvec[:,2] = 0.
0104         pvec[:,3] = 0.
0105
0106         xyzw = rvec*radius
0107
0108         # check orthogonality of the frames at every point
0109         check_pvec_tvec = np.abs(np.sum( pvec*tvec, axis=1 ))
0110         check_pvec_rvec = np.abs(np.sum( pvec*rvec, axis=1 ))
0111         check_tvec_rvec = np.abs(np.sum( tvec*rvec, axis=1 ))
0112
0113         assert check_pvec_tvec.max() < 1e-6
0114         assert check_pvec_rvec.max() < 1e-6
0115         assert check_tvec_rvec.max() < 1e-6
0116
0117         self.rtp = rtp
0118         self.xyzw = xyzw
0119         self.rvec = rvec
0120         self.pvec = pvec
0121         self.tvec = tvec
0122
0123     def __repr__(self):
0124         return "rtp %s rvec %s tvec %s pvec %s " % (str(self.rtp), str(self.rvec), str(self.tvec), str(self.pvec))
0125
0126     def pvplot(self, pl, mag=1):
0127         s = self
0128         pl.add_points( s.xyzw[:,:3], color=_white )
0129         pl.add_arrows( s.xyzw[:,:3], s.rvec[:,:3], mag=mag, color=_red )
0130         pl.add_arrows( s.xyzw[:,:3], s.tvec[:,:3], mag=mag, color=_green )
0131         pl.add_arrows( s.xyzw[:,:3], s.pvec[:,:3], mag=mag, color=_blue )
0132
0133
0134 if __name__ == '__main__':
0135
0136     sg = Spherical.Grid()
0137     s1 = Spherical.One(radius=22.)
0138
0139     pl = pv.Plotter(window_size=SIZE*2 )
0140
0141     sg.pvplot(pl)
0142     s1.pvplot(pl)
0143
0144     cp = pl.show()
0145
0146
0147