eic-opticks/analytic/cubic.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 import logging
0024 log = logging.getLogger(__name__)
0025 from opticks.ana.base import opticks_main
0026 from opticks.analytic.csg import CSG
0027
0028 import matplotlib.pyplot as plt
0029 import mpl_toolkits.mplot3d.axes3d as axes3d
0030 #import pylab as plt
0031
0032
0033 class Cubic(object):
0034     def __init__(self, A, B, C, D, z1, z2):
0035         self.A = A
0036         self.B = B
0037         self.C = C
0038         self.D = D
0039         self.z1 = z1
0040         self.z2 = z2
0041
0042     def __repr__(self):
0043         return "Cubic(A=%s,B=%s,C=%s,D=%s,z1=%s,z2=%s)" % (self.A, self.B, self.C, self.D, self.z1, self.z2 )
0044
0045     def as_csg(self):
0046         return CSG.MakeCubic(A=self.A, B=self.B, C=self.C, D=self.D, z1=self.z1, z2=self.z2)
0047
0048     def rrmax(self):
0049         """
0050         ::
0051
0052              v = A z**3 + B z**2 + C z + D
0053
0054              dv/dz =   3 A z**2 + 2 B z + C
0055
0056         """
0057         A,B,C,D,z1,z2 = self.A, self.B, self.C, self.D, self.z1, self.z2
0058
0059         d = 3*A
0060         disc = B*B - d*C
0061
0062         vals = np.zeros( 4 )
0063         vals[0] = self.rrval(z1)
0064         vals[1] = self.rrval(z2)
0065
0066         if disc > 0:
0067             sdisc = np.sqrt(max(disc,0.))
0068             q = -(B + sdisc) if B > 0 else -(B - sdisc)
0069
0070             e1 = q/d
0071             e2 = C/q
0072
0073             vals[2] = self.rrval(e1) if e1 > z1 and e1 < z2 else 0.
0074             vals[3] = self.rrval(e2) if e2 > z1 and e2 < z2 else 0.
0075         pass
0076         return vals.max()
0077
0078
0079     def rrval(self, z ):
0080         A,B,C,D = self.A, self.B, self.C, self.D
0081         return (((A*z+B)*z + C)*z) + D
0082
0083     def rval(self, z ):
0084         return np.sqrt(self.rrval(z))
0085
0086
0087     def plot_profile(self, ax):
0088         z = np.linspace(self.z1,self.z2,100)
0089
0090         zlen = self.z2 - self.z1
0091         rrmx = self.rrmax()
0092
0093         rmi = 0
0094         rmx = np.sqrt(rrmx)
0095
0096         rlen = rmx - rmi
0097         rlenp = rlen+rlen/10.
0098
0099         ax.set_ylim(self.z1-zlen/10.,self.z2+zlen/10.)
0100         ax.set_xlim( -rlenp , rlenp )
0101
0102         rv = self.rval(z)
0103
0104         ax.plot(rv, z, 'k')
0105         ax.plot(-rv, z, 'k')
0106
0107         ax.plot((-rv[0],rv[0]),(z[0],z[0]), 'r-')
0108         ax.plot((-rv[-1],rv[-1]),(z[-1],z[-1]), 'r-')
0109
0110
0111
0112     def plot_sor(self, ax):
0113         # https://stackoverflow.com/questions/36464982/ploting-solid-of-revolution-in-python-3-matplotlib-maybe
0114
0115         A,B,C,D,z1,z2 = self.A, self.B, self.C, self.D, self.z1, self.z2
0116
0117         u = np.linspace(z1,z2,100)
0118         v = np.linspace(0, 2*np.pi, 60)
0119
0120         U, V = np.meshgrid(u, v)
0121
0122         R = self.rval(U)
0123
0124         X = R*np.cos(V)
0125         Y = R*np.sin(V)
0126         Z = U
0127
0128         ax.plot_surface(X, Y, Z, alpha=0.3, color='red', rstride=6, cstride=12)
0129
0130
0131
0132
0133
0134 if __name__ == '__main__':
0135
0136     args = opticks_main(csgpath="$TMP/cubic_py")
0137
0138     CSG.boundary = args.testobject
0139     CSG.kwa = dict(poly="IM", resolution="50")
0140
0141     container = CSG("box", param=[0,0,0,200], boundary=args.container, poly="MC", nx="20" )
0142
0143     cubic = Cubic(A=1, B=10, C=100, D=400, z1=-10, z2=10)
0144     log.info("cubic:%r " % cubic )
0145     a = cubic.as_csg()
0146
0147     #zrrs = [[-100,30],[-50,80],[50,30],[100,100]]
0148     #a = CSG.MakeCubicBezier(zrrs)
0149
0150     CSG.Serialize([container, a], args.csgpath )
0151
0152
0153     plt.ion()
0154     plt.close()
0155
0156     fig = plt.figure()
0157
0158     ax = fig.add_subplot(121)
0159     cubic.plot_profile(ax)
0160
0161     ax = fig.add_subplot(122, projection='3d')
0162     cubic.plot_sor(ax)
0163
0164
0165     plt.show()
0166
0167
0168
0169