eic-opticks/CSG/sympy_cylinder.py

0001 #!/usr/bin/env python
0002 """
0003 sympy_cylinder.py
0004 ===================
0005
0006 Recast ray intersection with cylinder to requiring that the distance
0007 from cylinder axis line AB to a point C must be the cylinder radius.::
0008
0009
0010          [0,0,z2]  B +--r-+ C  (o + t v)
0011                      |   /
0012                      |  /
0013                      | /
0014                      |/
0015                    A +   [ 0,0,z1 ]
0016
0017
0018
0019                    B + [0,0,z2]
0020                      |\
0021                      | \
0022                      |r P  (o + t v)
0023                      | /
0024                      |/
0025                    A + [ 0,0,z1 ]
0026
0027 * Area of triangle ABP with height r, |BA|d/2
0028 * Area spanned by vectors u (BA) and v (BP), |u x v|/2
0029 * equating areas gives: r = |BAxBP| / |BA|
0030
0031 But the area can be given by another choice of side vectors too, so:
0032
0033 * r = |APxBP|/|BA|
0034
0035
0036 * AP = o+tv - [0,0,z1]
0037 * BP = o+tv - [0,0,z2]
0038 * |BA| = (z2-z1)    (z2 > z1 by definition)
0039
0040 * r*r = APxBP.dot(APxBP) / (z2-z1)**2
0041
0042 * fn(t) = APxBP.dot(APxBP) / (z2-z1)**2 - r*r
0043
0044
0045 Comparing the a,b,c coeff from nmskTailOuterIITube with the old ones get factor of 20800.7383 in all three.
0046 That is dz*dz::
0047
0048     In [12]: (72.112*2)*(72.112*2)
0049     Out[12]: 20800.562175999996
0050
0051
0052
0053 Endcap intersects of the ray, p  = o + t v
0054 Plane eqn, n = [0,0,1]     p0 = [0,0,z1]  OR [0,0,z2]
0055
0056     (p - p0).n = 0
0057
0058 This says that a vector in the plane which has no z component::
0059
0060     oz + t vz  - z1 = 0
0061
0062     t = (z1 - oz)/vz
0063     t = (z2 - oz)/vz
0064
0065
0066
0067
0068
0069 """
0070
0071 import numpy as np, sympy as sp
0072 from collections import OrderedDict as odict
0073
0074 def pp(q, q_label="?", note=""):
0075
0076     tq = type(q)
0077     if tq is np.ndarray:
0078         q_type = "np"
0079     elif q.__class__.__name__  == 'MutableDenseMatrix':
0080         q_type = "sp.MutableDenseMatrix"
0081     else:
0082         q_type = "%s.%s" % ( tq.__module__ , tq.__name__ )
0083     pass
0084
0085     sh = str(q.shape) if hasattr(q,"shape") else ""
0086
0087     print("\n%s : %s : %s : %s \n" % (q_label, q_type, sh, note) )
0088
0089     if q_type.startswith("sp"):
0090         sp.pprint(q)
0091     elif q_type.startswith("np"):
0092         print(q)
0093     else:
0094         print(q)
0095     pass
0096
0097
0098 #import sympy.physics.vector as spv
0099
0100 t = sp.symbols("t")
0101 ox, oy, oz = sp.symbols("ox oy oz")
0102 vx, vy, vz = sp.symbols("vx vy vz")
0103 r,z1,z2 = sp.symbols("r z1 z2")
0104
0105 A = sp.Matrix([ [0], [0], [z1] ])
0106 B = sp.Matrix([ [0], [0], [z2] ])
0107 O = sp.Matrix([ [ox], [oy], [oz] ])  # column vector
0108 V = sp.Matrix([ [vx], [vy], [vz] ])
0109
0110 P = O + t*V   # parametric ray eqn
0111
0112 AP = P - A
0113 BP = P - B
0114 AB = B - A
0115
0116 pp(P, "P")
0117 pp(AP, "AP")
0118 pp(BP, "BP")
0119 pp(AB, "AB")
0120
0121 APBP = AP.cross(BP)
0122 pp(APBP, "APBP", "AP.cross(BP) is perpendicular to axis, ie with zero z-component" )
0123
0124 apbp = sp.sqrt(APBP.dot(APBP))
0125 pp(apbp, "apbp", "magnitude of the cross product, cross product proportional to ABC triangle area" )
0126
0127 ab = z2 - z1
0128
0129 fn = (apbp*apbp)/(ab*ab) - r*r
0130 #fn = apbp  - r*r*ab*ab
0131
0132 pp(fn, "fn")
0133
0134
0135 collect_expand_fn = sp.collect( sp.expand(fn), t )
0136 pp(collect_expand_fn, "collect_expand_fn")
0137
0138 tfn = odict()
0139 sfn = odict()
0140
0141 for i in range(3):
0142     tfn[i] = collect_expand_fn.coeff(t, i)
0143     sfn[i] = sp.simplify(tfn[i])
0144 pass
0145
0146 for i in range(3): pp(tfn[i], "tfn[%d]"% i )
0147 for i in range(3): pp(sfn[i], "sfn[%d]"% i )
0148