eic-opticks/ana/tangential.py

0001 #!/usr/bin/env python
0002 """
0003 tangential.py
0004 ==============
0005
0006 Calculates positions of points on tangent planes to sphere
0007 using a tangential frame.
0008
0009 See also:
0010
0011 * ana/spherical.py
0012 * tangential.cc
0013 * https://mathworld.wolfram.com/SphericalCoordinates.html
0014 * https://en.wikipedia.org/wiki/Vector_fields_in_cylindrical_and_spherical_coordinates
0015
0016 """
0017
0018 import numpy as np
0019 from collections import OrderedDict as odict
0020 from opticks.ana.spherical import Spherical
0021
0022 import pyvista as pv
0023 pv.set_plot_theme("dark")
0024
0025 _white = "ffffff"
0026 _red = "ff0000"
0027 _green = "00ff00"
0028 _blue = "0000ff"
0029
0030 DTYPE = np.float64
0031 FLOAT = [float, np.float64, np.float32]
0032
0033 SIZE = np.array([1280, 720])
0034
0035 def spherical_to_cartesian( rtp ):
0036     """
0037        | x |       |  r*sin(theta)*cos(phi)  |
0038        |   |       |                         |
0039        | y |   =   |  r*sin(theta)*sin(phi)  |
0040        |   |       |                         |
0041        | z |       |      r*cos(theta)       |
0042     """
0043
0044     r,t,p = rtp
0045
0046     if type(t) in FLOAT and type(p) in FLOAT:
0047         num = 1
0048     else:
0049         assert len(t) == len(p)
0050         num = len(t)
0051     pass
0052
0053     radius = r
0054     theta = t*np.pi
0055     phi = p*np.pi
0056
0057     x = radius*np.sin(theta)*np.cos(phi)
0058     y = radius*np.sin(theta)*np.sin(phi)
0059     z = radius*np.cos(theta)
0060     w = DTYPE(1.)
0061
0062     xyzw = np.zeros([num, 4], dtype=DTYPE)
0063     xyzw[:,0] = x.ravel()
0064     xyzw[:,1] = y.ravel()
0065     xyzw[:,2] = z.ravel()
0066     xyzw[:,3] = w.ravel()
0067
0068     return xyzw
0069
0070 def cartesian_to_spherical( xyzw ):
0071     """
0072     """
0073     x,y,z,w = xyzw
0074
0075     num_x = 1 if type(x) in FLOAT else len(x)
0076     num_y = 1 if type(y) in FLOAT else len(y)
0077     num_z = 1 if type(z) in FLOAT else len(z)
0078     num_w = 1 if type(w) in FLOAT else len(w)
0079
0080     num = num_x*num_y*num_z*num_w
0081     assert num_w == 1
0082
0083     r = np.sqrt( x*x + y*y + z*z )
0084     t = np.arccos(z/r)
0085     p = np.arctan(y/x)
0086
0087     rtp = np.zeros( [num,3], dtype=DTYPE )
0088     rtp[:,0] = r
0089     rtp[:,1] = t
0090     rtp[:,2] = p
0091
0092     return rtp
0093
0094
0095 def spherical_to_cartesian_unit_vectors( r, theta, phi ):
0096     """
0097     Cartesian unit vectors (xu yu zu) in terms of spherical unit vectors related by
0098     the inverse of the above transform which is its transpose::
0099
0100
0101           | xu |     |   sin(theta)cos(phi)    cos(theta)cos(phi)      -sin(phi)    |  | ru |
0102           |    |     |                                                              |  |    |
0103           | yu | =   |  sin(theta)sin(phi)    cos(theta)sin(phi)        cos(phi)    |  | tu |
0104           |    |     |                                                              |  |    |
0105           | zu |     |   cos(theta)               -sin(theta)              0        |  | pu |
0106
0107     """
0108     pass
0109
0110 def cartesian_to_spherical_unit_vectors( r, theta, phi ):
0111     """
0112     Spherical unit vectors (ru tu pu) related to cartesian unit vectors (xu yu zu)
0113     via orthogonal rotation matrix::
0114
0115           | ru |     |   sin(theta)cos(phi)    sin(theta)sin(phi)      cos(theta)   |  | xu |
0116           |    |     |                                                              |  |    |
0117           | tu | =   |  cos(theta)cos(phi)    cos(theta)sin(phi)     -sin(theta)    |  | yu |
0118           |    |     |                                                              |  |    |
0119           | pu |     |  -sin(phi)                 cos(phi)              0           |  | zu |
0120
0121     """
0122     c2s = np.zeros( [4,4], dtype=DTYPE )
0123     c2s[0,0] = np.sin(theta)*np.cos(phi)   ; c2s[0,1] = np.sin(theta)*np.sin(phi)   ; c2s[0,2] = np.cos(theta)    ; c2s[0,3] = 0.
0124     c2s[1,0] = np.cos(theta)*np.cos(phi)   ; c2s[1,1] = np.cos(theta)*np.sin(phi)   ; c2s[1,2] = -np.sin(theta)   ; c2s[1,3] = 0.
0125     c2s[2,0] = -np.sin(phi)                ; c2s[2,1] = np.cos(phi)                 ; c2s[2,2] =  0.              ; c2s[2,3] = 0.
0126     c2s[3,0] = 0.                          ; c2s[3,1] = 0.                          ; c2s[3,2] =  0.              ; c2s[3,3] = 1.
0127     return c2s
0128
0129
0130 class Tangential(object):
0131     """
0132     Consider coordinate systems describing a point on a sphere:
0133
0134     1. spherical (r,t,p) with origin at center of sphere
0135     2. cartesian (x,y,z,1) with origin at center of sphere
0136     3. cartesian using
0137
0138        * origin : fixed point (r,t,p) on the surface of the sphere
0139        * unit vectors : normal-to-sphere and a conventional choice of "theta" and "phi" tangents-to-sphere
0140          at the fixed point (r,t,p) [see spherical.py for a demonstration of such frames]
0141
0142     What are the 4x4 matrices to transform frames 2 -> 3 and vice-versa ?
0143     They will be some combination of the below rotation and translation matrices or their inverses::
0144
0145           | ru |     |   sin(theta)cos(phi)    sin(theta)sin(phi)      cos(theta)        0  |  | xu |
0146           |    |     |                                                                      |  |    |
0147           | tu |     |  cos(theta)cos(phi)    cos(theta)sin(phi)     -sin(theta)         0  |  | yu |
0148           |    |     |                                                                      |  |    |
0149           | pu |     |  -sin(phi)                 cos(phi)              0                0  |  | zu |
0150           |    |     |                                                                      |  |    |
0151           |    |     |   0                         0                    0                1  |  |    |
0152
0153      Translation from center to the fixed point
0154
0155           | ru |     |   1                     0                       0                 0  |  | xu |
0156           |    |     |                                                                      |  |    |
0157           | tu |     |   0                     1                       0                 0  |  | yu |
0158           |    |     |                                                                      |  |    |
0159           | pu |     |   0                     0                       1                 0  |  | zu |
0160           |    |     |                                                                      |  |    |
0161           |    |     |   r sin(theta)cos(phi)    r sin(theta)sin(phi)    r cos(theta)    1  |  |    |
0162
0163
0164     """
0165
0166     def __init__(self, rtp):
0167         r,t,p = rtp
0168         theta = t*np.pi
0169         phi = p*np.pi
0170
0171         rot = np.zeros( [4,4], dtype=DTYPE )
0172         iro = np.zeros( [4,4], dtype=DTYPE )
0173         tra = np.zeros( [4,4], dtype=DTYPE )
0174         itr = np.zeros( [4,4], dtype=DTYPE )
0175
0176         rot[0,0] = np.sin(theta)*np.cos(phi)    ; rot[0,1] = np.sin(theta)*np.sin(phi)    ; rot[0,2] = np.cos(theta)     ; rot[0,3] = 0.
0177         rot[1,0] = np.cos(theta)*np.cos(phi)    ; rot[1,1] = np.cos(theta)*np.sin(phi)    ; rot[1,2] = -np.sin(theta)    ; rot[1,3] = 0.
0178         rot[2,0] = -np.sin(phi)                 ; rot[2,1] = np.cos(phi)                  ; rot[2,2] =  0.               ; rot[2,3] = 0.
0179         rot[3,0] = 0.                           ; rot[3,1] = 0.                           ; rot[3,2] =  0.               ; rot[3,3] = 1.
0180         iro = rot.T
0181
0182         tra[0,0] = 1.                           ; tra[0,1] = 0.                           ; tra[0,2] = 0.                ; tra[0,3] = 0.
0183         tra[1,0] = 0.                           ; tra[1,1] = 1.                           ; tra[1,2] = 0.                ; tra[1,3] = 0.
0184         tra[2,0] = 0.                           ; tra[2,1] = 0.                           ; tra[2,2] = 1.                ; tra[2,3] = 0.
0185         tra[3,0] = r*np.sin(theta)*np.cos(phi)  ; tra[3,1] =  r*np.sin(theta)*np.sin(phi) ; tra[3,2] = r*np.cos(theta)   ; tra[3,3] = 1.
0186
0187         itr[0,0] = 1.                           ; itr[0,1] = 0.                           ; itr[0,2] = 0.                ; itr[0,3] = 0.
0188         itr[1,0] = 0.                           ; itr[1,1] = 1.                           ; itr[1,2] = 0.                ; itr[1,3] = 0.
0189         itr[2,0] = 0.                           ; itr[2,1] = 0.                           ; itr[2,2] = 1.                ; itr[2,3] = 0.
0190         itr[3,0] = -r*np.sin(theta)*np.cos(phi) ; itr[3,1] = -r*np.sin(theta)*np.sin(phi) ; itr[3,2] = -r*np.cos(theta)  ; itr[3,3] = 1.
0191
0192         itr_rot = np.dot(itr, rot)
0193         iro_tra = np.dot(iro, tra)
0194         rot_tra = np.dot(rot, tra)
0195
0196         self.rot = rot
0197         self.iro = iro
0198         self.tra = tra
0199         self.itr = itr
0200         self.itr_rot = itr_rot
0201         self.iro_tra = iro_tra
0202         self.rot_tra = rot_tra
0203
0204     def conventional_to_tangential(self, xyzw):
0205         """
0206         :param xyzw: conventional cartesian coordinate in frame with  origin at center of sphere
0207         :return rtpw: tangential cartesian coordinate in frame with origin at point on surface of sphere
0208                       and with tangent-to-sphere and normal-to-sphere as "theta" "phi" unit vectors
0209         """
0210         return np.dot( xyzw, self.itr_rot )
0211
0212     def tangential_to_conventional(self, rtpw):
0213         """
0214
0215                    z:p
0216                    |
0217                    |  y:t
0218                    | /
0219                    |/
0220                    +-----x:r
0221
0222
0223         :param rtpw: tangential cartesian coordinate in frame with origin at point on surface of sphere
0224                       and with tangent-to-sphere and normal-to-sphere as "theta" "phi" unit vectors
0225         :return  xyzw: conventional cartesian coordinate in frame with  origin at center of sphere
0226         """
0227         #return np.dot( rtpw, self.iro_tra )
0228         return np.dot( rtpw, self.rot_tra )
0229
0230     def get_plane_rtpw(self, side=10, num=20 ):
0231         v_t = np.linspace( -side, side, num )
0232         v_p = np.linspace( -side, side, num )
0233         t, p = np.meshgrid(v_t, v_p)
0234         r = np.zeros_like(t)
0235         w = np.ones_like(t)
0236         _rtpw = np.dstack( (r, t, p, w))
0237         return _rtpw
0238
0239     def get_plane_xyzw(self, side=10, num=20):
0240         _rtpw = self.get_plane_rtpw(side=side, num=num )
0241         _xyzw = self.tangential_to_conventional( _rtpw )
0242         return _xyzw
0243
0244     def __repr__(self):
0245         return "\n".join(map(str, ["itr", self.itr, "rot", self.rot, "tra", self.tra, "itr_rot", self.itr_rot, "iro_tra", self.iro_tra ]))
0246
0247     def pvplot(self, pl, color=_white, side=10, num=20 ):
0248         _xyzw = self.get_plane_xyzw(side=side, num=num)
0249         pos =  _xyzw[:,:,:3].reshape(-1,3)
0250         pl.add_points( pos,  color=color )
0251
0252
0253 if __name__ == '__main__':
0254
0255     rtp = odict()
0256     xyzw = odict()
0257     xyzw2 = odict()
0258     rtpw = odict()
0259     ta = odict()
0260
0261     radius = 20.
0262
0263     #rtp["center"]           = np.array( [0., 0., 0.],       dtype=DTYPE )
0264     rtp["north_pole"]       = np.array( [radius, 0,   0.],  dtype=DTYPE )
0265     rtp["null_island"]      = np.array( [radius, 0.5, 0.],  dtype=DTYPE )
0266     rtp["mid_island"]       = np.array( [radius, 0.5, 0.5], dtype=DTYPE )
0267     rtp["anti-null_island"] = np.array( [radius, 0.5, 1.],  dtype=DTYPE )
0268     rtp["south_pole"]       = np.array( [radius, 1,   0.],  dtype=DTYPE )
0269     rtp["midlat"]           = np.array( [radius, 0.25, 0.], dtype=DTYPE )
0270     rtp["midlat2"]          = np.array( [radius, 0.25, 0.5], dtype=DTYPE )
0271
0272     for k in rtp: xyzw[k] = spherical_to_cartesian( rtp[k] )
0273     for k in rtp: ta[k] = Tangential( rtp[k] )
0274     for k in rtp: rtpw[k] = ta[k].conventional_to_tangential(xyzw[k])
0275     for k in rtp: xyzw2[k] = ta[k].tangential_to_conventional(rtpw[k])
0276
0277     fmt = "%20s %30s %30s %30s %30s"
0278     print(fmt % ("name", "rtp", "xyzw", "rtpw", "xyzw2"))
0279
0280     for k in rtp: print(fmt % (k, rtp[k], xyzw[k], rtpw[k], xyzw2[k] ))
0281
0282     all_rtp = np.vstack(rtp.values())
0283     all_xyzw = spherical_to_cartesian( all_rtp.T )
0284
0285     sg = Spherical.Grid(radius=radius, n_theta=24, n_phi=24 )
0286
0287
0288     sp = pv.Sphere(radius=radius)
0289
0290
0291     pl = pv.Plotter(window_size=SIZE*2 )
0292     pl.show_grid()
0293
0294     sg.pvplot(pl, mag=2.5)
0295
0296     pl.add_mesh(sp)
0297
0298     #ta["midlat2"].pvplot(pl, color="red", side=10, num=20 )
0299
0300     #ta["north_pole"].pvplot(pl, color="blue" )
0301     #ta["south_pole"].pvplot(pl, color="yellow" )
0302     #ta["null_island"].pvplot(pl, color="cyan" )
0303     #ta["anti-null_island"].pvplot(pl, color="magenta" )
0304
0305
0306     look = np.array( [0,0,0])
0307     up = np.array( [0,0,1])
0308     eye =  3.0*radius*np.array([1.,1.,1.])/np.sqrt(3.)
0309
0310     pl.set_focus(    look )
0311     pl.set_viewup(   up )
0312     pl.set_position( eye, reset=True )
0313
0314     pl.camera.Zoom(2.0)
0315
0316     outpath = "/tmp/tangential.png"
0317     cp = pl.show(screenshot=outpath)
0318
0319
0320