eic-opticks/ana/x018_torus_hyperboloid_plt.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 Continuing from tboolean-12
0023
0024 See also
0025
0026 x018.py x019.py x020.py x021.py
0027     these were manually created from the generated g4code
0028     but could be generated too
0029
0030 shape.py
0031     rationalizations replacing tubs-torus with cons
0032
0033 notes/issues/torus_replacement_on_the_fly.rst
0034
0035
0036 // arghh : G4Hype cannot be z-cut : so would have to intersect with cylinder
0037 // this makes me want to use a cone for the neck instead
0038
0039 """
0040
0041 import logging
0042 log = logging.getLogger(__name__)
0043 import numpy as np, math
0044
0045
0046 import matplotlib.pyplot as plt
0047 from opticks.ana.torus_hyperboloid import Tor, Hyp
0048
0049 from opticks.ana.shape import Shape, ellipse_closest_approach_to_point
0050
0051 class Ellipsoid(Shape):pass
0052 class Tubs(Shape):pass
0053 class Torus(Shape):pass
0054 class Cons(Shape):pass
0055 class Hype(Shape):pass
0056 class UnionSolid(Shape):pass
0057 class SubtractionSolid(Shape):pass
0058 class IntersectionSolid(Shape):pass
0059
0060
0061
0062 class X018(object):
0063     """
0064     AIMING TO PHASE THIS OUT, use instead x018.py
0065
0066     G4VSolid* make_solid()
0067     {
0068         G4ThreeVector A(0.000000,0.000000,-23.772510);
0069         G4ThreeVector B(0.000000,0.000000,-195.227490);
0070         G4ThreeVector C(0.000000,0.000000,-276.500000);
0071         G4ThreeVector D(0.000000,0.000000,92.000000);
0072
0073
0074         G4VSolid* d = new G4Ellipsoid("PMT_20inch_inner_solid_1_Ellipsoid0x4c91130", 249.000000, 249.000000, 179.000000, -179.000000, 179.000000) ; // 3
0075         G4VSolid* g = new G4Tubs("PMT_20inch_inner_solid_2_Tube0x4c91210", 0.000000, 75.951247, 23.782510, 0.000000, CLHEP::twopi) ; // 4
0076         G4VSolid* i = new G4Torus("PMT_20inch_inner_solid_2_Torus0x4c91340", 0.000000, 52.010000, 97.000000, -0.000175, CLHEP::twopi) ; // 4
0077
0078         G4VSolid* f = new G4SubtractionSolid("PMT_20inch_inner_solid_part20x4cb2d80", g, i, NULL, A) ; // 3
0079         G4VSolid* c = new G4UnionSolid("PMT_20inch_inner_solid_1_20x4cb30f0", d, f, NULL, B) ; // 2
0080
0081         G4VSolid* k = new G4Tubs("PMT_20inch_inner_solid_3_EndTube0x4cb2fc0", 0.000000, 45.010000, 57.510000, 0.000000, CLHEP::twopi) ; // 2
0082
0083         G4VSolid* b = new G4UnionSolid("PMT_20inch_inner_solid0x4cb32e0", c, k, NULL, C) ; // 1
0084         G4VSolid* m = new G4Tubs("Inner_Separator0x4cb3530", 0.000000, 254.000000, 92.000000, 0.000000, CLHEP::twopi) ; // 1
0085
0086         G4VSolid* a = new G4IntersectionSolid("PMT_20inch_inner1_solid0x4cb3610", b, m, NULL, D) ; // 0
0087         return a ;
0088     }
0089
0090
0091                            a
0092
0093                       b             m(D)
0094
0095                 c          k(C)
0096
0097             d     f(B)
0098
0099                 g(B)  i(B+A)
0100
0101     """
0102     def __init__(self, mode=0):
0103         A = np.array( [0, -23.772510] )
0104         C = np.array( [0, -276.500000] )
0105         D = np.array( [0, 92.000000] )
0106
0107         B0 = np.array( [0, -195.227490] )
0108         z0 = B0[1] + A[1]    ## offset of torus (ellipsoid frame)
0109
0110         print("z0 %s offset of torus (ellipsoid frame)" % z0 )
0111
0112         self.z0 = z0
0113
0114         i = Torus("i", [ 52.010000, 97.000000] )
0115         # 97-51.01 = 44.99 ( slightly inside the tube radius 45.01 ? accident or CSG planning ?)
0116         r = i.param[0]
0117         R = i.param[1]
0118
0119         self.r = r
0120         self.R = R
0121
0122         d = Ellipsoid("d", [249.000, 179.000 ] )
0123         ex = d.param[0]
0124         ez = d.param[1]
0125
0126         self.ex = ex
0127         self.ez = ez
0128
0129         k = Tubs("k", [45.010000, 57.510000] )
0130         bx = k.param[0]
0131         self.bx = bx
0132
0133
0134         p = ellipse_closest_approach_to_point( ex, ez, [R,z0] ) # center of torus RHS circle
0135         print(" p %s " % repr(p) )
0136
0137         pr = p[0]
0138         pz = p[1]    # z of ellipse to torus-circle closest approach point (ellipsoid frame)
0139
0140         self.p = p
0141
0142
0143         if mode == 0:
0144             g = Tubs("g", [75.951247,23.782510] )
0145             f = SubtractionSolid("f", [g,i,A] )
0146             B = B0
0147         elif mode == 1:
0148             r2 = 83.9935              # at torus/ellipse closest approach point : huh same as pr no (maybe not quite)
0149             r1 = R - r
0150
0151             mz = (z0 + pz)/2.         # mid-z cone coordinate (ellipsoid frame)
0152             hz = (pz - z0)/2.         # cons half height
0153
0154             f = Cons("f", [r1,r2,hz ] )
0155             B = np.array( [0, mz] )
0156
0157             print(" replacment Cons %s offset %s " % (repr(f),repr(B)))
0158
0159         elif mode == 2:
0160
0161             r0 = R - r
0162             rw,zw = pr,pz-z0   # target the closest approach point
0163
0164             zf = Hyp.ZF( r0, zw, rw )
0165             hyp = Hyp( r0, zf )
0166
0167             stereo = hyp.stereo
0168             hhh = pz - z0    # hype-half-height
0169
0170             f = Hype("f", [r0,stereo,hhh] )
0171             B = np.array( [0, z0] )
0172
0173         else:
0174             assert 0
0175         pass
0176
0177         c = UnionSolid("c", [d,f,B] )
0178
0179         m = Tubs("m", [254.000000, 92.000000] )
0180         b = UnionSolid("b", [c,k,C] )
0181         a = IntersectionSolid("a", [b,m,D] )
0182
0183         self.root = a
0184         self.sz = 400
0185
0186         self.ellipse = d
0187         self.torus = i
0188         #self.cyneck = g
0189         self.f = f
0190     pass
0191
0192
0193
0194 if __name__ == '__main__':
0195
0196     #x = X018(mode=0)  # crazy cylinder-torus neck
0197     x = X018(mode=1)  # cone neck
0198     #x = X018(mode=2)  # hype neck
0199
0200     print "x.f ", x.f
0201
0202
0203     sz = x.sz
0204
0205     r = x.r
0206     R = x.R
0207     z0 = x.z0   # natural torus frame offset in ellipse frame
0208
0209
0210
0211
0212     #ch = x.cyneck.param[1]*2   # full-height of cylinder neck
0213     #cr = x.cyneck.param[0]     # half-width of cylinder neck, ie the radius
0214
0215     ex = x.ex
0216     ez = x.ez
0217
0218     p = x.p  # closest approach of RHS torus circle center to a point on the ellipse
0219
0220
0221     tor = Tor(R,r)
0222     assert tor.rz(0) == R - r
0223     assert tor.rz(r) == R
0224     assert np.all(tor.rz([0,r]) == np.asarray( [R-r, R] ) )
0225
0226
0227     tz0 = 0
0228     tz1 = r      # how far in z to draw the torus and hype parametric lines
0229     #tz1 = ch
0230
0231     tz = np.linspace( tz0, tz1, 100)
0232     tr = tor.rz(tz)
0233
0234
0235     r0 = R - r     # Hyperboloid waist radius at z=0 (natural frame)
0236     # defining parameter of Hyperboloid by targetting an rz-point (natural frame)
0237
0238     #rw,zw = R, r   # target the top of the torus : gives wider neck than  cy-tor
0239     rw,zw = p[0],p[1]-z0   # target the closest approach point
0240     halfZlen = p[1]-z0
0241
0242     zf = Hyp.ZF( r0, zw, rw )
0243     hyp = Hyp( r0, zf )
0244     print hyp
0245     print "hyp halfZLen ", halfZlen
0246
0247     hr = hyp.rz(tz)
0248
0249     tz += z0
0250
0251
0252     plt.ion()
0253     fig = plt.figure(figsize=(6,5.5))
0254     plt.title("x018_torus_hyperboloid_plt")
0255
0256     ax = fig.add_subplot(111)
0257
0258     zoom = False
0259     #zoom = True
0260
0261     if zoom:
0262         zmp = np.array( [R, z0, 1.5*r, 1.5*r] )
0263         ax.set_xlim([zmp[0]-zmp[2],zmp[0]+zmp[2]])
0264         ax.set_ylim([zmp[1]-zmp[3],zmp[1]+zmp[3]])
0265     else:
0266         ax.set_ylim([-350,200])
0267         ax.set_xlim([-300,300])
0268     pass
0269
0270     for pt in x.root.patches():
0271         print "pt ", pt
0272         ax.add_patch(pt)
0273     pass
0274
0275     #rhs = Circle( [R,z0], radius=r )
0276     #ax.add_patch(rhs)
0277     #ell = Ellipse( [0,0], 2*ex, 2*ez )
0278     #ax.add_patch(ell)
0279
0280     #ax.scatter( e[:,0], e[:,1], s=1 )
0281
0282     ax.scatter( p[0], p[1] , marker="*")
0283     ax.scatter( -p[0], p[1] , marker="*" )
0284
0285     ax.plot( [R, p[0]], [z0, p[1]] )
0286     ax.plot( [-R, -p[0]], [z0, p[1]] )
0287
0288     ax.plot( [R-r, p[0]], [z0, p[1]] )
0289     ax.plot( [-R+r, -p[0]], [z0, p[1]] )
0290
0291     ax.plot(  tr, tz )
0292     ax.plot( -tr, tz )
0293
0294     ax.plot(  hr, tz , linestyle="dashed")
0295     ax.plot( -hr, tz , linestyle="dashed")
0296
0297
0298
0299
0300     fig.show()
0301
0302