eic-opticks/ana/torus_hyperboloid.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 import numpy as np
0023
0024
0025 class Tor(object):
0026     def __init__(self, R, r):
0027         self.R = R
0028         self.r = r
0029
0030     def __repr__(self):
0031         return "Tor r:%s R:%s " % (self.r, self.R )
0032
0033     def rz(self, _z):
0034         z = np.asarray(_z)
0035         R = self.R
0036         r = self.r
0037         return R - np.sqrt(r*r-z*z)
0038
0039
0040 class Hyp(object):
0041     def __init__(self, r0, zf ):
0042         """
0043         :param r0: waist radius, ie radius at z=0
0044         :param zf: param as returned by ZF
0045         """
0046         self.r0 = r0
0047         self.zf = zf
0048
0049     def _stereo(self):
0050         """
0051         X4Solid::convertHype_::
0052
0053             1082     /*
0054             1083      Opticks CSG_HYPERBOLOID uses
0055             1084                 x^2 +  y^2  =  r0^2 * (  (z/zf)^2  +  1 )
0056             1085
0057             1086      G4Hype uses
0058             1087                 x^2 + y^2 = (z*tanphi)^2 + r^2
0059             1088                 x^2 + y^2 =  r0^2 * ( (z*tanphi/r0)^2 + 1 )
0060             1089
0061             1090      So
0062             1091                tanphi/r0 = 1/zf
0063             1092
0064             1093                zf = r0/tanphi
0065             1094
0066             1095                tanphi = r0/zf
0067             1096
0068             1097             stereo = phi = arctan(r0/zf)
0069             1098
0070             1099     */
0071         """
0072         return np.arctan2( self.r0, self.zf )
0073
0074     stereo = property(lambda self:self._stereo())
0075
0076
0077     @classmethod
0078     def ZF(cls, r0, zw, rw ):
0079         """
0080         :param r0: waist radius, ie radius at z=0
0081         :param zw: z at which to pin the radius
0082         :param rw: target radius at z=zw
0083         :return zf: hyperboloid zf param to hit target radius w, at z=zw
0084         """
0085         rr0 = r0*r0
0086         ww = rw*rw
0087         return zw*np.sqrt(rr0/(ww-rr0))
0088
0089     def __repr__(self):
0090         return "Hyp r0:%s zf:%s stereo(radians):%s  " % (self.r0, self.zf, self.stereo )
0091
0092     def rz(self, _z):
0093         z = np.asarray(_z)
0094         r0 = self.r0
0095         zf = self.zf
0096         zs = z/zf
0097         return r0*np.sqrt( zs*zs + 1 )
0098
0099
0100 if __name__ == '__main__':
0101
0102
0103     R,r = 97.000,52.010
0104
0105     r0 = R - r
0106     rr0 = r0*r0
0107
0108     tor = Tor(R,r)
0109     assert tor.rz(0) == R - r
0110     assert tor.rz(r) == R
0111     assert np.all(tor.rz([0,r]) == np.asarray( [R-r, R] ) )
0112
0113     print tor
0114
0115
0116     zw = r
0117     rw = R
0118     zf = Hyp.ZF( r0, zw, rw )
0119     hyp = Hyp( r0, zf )
0120
0121     assert hyp.rz(0) == r0
0122     assert hyp.rz(zw) == rw
0123
0124     z = np.linspace(0,zw, 100)
0125     print hyp.rz(z)
0126
0127
0128     print "zf", zf
0129
0130
0131
0132
0133
0134
0135