eic-opticks/ana/spherical_0.py

0001 #!/usr/bin/env python
0002 """
0003 spherical_0.py
0004 ================
0005
0006 See spherical.py and tangential.py for further development on this theme.
0007
0008 https://mathworld.wolfram.com/SphericalCoordinates.html
0009
0010 Draws points on a sphere together with a (r,phi,theta) derivative vectors
0011 which are normals and tangents to the sphere and form orthogonal
0012 frames at every point on the sphere.
0013
0014                          .
0015
0016                             .
0017
0018                               .     pvec (blue)   increasing phi : eg around equator or other lines of latitude, from West to East
0019                                    /
0020                                .  /
0021                                  /
0022       +                         +-------->
0023     center                      |         rvec (red)   radial normal vector
0024                                .|
0025                                 |
0026                              . tvec (green) increasing theta : eg down lines of longitude : from North to South
0027
0028                             .
0029
0030                          ,
0031
0032 """
0033 import numpy as np
0034 import pyvista as pv
0035 SIZE = np.array([1280, 720])
0036 DTYPE = np.float64
0037
0038 if __name__ == '__main__':
0039
0040      radius = 20.
0041      v_phi = np.pi*np.linspace(0,2,50, dtype=DTYPE)
0042      v_theta = np.pi*np.linspace(0,1,50, dtype=DTYPE )
0043
0044      num = len(v_theta)*len(v_phi)
0045
0046      theta, phi = np.meshgrid(v_theta, v_phi)
0047      print("theta.shape %s " % str(theta.shape))
0048      print("phi.shape %s " % str(phi.shape))
0049
0050      rvec = np.zeros( (num,3),  dtype=DTYPE )
0051      pvec = np.zeros( (num,3),  dtype=DTYPE )
0052      tvec = np.zeros( (num,3),  dtype=DTYPE )
0053
0054      # derivative with radius of definition of spherical coordinates : direction of increasing radius : normal vector
0055      x = np.cos(phi)*np.sin(theta)
0056      y = np.sin(phi)*np.sin(theta)
0057      z = np.cos(theta)
0058
0059      rvec[:,0] = x.ravel()
0060      rvec[:,1] = y.ravel()
0061      rvec[:,2] = z.ravel()
0062      pos = rvec*radius
0063
0064      r = np.sqrt(x*x + y*y + z*z)
0065
0066      # derivative with phi : direction of increasing phi : phi tangent vector
0067      pvec[:,0] = -np.sin(phi).ravel()
0068      pvec[:,1] =  np.cos(phi).ravel()
0069      pvec[:,2] = 0
0070
0071      # derivative with theta  : direction of increasing theta : theta tangent vector
0072      tvec[:,0] =  (np.cos(phi)*np.cos(theta)).ravel()
0073      tvec[:,1] =  (np.sin(phi)*np.cos(theta)).ravel()
0074      tvec[:,2] = -np.sin(theta).ravel()
0075
0076      # check orthogonality of the frames at every point
0077      check_pvec_tvec = np.abs(np.sum( pvec*tvec, axis=1 ))
0078      check_pvec_rvec = np.abs(np.sum( pvec*rvec, axis=1 ))
0079      check_tvec_rvec = np.abs(np.sum( tvec*rvec, axis=1 ))
0080
0081      assert check_pvec_tvec.max() < 1e-6
0082      assert check_pvec_rvec.max() < 1e-6
0083      assert check_tvec_rvec.max() < 1e-6
0084
0085      pl = pv.Plotter(window_size=SIZE*2 )
0086
0087      _white = "ffffff"
0088      _red = "ff0000"
0089      _green = "00ff00"
0090      _blue = "0000ff"
0091
0092      pl.add_points( pos,       color=_white )
0093      pl.add_arrows( pos, rvec, color=_red )
0094      pl.add_arrows( pos, tvec, color=_green )
0095      pl.add_arrows( pos, pvec, color=_blue )
0096
0097      cp = pl.show()
0098
0099