eic-opticks/ana/mm0prim2.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 Hmm need to make connection to the volume traversal index
0023 """
0024 import os, logging, numpy as np
0025 log = logging.getLogger(__name__)
0026
0027 from opticks.ana.key import keydir
0028 #from opticks.ana.prim import Solid
0029 from opticks.ana.geom2d import Geom2d
0030 from opticks.ana.flight import Flight
0031
0032 np.set_printoptions(suppress=True)
0033 dtype = np.float32
0034
0035 try:
0036     import matplotlib.pyplot as plt
0037     from mpl_toolkits.mplot3d import Axes3D
0038     import mpl_toolkits.mplot3d.art3d as art3d
0039 except ImportError:
0040     plt = None
0041 pass
0042
0043
0044 if __name__ == '__main__':
0045
0046     #fmt = '[%(asctime)s] p%(process)s {%(pathname)s:%(lineno)d} %(levelname)s - %(message)s'
0047     fmt = '{%(pathname)s:%(lineno)d} - %(message)s'
0048     logging.basicConfig(level=logging.INFO, format=fmt)
0049
0050     ok = os.environ["OPTICKS_KEY"]
0051     kd = keydir(ok)
0052     log.info(kd)
0053     assert os.path.exists(kd), kd
0054
0055     os.environ["IDPATH"] = kd    ## TODO: avoid having to do this, due to prim internals
0056
0057     mm0 = Geom2d(kd, ridx=0)     ## TODO: renamings/splittings more suited to new all-in-GNodeLib-not-mm0 geometry model
0058
0059     target = mm0.pvfind('pTarget')    ## pick a convenient frame of operation
0060     assert len(target) == 1
0061     target = target[0]
0062
0063     sc = mm0.ce[target][3]/1000.      ##  big radius in meters 17.760008
0064
0065     if plt:
0066         plt.ion()
0067         fig = plt.figure(figsize=(9,9))
0068         ax = fig.add_subplot(111,projection='3d')
0069         plt.title("mm0 geom2d")
0070         sz = 25
0071
0072         ax.set_xlim([-sz,sz])
0073         ax.set_ylim([-sz,sz])
0074         ax.set_zlim([-sz,sz])
0075
0076         ax.set_xlabel("x")
0077         ax.set_ylabel("y")
0078         ax.set_zlabel("z")
0079
0080         mm0.render(ax, art3d=art3d)
0081     pass
0082
0083
0084     cmd = False
0085     #################### x-yoyo at (y=0,z=0)
0086     n9 = 7
0087     ezz = np.zeros( (n9,4), dtype=dtype )
0088     lzz = np.zeros( (n9,4), dtype=dtype )
0089     uzz = np.zeros( (n9,4), dtype=dtype )
0090     czz = np.zeros( (n9,4), dtype=dtype )
0091
0092     ezz[:,0] = [-3,-2,-1,0,-1,-2,-3]
0093     ezz[:,1] = np.zeros( n9, dtype=dtype )
0094     ezz[:,2] = np.zeros( n9, dtype=dtype )
0095     ezz[:,3] = np.ones(  n9, dtype=dtype )
0096
0097     lzz[:n9] = [1,0,0,0]
0098     uzz[:n9] = [0,0,1,0]
0099
0100     if cmd:
0101         czz[0].view("|S2")[0:8] = ["AA", "C0", "T5", "B2", "Q1", "", "N2", "" ]
0102         czz[2].view("|S2")[0:8] = ["","","","", "A1","","","" ]
0103         czz[n9//2].view("|S2")[0:8] = ["T6", "", "", "", "", "", "", "" ]
0104         # AA: home the animation as ascii A - 0 is greater than the number of modes
0105         # Q1:invis GlobalStyle
0106     pass
0107
0108     fzz = Flight(ezz,lzz,uzz,czz)
0109
0110     ####################  x-scan at (y=0,z=1)
0111     n0 = 3
0112     eaa = np.zeros( (n0,4), dtype=dtype )
0113     laa = np.zeros( (n0,4), dtype=dtype )
0114     uaa = np.zeros( (n0,4), dtype=dtype )
0115     caa = np.zeros( (n0,4), dtype=dtype )
0116
0117     eaa[:,0] = np.linspace(-n0, -1, n0, dtype=dtype )
0118     eaa[:,1] = np.zeros( n0, dtype=dtype )
0119     eaa[:,2] = np.ones( n0, dtype=dtype )
0120     eaa[:,3] = np.ones( n0, dtype=dtype )
0121
0122     uaa[:n0] = [0,0,1,0]
0123
0124     if cmd:
0125         caa[0].view("|S2")[0:8] = ["C1", "T6", "B2", "Q0", "A1", "P1", "N1", "" ]
0126         caa[2].view("|S2")[0:8] = ["Q1","","","Q0","C0","T5","",""]
0127     pass
0128
0129     laa[:-1] = eaa[1:]
0130     laa[-1] = [0,0,1,0]
0131
0132     faa = Flight(eaa,laa,uaa,caa)
0133
0134     #################  x-z circle
0135
0136     pz = 1.0
0137     pr = 1.05
0138
0139     phase0 = np.arccos(pz)
0140     ta = np.linspace( 0, 2*np.pi, 32 )[:-1]
0141     za = np.cos(ta+phase0)
0142     m = np.argmin(np.abs(za[2:]-pz))+2   # index of za closest to that pz value going around again, excluding 0
0143     t0 = ta[:m+1]
0144     n1 = len(t0)
0145     st0 = np.sin(t0+phase0)
0146     ct0 = np.cos(t0+phase0)
0147
0148     ebb = np.zeros( (n1,4) , dtype=dtype )  # eye position
0149     lbb = np.zeros( (n1,4) , dtype=dtype ) # up direction
0150     ubb = np.zeros( (n1,4) , dtype=dtype ) # up direction
0151     cbb = np.zeros( (n1,4),  dtype=dtype )  # ctrl
0152
0153     ebb[:,0] = st0
0154     ebb[:,1] = 0
0155     ebb[:,2] = ct0
0156     ebb *= pr
0157     ebb[:,3] = np.ones( n1, dtype=dtype )
0158
0159     ubb[:,0] = st0
0160     ubb[:,1] = 0
0161     ubb[:,2] = ct0
0162
0163     if cmd:
0164         cbb[0].view("|S2")[0:8] = ["C0", "","","","T6","","T7","" ]
0165         cbb[1].view("|S2")[0:1] = ["T8"]
0166         cbb[n1//2+0].view("|S2")[0:8] = ["E3","","","","","","","" ]
0167         cbb[n1//2+1].view("|S2")[0:8] = ["","","","","","","","" ]
0168         cbb[n1//2+2].view("|S2")[0:8] = ["","","","","","","","" ]
0169         cbb[n1//2+3].view("|S2")[0:8] = ["","","","","","","","" ]
0170         cbb[n1//2+4].view("|S2")[0:8] = ["E1","","","","","","","" ]
0171
0172         cbb[n1-2].view("|S2")[0:1] = ["T5"]
0173     pass
0174
0175     lbb[:-1] = ebb[1:]
0176     lbb[-1] = [0,0,1,1]
0177
0178     fbb = Flight(ebb,lbb,ubb,cbb)
0179
0180     ##############################
0181
0182     # take the last point x value (close to pz) and make xy loop
0183     r2 = np.abs(ebb[-1,0])
0184     tb = np.linspace( 0, 2*np.pi, 16)[:-1]
0185     n2 = len(tb)
0186
0187     ecc = np.zeros( (n2,4), dtype=dtype )
0188     lcc = np.zeros( (n2,4), dtype=dtype )
0189     ucc = np.zeros( (n2,4), dtype=dtype )
0190     ccc = np.zeros( (n2,4), dtype=dtype )
0191
0192     ecc[:,0] = r2*np.cos(tb)
0193     ecc[:,1] = r2*np.sin(tb)
0194     ecc[:,2] = ebb[-1,2]
0195     ecc[:,3] = np.ones( n2, dtype=dtype )
0196
0197     ucc[:,0] = np.zeros(n2, dtype=dtype)
0198     ucc[:,1] = np.zeros(n2, dtype=dtype)
0199     ucc[:,2] = np.ones(n2, dtype=dtype)
0200     ucc[:,3] = np.zeros(n2, dtype=dtype)
0201
0202     if cmd:
0203         ccc[0].view("|S2")[0:8] = ["T7","","","" ,"","","",""]
0204         ccc[n2//2].view("|S2")[0:8] = ["","","","" ,"","","",""]
0205         ccc[n2//2+2].view("|S2")[0:8] = ["","","","" ,"","","",""]
0206     pass
0207
0208     lcc[:-1] = ecc[1:]
0209     lcc[-1] = [0,0,1,1]
0210
0211     fcc = Flight(ecc,lcc,ucc,ccc)
0212
0213     ##################################################
0214
0215     f = Flight.combine( [fzz, faa, fbb, fcc] )
0216     f.save()
0217     #elu = f.elu
0218     #np.save("/tmp/flightpath.npy", elu )
0219
0220     #print(elu[:,3,:4].copy().view("|S2"))
0221
0222     if plt:
0223         f.quiver_plot(ax, sc=sc)
0224         fig.show()
0225     pass
0226
0227