## PointPlane3D.py
## Version 1.00 (
## Copyright (c)
2007 Bruce Vaughan, BV Detailing & Design, Inc.
## All rights
reserved.
## NOT FOR
################################################################################
import math
'''
Classes:
Point
dot()
cross()
mag(), uv()
dist(), format()
fformat()
equiv()
eq()
Plane3D
DistToPlane()
PointRotate3D()
ThreePtCircle()
data1()
data2()
data3()
LineLineIntersect3D
DistancePointLine3D
BasisTransToGlobal
translate()
BasisTransToLocal
trans_to_local()
Functions:
simplify
gcd
fifDim
cross_product
dot_product
mag
unitV
polarPt
polarPtXY
polarPtXZ
polarPtRotate
midPt
spaPts
fixSpaPts
autoSpaPts
polarSpaPts
polarFixSpaPts
polarAutoSpaPts
polarEllipse
sortPts
angBetweenLines
angBetweenLinesProjZ
endWP
memWPs
memPlane
Changes:
11/14/07 - Modify
Point() method __cmp__()
'''
epsilon = 0.000001
epsilon1 = 0.001
## Point
########################################
class Point(object):
'''
p1 = Point(1,2,3)
Point(1.000000, 2.000000, 3.000000)
attributes: x, y, z, data
methods: dot(), cross(), mag(), uv(),
dist(), format(), fformat(), equiv(), eq()
'''
def __init__(self, x=0.0, y=0.0, z=0.0):
self.x = float(x)
self.y = float(y)
self.z = float(z)
self.data = (self.x, self.y, self.z)
# cannot pickle bool objects in Python
2.3 (True, False)
self.__init = 1
def __str__(self):
return '(%0.4f, %0.4f, %0.4f)' %
self.data
def __repr__(self):
return 'Point(%f, %f, %f)' % self.data
def format(self):
return 'Point(%s, %s, %s)' %
tuple(map(fifDim, self))
def fformat(self):
return '%0.4f, %0.4f, %0.4f' %
tuple(self)
'''
Operator overloading
'other' can be a Point object, tuple or
list.
'f' can be an integer or float except as
noted.
'''
def __add__(self, other):
# return Point(*[a+b for a,b in
zip(self,other)])
# following code is more efficient
return Point(self.x+other[0],
self.y+other[1], self.z+other[2])
def __iadd__(self, other):
return Point(self.x+other[0], self.y+other[1],
self.z+other[2])
def __radd__(self, other):
return Point(other[0]+self.x,
other[1]+self.y, other[2]+self.z)
def __sub__(self, other):
return Point(self.x-other[0],
self.y-other[1], self.z-other[2])
def __isub__(self, other):
return Point(self.x-other[0],
self.y-other[1], self.z-other[2])
def __rsub__(self, other):
return Point(other[0]-self.x,
other[1]-self.y, other[2]-self.z)
def __mul__(self, f):
return Point(self.x*f, self.y*f,
self.z*f)
def __imul__(self, f):
return Point(self.x*f, self.y*f,
self.z*f)
def __div__(self, f):
return Point(self.x/f, self.y/f,
self.z/f)
def __imul__(self, f):
return Point(self.x/f, self.y/f,
self.z/f)
def __neg__(self):
return Point(-self.x, -self.y, -self.z)
def __pow__(pt, pwr):
return Point(pt.x**pwr, pt.y**pwr,
pt.z**pwr)
def __iter__(self):
'''
Point iteration:
for a in pt:
print a,
Prints:
pt.x pt.y pt.z
'''
for i in self.data:
yield i
def __getitem__(self, i):
return self.data[i]
def __setitem__(self, i, value):
try:
object.__setattr__(self,
{0:'x',1:'y',2:'z'}[i], float(value))
object.__setattr__(self, 'data',
(float(self.x), float(self.y), float(self.z)))
except KeyError:
raise IndexError, "Index
argument out of range in '%s' __setitem__" % (type(self).__name__)
except ValueError:
raise ValueError, "Invalid
literal in '%s' __setitem__" % (type(self).__name__)
def __setattr__(self, name, value):
if not self.__dict__.has_key('_Point__init'):
return object.__setattr__(self,
name, value)
elif self.__dict__.has_key(name):
object.__setattr__(self, name,
value)
object.__setattr__(self, 'data',
(float(self.x), float(self.y), float(self.z)))
else:
raise AttributeError, "'%s'
object has no attribute '%s'" % (type(self).__name__, name)
def __cmp__(self, other, epsilon=0.000001):
x = abs(self.x-other.x)
if x < epsilon:
y = abs(self.y-other.y)
if y < epsilon:
return cmp(self.z, other.z)
return cmp(self.y, other.y)
return cmp(self.x, other.x)
def dot(self, other):
'''
Return the dot product of two Points
(self, other)
'other' must be a Point object
'''
# return sum([a*b for a,b in
zip(self,other)])
return (self.x*other.x + self.y*other.y
+ self.z*other.z)
def uv(self):
'''
Return the unit vector of a Point
object.
If vector has no magnitude, returns
Point(0,0,0).
'''
m = self.mag()
try: return Point(self.x/m, self.y/m,
self.z/m)
except: return Point()
def cross(self, other):
'''
Return the cross product of two Points
(self, other)
'other' must be a Point object
'''
return Point(self.y*other.z -
self.z*other.y, self.z*other.x - self.x*other.z, self.x*other.y -
self.y*other.x)
def mag(self):
'''
Return the scalar length of a Point
instance vector
'''
return (self.x**2 + self.y**2 +
self.z**2)**0.5
def dist(self, other):
'''
Return the scalar distance between two
Points (self, other)
'other' can be a Point object, list or
tuple
'''
# return sum((self-other)**2)**0.5
p = self-other
return (p.x**2 + p.y**2 + p.z**2)**0.5
def equiv(self, other, epsilon=0.000001):
'''
Compare two vectors, return a tuple of
1, 0, or -1 values
Interactive example:
>>>
Point(1,1,5).equiv(Point(2,0,5))
(-1, 1, 0)
>>>
'''
def comp(x,y):
if abs(x-y) < epsilon: return 0
elif x > y: return 1
else: return -1
return tuple(map(comp, self, other))
def eq(self, other):
'''
Wrapper for equiv() method
If the points/vectors are equal or
equivalent, return True
Otherwise, return False
>>> p11 = Point(10,10,10)
>>> p12 = Point(1,1,1)
>>> p11.uv().eq(p12.uv())
True
>>> p11.eq(p12)
False
>>>
'''
a = self.equiv(other)
if 1 in a or -1 in a:
return False
return True
## Plane3D
######################################
class
Plane3D(object):
'''
Define a plane in 3D space from 3
non-collinear point objects (Plane 1)
Creates a Plane3D instance:
a = Plane3D(p1, p2, p3)
Methods:
ThreePtCircle(self) -
Calculate center point and radius
of an arc passing through
instance points p1, p2, p3
PointRotate3D(self, p, theta,
ctr_pt=False) -
Given a center point object
(ctr_pt) and angle (theta) in radians,
rotate point object 'p' about
ctr_pt and ctr_pt+self.N0 by theta
and return a point object
representing the new location. 'ctr_pt'
defaults to instance point p1.
DistToPlane(pt) -
Returns a signed distance from
point object 'pt' to instance Plane 1
data1(), data2() and data3() -
Each returns a string for display
of instance data
'''
def __init__(self, p1, p2, p3, theta=0):
self.p1 = p1
self.p2 = p2
self.p3 = p3
self.pts = (p1,p2,p3)
if False not in [isinstance(p, Point)
for p in [p1,p2,p3]]:
'''
p3
.
/ \\
/ \\
e,e0 G. \\
/ \\
/ \\
p1.-----.-----. p2
F
d,d0
Plane 1
Define a plane from 3 non-collinear
points
Ax + By + Cz + D = 0
The normal 'N' to the plane is the
vector (A, B, C)
'''
self.N, A, B, C, self.D =
self.__plane_def(p1, p2, p3)
if mag(self.N) > epsilon:
self.N0 = unitV(self.N)
else:
self.N0 = Point(0, 0, 0)
raise ValueError, 'INVALID -
The points are collinear.'
'''
If vector N is the normal to the
plane then all points 'p' on the
plane satisfy the following:
N dot p = k where 'dot' is the dot
product
N dot p = N.x*p.x + N.y*p.y + N.z*p.z
'''
# Plane constant 'k'
self.k = self.N.dot(p1)
# Displacement of the plane from
the origin (Point(0,0,0))
self.k0 = self.N0.dot(p1)
'''
Instance point object p1, unit
vector normal to the plane, and
the plane displacement from the
origin of the coordinate system
'''
self.data = (self.p1, self.N0,
self.k0)
'''
Vector e (p1 to p3) and unit vector
e0
Vector d (p1 to p2) and unit vector
d0
Point F, midpoint on vector d
Point G, midpoint on vector e
'''
self.e = p3 - p1
self.e0 = self.e.uv()
self.d = p2 - p1
self.d0 = self.d.uv()
self.F = p1 + self.d / 2
self.G = p1 + self.e / 2
'''
self.Ra = Distance between points
p1 and p2 (arc radius)
self.Q = Net angle between vectors
d0 and e0 (radius sector) in radians
self.pp = Point to point distance
along arc with radius self.Ra in sector self.Q
'''
self.Ra = p2.dist(p1)
if abs(theta - math.pi) < epsilon:
self.Q = math.pi
else:
self.Q =
math.acos(self.d0.dot(self.e0))
self.pp = abs(self.Q *
self.Ra)
else:
raise TypeError, "All
arguments passed to Plane3D must be <class '__main__.Point'>"
def __plane_def(self, p1, p2, p3):
'''
Plane definition method
'''
N = cross_product(p2-p1, p3-p1)
A = N.x
B = N.y
C = N.z
D = dot_product(-N, p1)
return N, A, B, C, D
def __str__(self):
return '<%s.Plane3D object>' %
(__name__)
def __repr__(self):
return 'Plane3D(%s, %s, %s,
theta=%0.6f)' % (repr(self.p1), repr(self.p2), repr(self.p3), self.Q)
def ThreePtCircle(self):
'''
The intersection of three planes is
either a point, a line, or there is no intersection.
Three planes can be written as:
N1 dot p = d1
N2 dot p = d2
N3 dot p = d3
'Nx' is the normal vector
'p' is a point on the plane
'dot' signifies the dot product of 'Nx'
and 'p'
'dx' = plane constant (displacement of
the plane from the origin if Nx is a unit vector)
The intersection point of the three
planes "M" is given by:
M = (d1*(N2 cross N3) + d2(N3 cross N1)
+ d3*(N1 cross N2)) / (N1 dot (N2 cross N3))
'cross' indicates the cross product and
'dot' indicates the dot product
Calculate the center point of the
circle (intersection point of three planes) and the radius.
Plane 1 is defined in __init__ method.
Plane 2 passes through self.G (midpoint
on vector 'e'), self.G + instance unit vector
and the cross product of instance unit
vector and the unit vector of 'e'
Plane 3 passes through self.F (midpoint
on vector 'd'), self.F + instance unit vector
and the cross product of instance unit
vector and the unit vector of 'f'
N2 = self.d0
N3 = self.e0
'''
self.d2 = dot_product(self.e0, self.G)
self.d3 = dot_product(self.d0, self.F)
N23 = cross_product(self.e0, self.d0)
N31 = cross_product(self.d0, self.N0)
N12 = cross_product(self.N0, self.e0)
NdN23 = dot_product(self.N0, N23)
numer = N23*self.k0 + N31*self.d2 +
N12*self.d3
if abs(NdN23) > epsilon:
'''Center point of 3D arc'''
self.M = numer/NdN23
'''Radius of 3D arc'''
self.R = self.M.dist(self.p1)
else:
self.M = Point(0.0, 0.0, 0.0)
self.R = 0.0
def PointRotate3D(self, p, theta,
ctr_pt=False):
'''
Rotate point p about a line passing
through 'ctr_pt' and normal
to instance Plane 1 by the angle
'theta' in radians.
A right hand coordinate system is
assumed. Positive angles are
counter-clockwise when looking down the
axis toward the origin.
Return the new point. 'ctr_pt' defaults
to 'self.p1'
Interactive example:
>>> c = Plane3D(Point(144, 256,
300), Point(324, 587, 398), Point(645, 400, 130))
>>>
c.PointRotate3D(Point(4003.4, 7645.3254, 18745.48867), 1.0, Point(5,5,5))
Point(16840.039750, 5048.060568,
10808.548819)
>>>
c.PointRotate3D(Point(4003.4, 7645.3254, 18745.48867), 1.0)
Point(16590.661551, 5193.988794,
11017.123729)
>>> c.N0
Point(-0.400516, 0.453529, -0.796177)
>>>
c.N0 is the rotation axis unit
vector
'''
return polarPtRotate((ctr_pt or
self.p1), (ctr_pt or self.p1)+self.N0, p, theta)
def DistToPlane(self, p):
'''
Given any point 'a' on a plane: N dot
(a-p) = 0
Return a signed distance from 'p' to
Plane 1
A positive distance is measured in the
same direction as self.N0 (the
unit vector normal to the plane), and a
negative distance is measured
in the opposite direction.
Alternatively:
s = A*p.x + B*p.y + C*p.z + D
If s>0 the point p lies on the
same side as the plane normal 'N'
Interactive example:
>>> d = Plane3D(Point(0,0,20),
Point(20,0,20), Point(0,20,15))
>>>
d.DistToPlane(Point(0,0,0))
19.402850002906639
>>> d.N0
Point(0.000000, 0.242536, 0.970143)
>>>
'''
return self.k0 - dot_product(self.N0,
p)
def data1(self):
s1 = 'Instance of Class %s - Plane
Definition Data' % (type(self).__name__)
s2 = '\n Normal Unit Vector: %s' % self.N0.fformat()
s3 = '\n Displacement from Model Origin: %0.4f' %
(self.k0)
s4 = '\nPlane Origin and Unit Vector
Points'
s5 = '\n Origin point: %s' % (p1.format())
s6 = '\n Unit Vector d0: %s' % self.d0.fformat()
s7 = '\n Unit Vector e0: %s' % self.e0.fformat()
return '%s%s%s%s%s%s%s' % (s1, s2, s3,
s4, s5, s6, s7)
def data2(self):
return
'Plane3D(%s,\n%s%s,\n%s%s,\n%sSector (radians) = %0.6f)' % \
(p1.format(), ' '*8,
p2.format(), ' '*8, p3.format(), ' '*8, self.Q)
def data3(self):
s8 = 'Three Point Circle'
s9 = '\n Center point M: %s' % (self.M.format())
s10 = '\n Radius R: %s' % (fifDim(self.R))
s11 = '\n Arc length pp: %s' % (fifDim(self.pp))
s12 = '\n Sector Q: %0.6f' % (self.Q)
return '%s%s%s%s%s' % (s8, s9, s10,
s11, s12)
##
LineLineIntersect3D ##########################
class
LineLineIntersect3D(object):
"""
Calculate the points in 3D space Pa and Pb
that define the line segment which
<--> <-->
is the shortest route between two line
segments p1p2 and p3p4.
Each point occurs at the apparent
intersection of the 3D lines.
The apparent intersection is the location
where the two lines 'appear' to
intersect when viewed along the line
segment PaPb.
Pa lies on the line passing through p1p2.
Pb lies on the line passing through p3p4.
Equation for each line:
Pa = p1 + ma(p2-p1)
Pb = p3 + mb(p4-p3)
The shortest line segment is perpendicular
to both lines. Therefore:
(Pa-Pb).(p2-p1) = 0
(Pa-Pb).(p4-p3) = 0
Where:
'.' indicates the dot product
A = p1-p3
B = p2-p1
C = p4-p3
Substituting:
(A + ma(B) - mb(C)).B = 0 & (A + ma(B) - mb(C)).C = 0
-----------------------------------------------------------------
A.B + ma(B.B) - mb(C.B) = 0
A.B + ma(B.B) - (ma(C.B)-A.C)/C.C)(C.B) = 0
ma(B.B)(C.C) - ma(C.B)(C.B) =
(A.C)(C.B)-(A.B)(C.C)
ma = ((A.C)(C.B)-(A.B)(C.C))/((B.B)(C.C) -
(C.B)(C.B))
mb = (A.B + ma(B.B))/(C.B)
If the cross product of the two lines has
no magnitude, the lines are parallel.
A line extends forever in both
directions.
<-->
The name of a line passing through two
different points p1 and p2 would be "line p1p2" or p1p2.
The two-headed arrow over p1p2 signifies a
line passing through points p1 and p2.
Two lines which have no actual intersection
but are not parallel are called 'skew' or 'agonic' lines.
Skew lines can only exist in three or more
dimensions.
p4
.
/ <---member 2
/
Pb .
/ \\
/
\\<-----Shortest line segment
. \\
p3
.---.---------. <---member 1
p1 Pa
p2
"""
def __init__(self, p1, p2, p3, p4):
self.p1 = p1
self.p2 = p2
self.p3 = p3
self.p4 = p4
if False not in [isinstance(p, Point)
for p in [p1,p2,p3,p4]]:
A = p1-p3
B = p2-p1
C = p4-p3
'''
Line p1p2 and p3p4 unit vectors
'''
self.uv1 = (p2-p1).uv()
self.uv2 = (p4-p3).uv()
'''
Check for parallel lines
'''
if self.uv1.cross(self.uv2).mag()
>= epsilon:
ma = ((dot_product(A,
C)*dot_product(C, B)) - (dot_product(A, B)*dot_product(C, C)))/ \
((dot_product(B,
B)*dot_product(C, C)) - (dot_product(C, B)*dot_product(C, B)))
mb = (ma*dot_product(C, B) +
dot_product(A, C))/ dot_product(C, C)
self.ma = ma
self.mb = mb
'''
Calculate the point on line 1
that is the closest point to line 2
'''
Pa = p1 + B*ma
self.Pmem1 = Pa
self.Pa = Pa
'''
Calculate the point on line 2
that is the closest point to line 1
'''
Pb = p3 + C*mb
self.Pmem2 = Pb
self.Pb = Pb
'''
Calculate the length of line
segment Pa-Pb
'''
self.inters_dist = Pa.dist(Pb)
'''
Determine if point Pa is on the
line segment p1-p2
Determine point Pa position
with respect to points p1 and p2
p1 is LE and p2 is RE
'''
if ma >= 0-epsilon1 and ma
<= 1+epsilon1:
self.on_segment1 = 1
else:
self.on_segment1 = 0
if abs(ma) < epsilon1*10:
# apparent intersection is
at p1
self.position =
"LE"
elif abs(ma-1) <
epsilon1*10:
# apparent intersection is
at p2
self.position =
"RE"
elif ma < 0:
self.position =
"Beyond LE"
elif ma > 1:
self.position =
"Beyond RE"
elif ma <= 0.5:
self.position = "Not
Beyond LE"
else:
self.position = "Not
Beyond RE"
'''
Test the equivalence of unit
vectors p1->p2, p1->Pa and p2->Pa
'''
if
(p2-p1).uv().eq((Pa-p1).uv()): xl_dir = 1
else: xl_dir = -1
if (p2-p1).uv().eq((Pa-p2).uv()):
xr_dir = 1
else: xr_dir = -1
'''
Set the member 'X' direction
with respect to p1 and p2 - either '+' or '-'
'''
self.left_dist =
Pa.dist(p1)*xl_dir
self.right_dist =
Pa.dist(p2)*xr_dir
'''
Determine if point Pb is on the
line segment p3-p4
'''
if mb >= 0-epsilon1 and mb
<= 1+epsilon1:
self.on_segment2 = 1
else:
self.on_segment2 = 0
'''
Calculate the unit vector of
PaPb
'''
self.uv = unitV(Pb-Pa)
# Lines are parallel
else:
self.Pmem1, self.Pa = None,
None
self.Pmem2, self.Pb = None,
None
self.inters_dist = None
self.left_dist = None
self.right_dist = None
self.uv = None
else:
raise TypeError, "All
arguments passed to LineLineIntersect3D must be <class
'__main__.Point'>"
# Return False if lines are parallel, and
return True if lines are not parallel
def not_parallel(self):
if self.uv1.cross(self.uv2).mag() >=
epsilon:
return True
else:
return False
def __repr__(self):
return 'LineLineIntersect3D(%s, %s, %s,
%s)' % (repr(self.p1), repr(self.p2), repr(self.p3), repr(self.p4))
def __str__(self):
s1 = 'Instance of Class %s - End points
of Shortest Line Segment:' % (type(self).__name__)
s2 = '
Point on line 1 (self.Pa): %s' % (repr(self.Pa))
s3 = '
Point on line 2 (self.Pb): %s' % (repr(self.Pb))
s4 = '
Distance between points: %0.6f' % (self.inters_dist)
return '%s\n%s\n%s\n%s' % (s1, s2, s3,
s4)
def version(self):
"PointPlane3D.LineLineIntersect3D
Version 2.00"
##
DistancePointLine3D ##########################
class
DistancePointLine3D(object):
"""
Calculate the point in 3D space, Pa, that
defines the shortest line
<-->
segment between line p1p2 and point Pb.
Point Pb is closest to line p1p2 at an
intersecting perpendicular
line PaPb.
The dot product of two vectors, A and B,
will equal the cosine of
the angle between the vectors, times the
length of each vector.
A dot B = A.x * B.x + A.y * B.y + A.z * B.z
If the vectors are unit vectors, the dot
product is equal to the
cosine of the angle between the vectors.
Since the angle between lines p1p2 and PaPb
is 90 degrees, the dot
product of the normalized vector p1p2 and
normalized vector PaPb
is 0 (cos(90 deg)=0).
Determine location of point Pa and the
scalar distance from point
Pa to point Pb.
If the calculation result for 'u' is
between 0 and 1, Pa lies on
line segment p1p2.
Pb .
\\
\\<-----Shortest line
segment
\\
.---.---------.
p1 Pa
p2
"""
def __init__(self, p1, p2, Pb):
self.p1 = p1
self.p2 = p2
self.Pb = Pb
if False not in [isinstance(p, Point)
for p in [p1,p2,Pb]]:
u = (Pb-p1).dot(p2-p1) / sum((p2-p1)**2)
self.u = u
self.Pa = p1 + (p2-p1)*u
self.Pmem1 = self.Pa
'''
Calculate the length of line
segment PaPb
'''
self.dist = self.Pa.dist(Pb)
'''
Determine whether point Pa lies
between the line segment end points
or beyond one of the line segment
end points.
'''
if u >= 0-epsilon1 and u <=
1+epsilon1:
self.on_segment = 1
else:
self.on_segment = 0
'''
Set attribute 'position' indicating
point Pa position with respect
to the line segment end points as
follows:
'LE' indicates Pa occurs at p1
'RE' indicates Pa occurs at p2
'Beyond LE' indicates Pa occurs
beyond p1
'Beyond RE' indicates Pa occurs
beyond p2
'Not Beyond LE' indicates Pa
occurs in between p1 and p2 and
is closer to p1
'Not Beyond RE' indicates Pa
occurs in between p1 and p2 and
is closer to p2
'''
if abs(u) < epsilon1*10:
self.position = "LE"
elif abs(u-1) < epsilon1*10:
self.position = "RE"
elif u < 0:
self.position = "Beyond
LE"
elif u > 1:
self.position = "Beyond
RE"
elif u <= 0.5:
self.position = "Not
Beyond LE"
else:
self.position = "Not
Beyond RE"
'''
Test the equivalence of unit
vectors p1->p2, p1->Pa and p2->Pa
Calculate the scalar distance and
direction (beam member 'X'
distance) Pa occurs from each line
segment end point.
'''
if
sum((p2-p1).uv().equiv((self.Pa-p1).uv())) != 0: xl_dir = -1
else: xl_dir = 1
if
sum((p2-p1).uv().equiv((self.Pa-p2).uv())) != 0: xr_dir = -1
else: xr_dir = 1
self.left_dist =
self.Pa.dist(p1)*xl_dir
self.right_dist =
self.Pa.dist(p2)*xr_dir
'''
Calculate the unit vector of PaPb
'''
self.uv = unitV(Pb-self.Pa)
else:
raise TypeError, "All
arguments passed to DistancePointLine3D must be <class
'__main__.Point'>"
def __repr__(self):
return 'DistancePointLine3D(%s, %s, %s)'
% (repr(self.p1), repr(self.p2), repr(self.Pa))
def __str__(self):
s1 = 'Instance of Class %s - Closest
point on Line Segment:' % (type(self).__name__)
s2 = '
Point on line segment p1-p2 (self.Pa): %s' % (repr(self.Pa))
s3 = '
Distance between closest point on line segment and point: %0.6f' %
(self.dist)
return '%s\n%s\n%s' % (s1, s2, s3)
def version(self):
"PointPlane3D.DistancePointLine3D
Version 2.00"
##
BasisTransToGlobal ###########################
class
BasisTransToGlobal(Plane3D):
'''
Given 3 counter-clockwise non-collinear
points, define an orthonormal basis
in 3D space (the local basis).
Calculate the global point, a vector in the
standard basis set, given a
displacement vector in the local basis.
The object type of the displacement vector
'vR' is '<class '__main__.Point'>'.
An instance translate method can be called:
ref_pt + instance.translate(x, y, z)
--> global coordinate point object
Method 'translate' example usage:
pt = global point origin of translation
to be calculated (example: mem.left.location)
a = class instance
x, y, and z are offsets from pt in the
local basis
pt + a.translate(x, y, z)
Y |
/<----local X
| |<--/--local Z
| pt.
| /
| | /
|
________|/
|
^local Y
|
|
Z._______________X
Interactive Example:
Point 'p1' is the local coordinate (0,0,0)
and 'p2' determines the local
basis 'X' axis. Point 'p3' is required to
define the plane of the local
basis. The 'Y' and 'Z' axes are calculated.
Point(2,2,2) is a coordinate
in the local basis. Instance attribute 'R'
is the global coordinate
of Point(2,2,2).
>>> p1
Point(5.000000, 5.000000, 5.000000)
>>> p2
Point(7.000000, 10.000000, 6.000000)
>>> p3
Point(1.000000, 15.000000, 7.000000)
>>>
a=BasisTransToGlobal(p1,p2,p3,Point(2,2,2))
>>> a.R
Point(3.868398, 7.149624, 7.469533)
>>> p1+a.translate(2,2,2)
Point(3.868398, 7.149624, 7.469533)
>>> p1+a.translate(6,6,6)
Point(1.605194, 11.448873, 12.408598)
>>>
'''
def __init__(self, vN, vA, vB,
vR=Point(0,0,0)):
'''
vN = local coordinate (0,0,0)
vA = determines local basis 'X' axis
vB = defines local basis 'XY' plane
Local basis 'Z' axis is Plane3D
instance 'N0', self.N
'Y' axis = self.B =
self.N.cross(self.A)
'''
self.vN = vN
self.vA = vA
self.vB = vB
self.vR = vR
'''
The local basis is aligned with plane
defined
'''
Plane3D.__init__(self, vN, vA, vB)
if isinstance(vR, Point):
if mag(self.N) > epsilon:
'''
Unit vector normal to defined
plane, local basis 'Z'
'''
self.N = self.N0
'''
Unit vector between vN and vA,
local basis 'X'
'''
self.A = self.d0
'''
Unit cross product vector,
local basis 'Y'
'''
self.B = self.N.cross(self.A)
'''
Calculate the global coordinate
vector 'R'
'vR' is the local coordinate
vector
'''
self.R =
self.translate(*vR.data) + vN
else:
Warning("The points are
collinear ***INVALID***")
self.R = None
else:
raise TypeError, "Arguments
must be <class '__main__.Point'>"
def translate(self, X, Y, Z):
def determinant3(a,b,c,m,n,k,u,v,w):
return a*n*w + b*k*u + m*v*c -
c*n*u - b*m*w - a*k*v
A, B, N = self.A, self.B, self.N
'''
Normalize vR
'''
M = (X*X + Y*Y + Z*Z)**0.5
try:
X1, Y1, Z1 = X/M, Y/M, Z/M
except:
X1, Y1, Z1 = 0.0, 0.0, 0.0
D = determinant3(A.x, A.y, A.z, N.x,
N.y, N.z, B.x, B.y, B.z)
Dx = determinant3(X1, A.y, A.z, Z1,
N.y, N.z, Y1, B.y, B.z)
Dy = determinant3(A.x, X1, A.z, N.x,
Z1, N.z, B.x, Y1, B.z)
Dz = determinant3(A.x, A.y, X1, N.x,
N.y, Z1, B.x, B.y, Y1)
'''
Calculate the resultant unit vector
'R1'
'''
R1 = Point(Dx/D, Dy/D, Dz/D)
'''
Return the global coordinate vector
'''
return R1*M
def __repr__(self):
return 'BasisTransToGlobal(%s, %s, %s,
%s)' % \
(repr(self.vN), repr(self.vA),
repr(self.vB), repr(self.vR))
def version(self):
"PointPlane3D.BasisTransToGlobal
Version 2.00"
## BasisTransToLocal
############################
class
BasisTransToLocal(Plane3D):
'''
Given 3 counter-clockwise non-collinear
points, define a local orthonormal
basis.
Given a fourth point in the standard basis
set (a global coordinate point
object), calculate the local coordinate in
the defined orthonormal basis.
Method 'trans_to_local' example usage:
a = class instance (defines local
basis)
x, y, and z are floats corresponding to
the global coordinate point
object attributes x, y, and z
a.trans_to_local(x, y, z) --> local
point in local basis defined
by class instance 'a'
Y |
/<----local X
| |<--/--local Z
| pt.
| /
| | /
|
________|/
|
^local Y
|
|
Z._______________X
Interactive Example:
Point 'p1' is the local coordinate (0,0,0)
and 'p2' determines the local
basis 'X' axis. Point 'p3' is required to
define the plane of the local
basis. The 'Y' and 'Z' axes are calculated.
Point(2,2,2) is a global
coordinate. Instance attribute 'P' is the
local coordinate of Point(2,2,2).
>>> p1
Point(5.000000, 5.000000, 5.000000)
>>> p2
Point(7.000000, 10.000000, 6.000000)
>>> p3
Point(1.000000, 15.000000, 7.000000)
>>> b =
BasisTransToLocal(p1,p2,p3,Point(2,2,2))
>>> b.P
Point(-4.381780, 1.503841, -2.353394)
>>> b.trans_to_local(3,3,3)
Point(-2.921187, 1.002561, -1.568929)
>>> b.trans_to_local(6,6,6)
Point(1.460593, -0.501280, 0.784465)
>>>
'''
def __init__(self, vN, vA, vB,
R1=Point(0,0,0)):
'''
The local basis is aligned with defined
plane
vN = local coordinate (0,0,0)
vA = determines local basis 'X' axis
vB = defines local basis 'XY' plane
Local basis 'Z' axis is Plane3D
instance 'N0', self.N
Unit Vectors:
'X' axis = self.A
'Y' axis = self.B =
self.N.cross(self.A)
'Z' axis = self.N
'''
Plane3D.__init__(self, vN, vA, vB)
self.vN = vN
self.vA = vA
self.vB = vB
self.R1 = R1
if isinstance(R1, Point):
if mag(self.N) > epsilon:
'''
Unit vector normal to defined
plane, local basis 'Z'
'''
self.N = self.N0
'''
Unit vector between vN and vA,
local basis 'X'
'''
self.A = self.d0
'''
Unit cross product vector,
local basis 'Y'
'''
self.B = self.N.cross(self.A)
'''
Instance attribute 'P' is the
local coordinate vector
'''
self.P =
self.trans_to_local(*R1.data)
else:
Warning("The points are
collinear ***INVALID***")
self.P = None
else:
raise TypeError, "Arguments
must be <class '__main__.Point'>"
def trans_to_local(self, X, Y, Z):
R = Point(X, Y, Z) - self.vN
'''
Vector projection of: R along vA, R
along vB, R along vN
'''
return Point(R.dot(self.A),
R.dot(self.B), R.dot(self.N))
def __repr__(self):
return 'BasisTransToLocal(%s, %s, %s,
%s)' % \
(repr(self.vN), repr(self.vA),
repr(self.vB), repr(self.R1))
def version(self):
"PointPlane3D.BasisTransToLocal
Version 2.00"
class
CircleCircleIntersection(object):
def __init__(self, p1, p2, r1, r2):
'''
Given circle center points p1 and p2
and circle radii r1 and r2
Calculate intersection points of
circles
Results are valid only when in the X-Y
plane at the present time
'''
self.p1 = p1
self.p2 = p2
self.r1 = r1
self.r2 = r2
self.d = p1.dist(p2)
if self.d > r1+r2:
self.Pa = None
self.Pb = None
elif self.d < abs(r1-r2):
self.Pa = None
self.Pb = None
else:
self.a =
(r1**2-r2**2+self.d**2)/(2*self.d)
self.b = self.d-self.a
self.P0 = polarPt(p1, p2, self.a,
p1)
self.h = (r1**2-self.a**2)**0.5
self.intPts()
def intPts(self):
self.Pa = Point()
self.Pb = Point()
self.Pa.x = self.P0.x + (self.h *
(self.p2.y - self.p1.y) / self.d)
self.Pb.x = self.P0.x - (self.h *
(self.p2.y - self.p1.y) / self.d)
self.Pa.y = self.P0.y - (self.h *
(self.p2.x - self.p1.x) / self.d)
self.Pb.y = self.P0.y + (self.h *
(self.p2.x - self.p1.x) / self.d)
self.theta0 =
math.atan2(self.p2.y+self.p1.y, self.p2.x+self.p1.x)
self.theta1 = self.theta0 -
math.atan2(self.h, self.a)
self.theta2 = self.theta0 +
math.atan2(self.h, self.a)
## Miscellaneous
functions ######################
def simplify(n, d):
'''
Return the simplified numerator and
denominator.
>>> simplify(144,256)
(9, 16)
'''
n, d = map(int, [n, d])
b = gcd(n, d)
return n/b, d/b
def gcd(m,n):
"""Return greatest common
divisor using Euclid's Algorithm."""
if n == 0: return m
return gcd(n,m%n)
def fifDim(d, a=16):
'''
Distance formatting function (float -->
feet, inches, fraction)
>>> fifDim(146.45, 16)
"12'-2 7/16"
>>> fifDim(146.45, 4)
"12'-2 1/2"
'''
if d < 0: sign = '-'
else: sign = ''
# if interger, convert to float
# round off to the nearest 'a'
feet, rem = divmod(abs(round(float(d) * a)
/ a), 12.0)
frac, inch = math.modf(rem)
num, den = simplify(frac*a, a)
if num > 0 and feet > 0:
return "%s%s'-%s %s/%s" %
(sign, int(feet), int(inch), num, den)
elif num > 0 and inch > 0:
return "%s%s %s/%s" % (sign,
int(inch), num, den)
elif feet > 0 and num == 0:
return "%s%s'-%s" % (sign,
int(feet), int(inch))
elif inch >= 0 and num == 0:
return "%s%s" % (sign,
int(inch))
else:
return "%s%s/%s" % (sign,
num, den)
def
cross_product(p1, p2):
'''
Point cross product calculation
The cross product of two vectors returns a
vector that is perpendicular to both vectors.
>>> cross_product(p1,p2)
Point(-20.000000, 5.000000, 15.000000)
'''
return Point(p1.y*p2.z - p1.z*p2.y,
p1.z*p2.x - p1.x*p2.z, p1.x*p2.y - p1.y*p2.x)
def dot_product(p1,
p2):
'''
Point dot product calculation
The dot product of two unit vectors returns
a scalar value equal to the cosine of
the angle between the vectors.
>>> p1
Point(5.000000, 5.000000, 5.000000)
>>> p2
Point(7.000000, 10.000000, 6.000000)
>>> dot_product(p1,p2)
115.0
>>> math.acos(dot_product(p1.uv(),
p2.uv()))
0.21816790834595173
>>> math.acos(dot_product(p1.uv(),
p2.uv())) * (180/math.pi)
12.500100373420004
>>>
'''
return (p1.x*p2.x + p1.y*p2.y + p1.z*p2.z)
def mag(p):
'''
Return the magnitude of a Point object vector
>>> mag(p2-p1)
5.4772255750516612
>>>
'''
return (p.x**2 + p.y**2 + p.z**2)**0.5
def unitV(p):
'''
Return the unit vector of a Point object
vector
>>> unitV(p2-p1)
Point(0.365148, 0.912871, 0.182574)
>>> unitV(p1)
Point(0.577350, 0.577350, 0.577350)
>>> unitV(p2)
Point(0.514650, 0.735215, 0.441129)
>>>
'''
m = mag(p)
try:
return Point(p.x/m, p.y/m, p.z/m)
except:
return Point(0,0,0)
def polarPt (p1, p2,
d, p3=False):
'''
Calculate the vector p1-->p2
Translate distance 'd' parallel to vector
p1-->p2 from point 'p3'
Return the new point
'p3' defaults to 'p2'
>>> polarPt(p1,p2,10,p3)
Point(4.651484, 24.128709, 8.825742)
>>> polarPt(p1,p2,10)
Point(10.651484, 19.128709, 7.825742)
>>>
'''
return (p3 or p2)+(p2-p1).uv()*d
def polarPtXY (p1,
d, a):
'''
Y+
|
/ positive angle
| /
|/______X+ ---> a = 0.0
Return a point in the X-Y plane at a given
distance and angle, in radians,
from a reference point
'''
return Point(p1.x+d*math.cos(a),
p1.y+d*math.sin(a), p1.z)
def polarPtXZ (p1,
d, a):
'''
Z+
| / positive angle
| /
|/______X+ ---> a = 0.0
Return a point in the X-Z plane at a given
distance and angle, in radians,
from a reference point
'''
return Point(p1.x+d*math.cos(a), p1.y,
p1.z+d*math.sin(a))
def polarPtRotate(p1,
p2, p, theta):
'''
Rotate point object p about an axis line
passing through p1 and p2 by
angle 'theta' in radians. Return the
rotated point object. A right hand
coordinate system is assumed. Positive
angles are counter-clockwise when
looking down the axis toward the origin.
>>> polarPtRotate(Point(0,0,0),
Point(10,0,0), Point(15,-3,0), math.pi)
Point(15.000000, 3.000000, -0.000000)
>>> polarPtRotate(Point(0,0,0),
Point(10,0,0), Point(15,-3,0), math.pi/2)
Point(15.000000, -0.000000, -3.000000)
>>> polarPtRotate(Point(0,0,0),
Point(10,0,0), Point(15,-3,0), math.pi/4)
Point(15.000000, -2.121320, -2.121320)
>>> polarPtRotate(Point(0,0,0),
Point(10,0,0), Point(15,-3,0), -math.pi/4)
Point(15.000000, -2.121320, 2.121320)
>>> polarPtRotate(Point(0,0,0),
Point(10,10,10), Point(15,-3,0), math.pi/3)
Point(11.000000, 8.000000, -7.000000)
>>> polarPtRotate(Point(10,10,0),
Point(10,10,10), Point(15,-3,0), math.pi/3)
Point(23.758330, 7.830127, 0.000000)
>>>
'''
# Initialize point q
q = Point(0.0,0.0,0.0)
# Rotation axis unit vector
n = (p2-p1).uv()
# Translate so axis is at origin
p = p - p1
# Matrix common factors
c = math.cos(theta)
t =
(1 - math.cos(theta))
s = math.sin(theta)
X = n.x
Y = n.y
Z = n.z
# Matrix 'M'
d11 = t*X**2 + c
d12 = t*X*Y - s*Z
d13 = t*X*Z + s*Y
d21 = t*X*Y + s*Z
d22 = t*Y**2 + c
d23 = t*Y*Z - s*X
d31 = t*X*Z - s*Y
d32 = t*Y*Z + s*X
d33 = t*Z**2 + c
#
|p.x|
# Matrix 'M'*|p.y|
#
|p.z|
q.x = d11*p.x + d12*p.y + d13*p.z
q.y = d21*p.x + d22*p.y + d23*p.z
q.z = d31*p.x + d32*p.y + d33*p.z
#
Translate axis and rotated point back to original location
return q + p1
def midPt(p1, p2):
'''
Return the midpoint between two Point
object vectors
'''
return (p2-p1)/2+p1
def spaPts(p1, p2,
spaces=2, ends=False):
'''
Return a list of evenly spaced points
between p1 and p2
If 'ends' is True, the point list includes
p1 and p2
>>> spaPts(p1,p2,5)
[Point(5.400000, 6.000000, 5.200000), \
Point(5.800000, 7.000000, 5.400000), \
Point(6.200000, 8.000000, 5.600000), \
Point(6.600000, 9.000000, 5.800000)]
>>>spaPts(p1,p2,5,1)
[Point(5.000000, 5.000000, 5.000000), \
Point(5.400000, 6.000000, 5.200000), \
Point(5.800000, 7.000000, 5.400000), \
Point(6.200000, 8.000000, 5.600000), \
Point(6.600000, 9.000000, 5.800000), \
Point(7.000000, 10.000000, 6.000000)]
>>>
>>> ptList = spaPts(Point(0,0,0),
Point(1000,1000,1000), 10, 1)
>>> print '\n'.join([repr(pt) for
pt in ptList])
Point(0.000000, 0.000000, 0.000000)
Point(100.000000, 100.000000, 100.000000)
Point(200.000000, 200.000000, 200.000000)
Point(300.000000, 300.000000, 300.000000)
Point(400.000000, 400.000000, 400.000000)
Point(500.000000, 500.000000, 500.000000)
Point(600.000000, 600.000000, 600.000000)
Point(700.000000, 700.000000, 700.000000)
Point(800.000000, 800.000000, 800.000000)
Point(900.000000, 900.000000, 900.000000)
Point(1000.000000, 1000.000000,
1000.000000)
>>>
'''
p0 = p2-p1
ptList = [p0*(i/float(spaces))+p1 for i in
xrange(1, spaces)]
if ends:
ptList.insert(0,p1)
ptList.append(p2)
return ptList
def fixSpaPts(p1,
p2, spacing, leftDist=False, rightDist=False, ends=False):
'''
Return a list of points at a fixed spacing
between p1 and p2
'leftDist' and 'rightDist' are minumum
distances from the
respective ends (p1=left end, p2=right end)
and default to
'spacing/2'.
If 'ends' is True, the point list includes
p1 and p2.
>>> fixSpaPts(p1,p2,3,0.5,1)
[Point(5.360990, 5.902476, 5.180495),
Point(6.456435, 8.641089, 5.728218)]
>>> Point(5.360990, 5.902476,
5.180495).dist(Point(6.456435, 8.641089, 5.728218))
3.0000002327538233
>>>
'''
spacing = float(spacing)
leftDist=leftDist or spacing/2
rightDist=rightDist or spacing/2
maxOutOut = p1.dist(p2)-leftDist-rightDist
if maxOutOut < 0.0:
if ends:
return [p1,p2]
return []
else:
noSpa =
int(math.floor(maxOutOut/spacing))
ptList =
[polarPt(p1,p2,(leftDist+(maxOutOut-spacing*noSpa)/2),p1), ]
for i in xrange(noSpa):
ptList.append(polarPt(p1,p2,spacing,ptList[-1]))
if ends:
ptList.insert(0,p1)
ptList.append(p2)
return ptList
def autoSpaPts(p1,
p2, spacing, leftDist=False, rightDist=False, ends=False):
'''
Return a list of points at a maximum
spacing between p1 and p2
'leftDist' and 'rightDist' are the
distances from the
respective ends (p1=left end, p2=right end)
and default to
'spacing/2'.
If 'ends' is True, the point list includes
p1 and p2.
>>> autoSpaPts(p1,p2,3,0.5,1)
[Point(5.182574, 5.456435, 5.091287), \
Point(5.908713, 7.271782, 5.454356), \
Point(6.634852, 9.087129, 5.817426)]
>>> Point(5.182574, 5.456435,
5.091287).dist(Point(5.908713, 7.271782, 5.454356))
1.9886130031987117
>>>
'''
spacing = float(spacing)
leftDist=leftDist or spacing/2
rightDist=rightDist or spacing/2
OutOut = p1.dist(p2)-leftDist-rightDist
if OutOut < 0.0:
if ends:
return [p1,p2]
return []
else:
noSpa = int(math.ceil(OutOut/spacing))
spacing = OutOut/noSpa
ptList = [polarPt(p1,p2,leftDist,p1), ]
if noSpa > 0:
for i in xrange(noSpa):
ptList.append(polarPt(p1,p2,spacing,ptList[-1]))
else:
ptList =
ptList.append(polarPt(p2,p1,rightDist,p2))
if ends:
ptList.insert(0,p1)
ptList.append(p2)
return ptList
def polarSpaPts(a,
spaces=6, ends=False):
'''
'a' is a Plane3D instance.
Return a list of evenly spaced points
starting at point a.p2, rotated
about point a.p1, between the angle a.Q and
along radius a.Ra.
If 'ends' is True, the point list includes
the start and end points.
Interactive example:
>>> ptList =
polarSpaPts(Plane3D(Point(), Point(10,0,0), Point(0,10,0), math.pi/2), 4, True)
>>> for p in ptList: print p
...
(10.0000, 0.0000, 0.0000)
(9.2388, 3.8268, 0.0000)
(7.0711, 7.0711, 0.0000)
(3.8268, 9.2388, 0.0000)
(0.0000, 10.0000, 0.0000)
>>>
'''
try:
spacing = a.pp/spaces
ptList = [a.PointRotate3D(a.p2,
spacing*i/a.Ra) for i in range(1,spaces)]
if ends: return [a.p2,] + ptList +
[a.PointRotate3D(a.p2, a.Q),]
else: return ptList
except:
raise TypeError, 'Invalid argument
passed to polarSpaPts'
def polarFixSpaPts(a,
spacing, startDist=False, endDist=False, ends=False):
'''
'a' is a Plane3D instance.
Return a list of points at a fixed spacing
starting at point a.p2, rotated
about point a.p1, between the angle a.Q and
along radius a.Ra.
'startDist' and 'endDist' are minumum
distances along the arc from the
respective ends (a.p2 is the start point,
the array progresses
counter-clockwise) and default to
'spacing/2'.
If 'ends' is True, the point list includes
the start and end points.
Interactive example:
>>> ptList =
polarFixSpaPts(Plane3D(Point(), Point(10,0,0), Point(0,10,0), math.pi/2), 4,
0.5, 2, True)
>>> for p in ptList: print p
...
(10.0000, 0.0000, 0.0000)
(9.9391, 1.1017, 0.0000)
(8.7255, 4.8852, 0.0000)
(6.1343, 7.8975, 0.0000)
(2.5747, 9.6629, 0.0000)
(0.0000, 10.0000, 0.0000)
>>> Point(8.7255, 4.8852,
0.0000).dist(Point(6.1343, 7.8975, 0.0000))
3.973445448222487
>>>
'''
startDist=startDist or spacing/2
endDist=endDist or spacing/2
maxOutOut = a.pp - startDist - endDist
spaces = int(math.floor(maxOutOut/spacing))
netDist = startDist +
((maxOutOut-(spaces*spacing))/2.0)
ptList = [a.PointRotate3D(a.p2,
(netDist+spacing*i)/a.Ra) for i in range(0,spaces+1)]
if ends: return [a.p2,] + ptList +
[a.PointRotate3D(a.p2, a.Q),]
else: return ptList
def
polarAutoSpaPts(a, spacing, startDist=False, endDist=False, ends=False):
'''
'a' is a Plane3D instance.
Return a list of points at a maximum
spacing starting at point a.p2,
rotated about point a.p1, between the angle
a.Q and along radius a.Ra.
'startDist' and 'endDist' are exact
distances along the arc from the
respective ends (a.p2 is the start point,
the array progresses
counter-clockwise) and default to
'spacing/2'.
If 'ends' is True, the point list includes
the start and end points.
Interactive example:
>>> ptList =
polarAutoSpaPts(Plane3D(Point(), Point(10,0,0), Point(0,10,0), math.pi/2), 4,
0.5, 2, True)
>>> for p in ptList: print p
...
(10.0000, 0.0000, 0.0000)
(9.9875, 0.4998, 0.0000)
(9.2859, 3.7111, 0.0000)
(7.5810, 6.5214, 0.0000)
(5.0571, 8.6271, 0.0000)
(1.9867, 9.8007, 0.0000)
(0.0000, 10.0000, 0.0000)
>>> Point(10.0000, 0.0000,
0.0000).dist(Point(9.9875, 0.4998, 0.0000))
0.49995628808926884
>>> Point(1.9867, 9.8007,
0.0000).dist(Point(0.0000, 10.0000, 0.0000))
1.996671575397416
>>> Point(7.5810, 6.5214,
0.0000).dist(Point(5.0571, 8.6271, 0.0000))
3.2869505168164617
>>>
'''
startDist=startDist or spacing/2
endDist=endDist or spacing/2
spaces = int(math.ceil((a.pp - startDist -
endDist) / spacing))
actualSpa = (a.pp - startDist - endDist) /
spaces
ptList = [a.PointRotate3D(a.p2,
(startDist+actualSpa*i)/a.Ra) for i in range(0,spaces+1)]
if ends: return [a.p2,] + ptList +
[a.PointRotate3D(a.p2, a.Q),]
else: return ptList
def polarEllipse(a,
res, quad=1, xAxis=False, yAxis=False, rot=0):
'''
Return a list of points that represent an
ellipse
'a' is a Plane3D instance.
'res' is the number of resolution points on
the ellipse
'quad' is the quadrant of the local basis.
'rot' is an optional rotation of the
ellipse with respect to the local
basis in radians.
One quadrant is laid out at a time.
Quadrants are 1, 2, 3, or 4
("0-90", "90-180", "180-270",
"270-360").
The center of the ellipse is at the origin
of the Plane3D instance (a.p1).
The semi-major axis (xAxis) defaults to
a.Ra (a.p1.dist(a.p2)).
The semi-minor axis (yAxis) defaults to
a.p1.dist(a.p3).
.(a.p3)
|
|
|<--semi-minor axis (h)
|
|
(a.p1)._____________________. (a.p2)
semi-major axis (b)
Equation of an ellipse: (x**2/b**2) +
(y**2/h**2) = 1
Equation of resolution ray: mx-y = 0 where
m = tan(A)
Substituting and solving for 'x':
x=(b**2*h**2/(h**2+b*m**2))**0.5
Solving for 'y':
y=(((b**2*h**2)-(h**2*x**2))/b**2)**0.5
Interactive examples:
>>> a = Plane3D(Point(5,5,5),
Point(105,5,5), Point(5,105,5))
>>> ptList = polarEllipse(a, 10)
>>> for pt in ptList: print pt
...
(105.0000, 5.0000, 5.0000)
(103.7688, 20.6434, 5.0000)
(100.1057, 35.9017, 5.0000)
(94.1007, 50.3990, 5.0000)
(85.9017, 63.7785, 5.0000)
(75.7107, 75.7107, 5.0000)
(63.7785, 85.9017, 5.0000)
(50.3990, 94.1007, 5.0000)
(35.9017, 100.1057, 5.0000)
(20.6434, 103.7688, 5.0000)
(5.0000, 105.0000, 5.0000)
>>>
>>> a = Plane3D(Point(5,5,5),
Point(105,5,25), Point(5,75,50))
>>> ptList = polarEllipse(a, 10,
quad=1, rot=1.5707963)
>>> for pt in ptList: print
repr(pt)
...
Point(-5.665271, 91.270216, 58.326370)
Point(-21.018146, 89.689677, 54.239735)
Point(-35.095564, 85.150003, 48.505889)
Point(-47.265563, 78.176263, 41.588770)
Point(-57.228173, 69.431933, 33.974894)
Point(-64.965609, 59.542390, 26.069843)
Point(-70.630978, 48.995608, 18.156695)
Point(-74.441970, 38.120370, 10.403272)
Point(-76.608608, 27.109452, 2.891497)
Point(-77.294520, 16.057820, -4.350306)
Point(-76.600575, 5.000002, -11.320114)
>>> a.p1.dist(Point(-5.665271,
91.270216, 58.326370))
101.98038983827723
>>> a.p1.dist(Point(-76.600575,
5.000002, -11.320114))
83.216584652962212
>>>
>>> a = Plane3D(Point(5,5,5),
Point(105,5,25), Point(5,75,50))
>>> ptList = polarEllipse(a, 10,
3, 90.0, 45.0, rot=0)
>>> for pt in ptList: print repr(pt)
...
Point(-83.252261, 5.000000, -12.650452)
Point(-77.710960, -6.495697, -18.932283)
Point(-66.435611, -15.742853, -22.621814)
Point(-53.453189, -22.170711, -24.157524)
Point(-41.154722, -26.359268, -24.390474)
Point(-30.258282, -29.048793, -23.940166)
Point(-20.709500, -30.779975, -23.143313)
Point(-12.226952, -31.889392, -22.159999)
Point(-4.506638, -32.575074, -21.056732)
Point(2.724405, -32.948898, -19.850839)
Point(9.706173, -33.067708, -18.530864)
>>> a.p1.dist(Point(-83.252261,
5.000000, -12.650452))
90.000000152313476
>>> a.p1.dist(Point(9.706173,
-33.067708, -18.530864))
45.000000191840989
>>>
'''
localRP_list = []
xAxis = float(xAxis or a.Ra)
yAxis = float(yAxis or a.p1.dist(a.p3))
for i in range(res + 1):
xOff =
(xAxis**2*yAxis**2/(yAxis**2+(xAxis*math.tan((math.pi/2.0)/res *
(i)))**2))**0.5
yOff =
(((xAxis**2*yAxis**2)-(yAxis**2*xOff**2))/xAxis**2)**0.5
localRP_list.append(Point(xOff, yOff,
0.0))
# Rotate points into quadrant 1-4 + rot
argument
theta = {1:0.0, 2:math.pi/2.0, 3:math.pi,
4:math.pi*1.5}[quad]+rot
b = BasisTransToGlobal(a.p1, a.p2, a.p3)
return [a.PointRotate3D(a.p1+b.translate(*p),
theta) for p in localRP_list]
def sortPts(ptlist,
key='x', reverse=False):
'''
Sort a point list on the 'x', 'y' or 'z'
attribute
'''
try:
def cmpItems(a,b):
if reverse:
return -(cmp(getattr(a,
key.lower()), getattr(b, key.lower())))
return cmp(getattr(a, key.lower()),
getattr(b, key.lower()))
ptlist.sort(cmpItems)
except AttributeError, e:
raise AttributeError, 'Error in sortPoints()
- %s' % (e)
def
angBetweenLines(p1,p2,p3,p4):
'''
<--> <-->
Return the true angle between lines p1p2
and p3p4 in radians.
>>> p1=Point(10,10,0)
>>> p2=Point(20,20,0)
>>> p3=Point(60,40,0)
>>> p4=Point(20,20,20)
>>>
angBetweenLines(p1,p2,p3,p4)*180.0/math.pi
150.00000000000003
>>>
angBetweenLines(p1,p2,p4,p3)*180.0/math.pi
30.000000000000011
>>>
SDS/2 members:
angBetweenLines(*(memWPs(mem1)+memWPs(mem2)))
'''
return math.acos((p2-p1).uv().dot((p4-p3).uv()))
def
angBetweenLinesProjZ(a,p3,p4):
'''
<-->
Given a BasisTransToLocal instance 'a' and
a line p3p4, return the
complement to the projected angle in the
'Z' plane of the local basis
defined by instance 'a', in radians.
Interactive example:
>>>
angBetweenLinesProjZ(BasisTransToLocal(Point(5,5,0), Point(10,10,0),
Point(5,5,10)), Point(6,6,8), Point(10,5,12))
0.54041950027058416
>>>
SDS/2 example:
from member import MemberLocate
mem1 = MemberLocate('Select a member')
mem2 = MemberLocate('Select another
member')
print angBetweenLinesProjZ(memPlane(mem1),
Point(75,175,225), Point(175,225,905))
print angBetweenLinesProjZ(memPlane(mem2),
Point(75,175,225), Point(175,225,905))
print angBetweenLinesProjZ(memPlane(mem1),
*memWPs(mem2))
Returned:
-1.03780930466
-1.32178325868
-0.0184012714524
'''
pt1 = a.trans_to_local(*p3)
pt2 = a.trans_to_local(*p4)
try: return math.atan((pt2.x - pt1.x) /
(pt2.z - pt1.z))
# if the local 'z' coordinate difference is
0, the complement is 0 radians.
except: return 0
def endWP(m1, m2):
'''
Return the WP on member 1 and the framing
end of member 2 (left or right)
Arguments can be SDS/2 'member' type or a
tuple of Point objects
Interactive example:
>>> p1
Point(5.000000, 5.000000, 5.000000)
>>> p2
Point(7.000000, 10.000000, 6.000000)
>>> p3
Point(1.000000, 15.000000, 7.000000)
>>> p4
Point(10.000000, 10.000000, 10.000000)
>>> endWP((p1,p2),(p3,p4))
(Point(7.760629, 11.901573, 6.380315),
'right end')
>>>
'''
try:
a = LineLineIntersect3D(*(memWPs(m1)+memWPs(m2)))
b =
LineLineIntersect3D(*(memWPs(m2)+memWPs(m1)))
except:
a = LineLineIntersect3D(*(m1+m2))
b = LineLineIntersect3D(*(m2+m1))
if b.ma <= 0.5:
return (a.Pa, 'left end')
return (a.Pa, 'right end')
def memWPs(m):
'''
Return a tuple of SDS/2 member work points
as Point objects
'''
try:
return
Point(*eval(repr(m.left.location))), Point(*eval(repr(m.right.location)))
except:
raise TypeError, 'Incorrect argument
type passed to memWPs'
def memPlane(mem,
XZ=False):
'''
Create a BasisTransToLocal instance given
an SDS/2 member object
The local basis equates to the local member
coordinate system
The X-Y plane is at member web CL for W
flange, W Tee, etc.
If 'XZ' is True, the Y and Z planes swap.
The X-Y plane of the
local basis will be in the top flange plane
of a beam member.
Instance attribute 'A' represents 'X' axis
unit vector
Instance attribute 'B' represents 'Y' axis
unit vector
Instance attribute 'N' represents 'Z' axis
unit vector
Instance attributes 'p1' and 'vN' represent
the local basis (0,0,0)
Return the BasisTransToLocal instance
'''
p = memWPs(mem)
if XZ: trans = (0.0, 0.0, -1.0)
else: trans = (0.0, 1.0, 0.0)
return BasisTransToLocal(p[0], p[1],
p[0]+eval(repr(mem.translate(*trans))))
if __name__ ==
'__main__':
'''
def SDSangBetweenLinesProjZ():
from member import MemberLocate
mem1 = MemberLocate('Select a member')
mem2 = MemberLocate('Select another
member')
print
angBetweenLinesProjZ(memPlane(mem1), Point(75,175,225), Point(175,225,905))
print
angBetweenLinesProjZ(memPlane(mem2), Point(75,175,225), Point(175,225,905))
print
angBetweenLinesProjZ(memPlane(mem1), *memWPs(mem2))
SDSangBetweenLinesProjZ()
'''
'''
-1.03780930466
-1.32178325868
-0.0184012714524
'''
'''
def SDSmemPlane():
from member import MemberLocate
mem1 = MemberLocate('Select a member')
a = memPlane(mem1)
b = memPlane(mem1, True)
print
a.trans_to_local(*Point(360,240,1308))
print b.trans_to_local(*Point(360,240,1308))
SDSmemPlane()
'''
'''
def SDSPtTest():
def p(*pt):
return pt
from cons_line import ConsLine
cl2 = ConsLine()
cl2.pt1 = Point(10,10,10).data
cl2.angle = 0.0
cl2.pen = 'Red'
cl2.add()
cl2 = ConsLine()
cl2.pt1 = Point(10,10,10).data
cl2.angle = 90.0
cl2.pen = 'Yellow'
cl2.add()
cl2 = ConsLine()
cl2.pt1 = p(*Point(15,15,10))
cl2.angle = 0.0
cl2.pen = 'Cyan'
cl2.add()
cl2 = ConsLine()
cl2.pt1 = tuple(Point(15,15,10))
cl2.angle = 90.0
cl2.pen = 'Green'
cl2.add()
SDSPtTest()
'''
'''
# SDS/2 test of memPlane() and
angBetweenLinesProjZ()
def test_memPlane():
from member import MemberLocate
mem1 = MemberLocate('Select a member')
mem2 = MemberLocate('Select another
member')
print memWPs(mem1)[0], memWPs(mem1)[1]
print memWPs(mem2)[0], memWPs(mem2)[1]
return
angBetweenLinesProjZ(memPlane(mem1), *memWPs(mem2))
print test_memPlane()
'''
'''
# SDS/2 test of angBetweenLines()
def test_angBetweenLines():
from member import MemberLocate
mem1 = MemberLocate('Select a member')
mem2 = MemberLocate('Select another
member')
return
angBetweenLines(*(memWPs(mem1)+memWPs(mem2)))
print test_angBetweenLines()*180/math.pi
'''
p1=Point(1,2,3)
p2=Point(3,5,4)
p3=Point(6,8,5)
p4=Point(20,20,20)
a=Plane3D(p1,p2,p3)
a.ThreePtCircle()
print a.M
print a.R
print a.PointRotate3D(p2, math.pi)
print
print a.DistToPlane(Point(1000,1000,1000))
print
print a
print
print repr(a)
print
print a.data1()
print
print a.data2()
print
print a.data3()
print
print
p1=Point(144,256,300)
p2=Point(324,587,398)
p3=Point(645,400,130)
c=Plane3D(p1,p2,p3)
c.ThreePtCircle()
print c.M
print c.R
print c.PointRotate3D(p2, math.pi)
print
print c.DistToPlane(Point(1000,1000,1000))
print
print c
print
print repr(c)
print
print c.data1()
print
print c.data2()
print
print c.data3()
print
p1=Point(5,5,5)
p2=Point(7,10,6)
p3=Point(1,15,7)
b=Plane3D(p1,p2,p3)
b.ThreePtCircle()
print b.M
print b.R
print b.PointRotate3D(p2, math.pi)
print
print b.DistToPlane(Point(1000,1000,1000))
print
print b
print
print repr(b)
print
print b.data1()
print
print b.data2()
print
print b.data3()
print
print
d = Plane3D(Point(), Point(10,0,0),
Point(0,10,0))
print d
print d.N0
e = Plane3D(Point(), Point(0,10,0),
Point(10,0,0))
print e.N0
'''
# Cannot pickle bool objects in Python 2.3
# Works fine in 2.4
# probably fails on
self.__dict__['_Point__init']
print
print 'pickle test of Plane3D instance'
import pickle
fn = r'H:\TEMP\temsys\Plane3D_test.txt'
f = open(fn, "w")
pickle.dump(a, f)
f.close()
z = pickle.load(open(fn))
# z = Plane3D(dd[p1], dd[p2], dd[p3],
dd[Q])
print repr(z)
'''
'''
>>>
(19.5000, -7.0000, 0.0000)
20.7906228863
(-1.0000, -1.0000,
2.0000)
630.241937672
<__main__.Plane3D
object>
Plane3D(Point(1.000000,
2.000000, 3.000000), Point(3.000000, 5.000000, 4.000000), Point(6.000000,
8.000000, 5.000000), theta=0.105021)
Instance of Class
Plane3D - Plane Definition Data
Normal Unit Vector: 0.0000, 0.3162, -0.9487
Displacement from Model Origin: -2.2136
Plane Origin and
Unit Vector Points
Origin point: Point(1, 2, 3)
Unit Vector d0: 0.5345, 0.8018, 0.2673
Unit Vector e0: 0.6202, 0.7442, 0.2481
Plane3D(Point(1, 2,
3),
Point(3, 5, 4),
Point(6, 8, 5),
Sector (radians) = 0.105021)
Three Point Circle
Center point M: Point(1'-7 1/2, -7, 0)
Radius R: 1'-8 13/16
Arc length pp: 3/8
Sector Q: 0.105021
(391.2541, 368.3665,
239.6270)
278.218851749
(-36.0000, -75.0000,
202.0000)
562.739755277
<__main__.Plane3D
object>
Plane3D(Point(144.000000,
256.000000, 300.000000), Point(324.000000, 587.000000, 398.000000),
Point(645.000000, 400.000000, 130.000000), theta=0.967078)
Instance of Class
Plane3D - Plane Definition Data
Normal Unit Vector: -0.4005, 0.4535, -0.7962
Displacement from Model Origin: -180.4240
Plane Origin and
Unit Vector Points
Origin point: Point(12'-0, 21'-4, 25'-0)
Unit Vector d0: 0.4624, 0.8502, 0.2517
Unit Vector e0: 0.9137, 0.2626, -0.3100
Plane3D(Point(12'-0,
21'-4, 25'-0),
Point(27'-0, 48'-11, 33'-2),
Point(53'-9, 33'-4, 10'-10),
Sector (radians) = 0.967078)
Three Point Circle
Center point M: Point(32'-7 1/4, 30'-8 3/8,
19'-11 5/8)
Radius R: 23'-2 1/4
Arc length pp: 31'-4 1/2
Sector Q: 0.967078
(1.2500, 9.3269,
5.8654)
5.79082497112
(3.0000, 0.0000,
4.0000)
-780.54221785
<__main__.Plane3D
object>
Plane3D(Point(5.000000,
5.000000, 5.000000), Point(7.000000, 10.000000, 6.000000), Point(1.000000,
15.000000, 7.000000), theta=0.747584)
Instance of Class
Plane3D - Plane Definition Data
Normal Unit Vector: 0.0000, -0.1961, 0.9806
Displacement from Model Origin: 3.9223
Plane Origin and
Unit Vector Points
Origin point: Point(5, 5, 5)
Unit Vector d0: 0.3651, 0.9129, 0.1826
Unit Vector e0: -0.3651, 0.9129, 0.1826
Plane3D(Point(5, 5,
5),
Point(7, 10, 6),
Point(1, 1'-3, 7),
Sector (radians) = 0.747584)
Three Point Circle
Center point M: Point(1 1/4, 9 5/16, 5 7/8)
Radius R: 5 13/16
Arc length pp: 4 1/8
Sector Q: 0.747584
<__main__.Plane3D
object>
(0.0000, 0.0000,
1.0000)
(0.0000, 0.0000,
-1.0000)
>>> p1
Point(5.000000,
5.000000, 0.000000)
>>> p1.x=7
>>> p1.y=8
>>> p1.z=15
>>> p1
Point(7.000000,
8.000000, 15.000000)
>>> a
Plane3D(Point(1.000000,
2.000000, 3.000000), Point(3.000000, 5.000000, 4.000000), Point(6.000000,
8.000000, 5.000000), theta=0.105021)
>>>
a.DistToPlane(Point(300,500,250))
76.843347142091616
>>> a.N0
Point(0.000000,
0.316228, -0.948683)
>>> b
Plane3D(Point(5.000000,
5.000000, 5.000000), Point(7.000000, 10.000000, 6.000000), Point(1.000000,
15.000000, 7.000000), theta=0.747584)
>>>
b.DistToPlane(Point(300,500,250))
-143.16477865087435
>>> b.N0
Point(0.000000,
-0.196116, 0.980581)
>>> c
Plane3D(Point(144.000000,
256.000000, 300.000000), Point(324.000000, 587.000000, 398.000000),
Point(645.000000, 400.000000, 130.000000), theta=0.967078)
>>>
c.DistToPlane(Point(300,500,250))
-87.989570121985452
>>> c.N0
Point(-0.400516,
0.453529, -0.796177)
>>>
>>> p1 +=
p2
>>> p1
Point(14.000000,
18.000000, 15.000000)
>>> p1**2
Point(196.000000,
324.000000, 225.000000)
>>>
p1.mag()
27.294688127912362
>>>
p1.dist(p2)
18.384776310850235
>>>
p1.dot(p2)
278.0
>>>
p1.cross(p1)
Point(0.000000,
0.000000, 0.000000)
>>>
p1.cross(p2)
Point(-150.000000,
105.000000, 14.000000)
>>> p1.uv()
Point(0.512920,
0.659469, 0.549557)
>>>
p1.uv().dot(p2.uv())
0.83439852046300866
>>> import
math
>>>
math.acos(p1.uv().dot(p2.uv()))
0.58375573905810851
>>> b =
eval(repr(b))
>>>
math.acos(p3.uv().dot(p2.uv()))
0.54415780061338503
>>>
p1=Point(7.000000, 8.000000, 15.000000)
>>>
math.acos((p3-p1).uv().dot((p2-p1).uv()))
0.4578290164982331
>>> b
Plane3D(Point(7.000000,
8.000000, 15.000000), Point(7.000000, 10.000000, 0.000000), Point(1.000000,
15.000000, 0.000000), theta=0.457829)
>>>
repr(p1)
'Point(7.000000,
8.000000, 15.000000)'
>>> repr(a)
'Plane3D(Point(1.000000,
2.000000, 3.000000), Point(3.000000, 5.000000, 4.000000), Point(6.000000,
8.000000, 5.000000), theta=0.105021)'
>>> print a
<__main__.Plane3D
object>
>>> p1*10
Point(70.000000,
80.000000, 150.000000)
>>> -p1
Point(-7.000000,
-8.000000, -15.000000)
>>> [r for
r in p1]
[7.0, 8.0, 15.0]
>>>
p1+(1,2,3)
Point(8.000000,
10.000000, 18.000000)
>>> p1[1]
8.0
>>>
a.DistToPlane(Point(1000,1000,1000))
630.24193767155793
>>>
>>> p1
Point(5.000000,
5.000000, 0.000000)
>>> p2
Point(7.000000,
10.000000, 0.000000)
>>> p11
Point(5.000000,
5.000000, 1.000000)
>>>
p1<p11
True
>>>
p2>p11
True
>>> p2==p11
False
>>>
'''
''' Interactive
using LineLineIntersect3D
>>> m =
LineLineIntersect3D(Point(14,9,10), Point(20,1,20), Point(17,11,15),
Point(10,9,15))
>>> m.Pb
Point(14.426844,
10.264813, 15.000000)
>>> m.Pa
Point(15.273277,
7.302297, 12.122128)
>>>
m.inters_dist
4.2160515520138286
>>> m.uv
Point(-0.200764,
0.702675, 0.682599)
>>>
m.left_dist
3.0011424400000002
>>>
m.right_dist
-11.140993180000001
>>> m.ma
0.21221281741233375
>>> m.mb
0.36759371221281739
>>>
m.not_parallel()
True
>>>
m.on_segment1
1
>>>
m.on_segment2
1
>>>
m.position
'Not Beyond LE'
>>> m.uv1
Point(0.424264,
-0.565685, 0.707107)
>>> m.uv2
Point(-0.961524,
-0.274721, 0.000000)
>>>
'''
''' Interactive with
DistancePointLine3D
>>> r =
DistancePointLine3D(m.p1, m.p2, m.Pb)
>>> r
DistancePointLine3D(Point(14.000000,
9.000000, 10.000000), Point(20.000000, 1.000000, 20.000000), Point(15.273277,
7.302297, 12.122128))
>>> r.dist
4.2160515520138295
>>>
m.inters_dist
4.2160515520138286
>>> r.Pa
Point(15.273277,
7.302297, 12.122128)
>>> m.Pa
Point(15.273277,
7.302297, 12.122128)
>>> r.uv
Point(-0.200764,
0.702675, 0.682599)
>>> m.uv
Point(-0.200764,
0.702675, 0.682599)
>>>
r.left_dist
3.0011424449392754
>>>
m.left_dist
3.0011424400000002
>>>
r.right_dist
-11.140993178791675
>>>
m.right_dist
-11.140993180000001
>>>
r.position
'Not Beyond LE'
>>>
r.on_segment1
1
>>>
p1.equiv(p1)
(0, 0, 0)
>>>
p2.equiv(p2)
(0, 0, 0)
>>>
p1.equiv(p2)
(-1, -1, -1)
>>> p1
Point(5.000000,
5.000000, 5.000000)
>>> p2
Point(7.000000,
10.000000, 6.000000)
>>>
p2.equiv(p1)
(1, 1, 1)
>>>
'''
'''
BasisTransToGlobal and BasisTransToLocal interactive
>>> b =
BasisTransToGlobal(p1,p2,p3,Point(2,2,8))
>>> b.R
Point(3.868398,
5.972928, 13.353017)
>>> p1
Point(5.000000,
5.000000, 5.000000)
>>> p2
Point(7.000000, 10.000000,
6.000000)
>>> p3
Point(1.000000,
15.000000, 7.000000)
>>> c =
BasisTransToLocal(p1,p2,p3,Point(10,5,6))
>>> c.P
Point(2.008316,
-4.583135, 0.980581)
>>>
>>> ptList
= polarEllipse(Plane3D(Point(12,16.3,99.6836), Point(72.64,45.86499, 122.4663),
Point(-14.2664,-4.3245, 345.36645)), 15, 1)
>>> for p
in ptList:
... print p.format()
...
Point(6'-0 5/8, 3'-9
7/8, 10'-2 7/16)
Point(5'-10 9/16,
3'-8 5/8, 10'-9 9/16)
Point(5'-8 3/8, 3'-7
5/16, 11'-4 3/4)
Point(5'-6 1/16,
3'-5 15/16, 12'-0 3/16)
Point(5'-3 1/2, 3'-4
3/8, 12'-8 1/16)
Point(5'-0 11/16,
3'-2 3/4, 13'-4 9/16)
Point(4'-9 7/16,
3'-0 13/16, 14'-1 15/16)
Point(4'-5 5/8,
2'-10 5/8, 15'-0 1/2)
Point(4'-1 1/16,
2'-7 15/16, 16'-0 3/4)
Point(3'-7 3/8, 2'-4
5/8, 17'-3 3/16)
Point(2'-11 15/16,
2'-0 7/16, 18'-8 11/16)
Point(2'-2, 1'-6
13/16, 20'-6 1/8)
Point(1'-0 5/16, 11
3/16, 22'-8 1/4)
Point(-6 11/16,
13/16, 25'-2)
Point(-2'-7 1/8,
-1'-0 3/16, 27'-4)
Point(-4'-8 1/8,
-2'-0 7/8, 27'-10 1/2)
>>>
'''
'''
Problem: Find the
angle or slope of an object (member, etc.) or vector between two points.
The following code
finds the angle with respect to the horizontal plane.
Point p1 and p2 are
the end points of a vertical brace.
Module PointPlane3D
is used for Point objects.
[code]>>>
p1
Point(5.000000,
5.000000, 5.000000)
>>> p2
Point(7.000000,
10.000000, 6.000000)
>>> p3 =
p2-p1
>>> p3
Point(2.000000,
5.000000, 1.000000)
>>> import
math
>>>
math.atan2(p3.z, (p3.x**2+p3.y**2)**0.5)
0.18360401027891859
>>>
math.atan2(p3.z, (p3.x**2+p3.y**2)**0.5)*180/math.pi
10.51973489065862
>>> p3=-p3
>>>
math.atan2(p3.z, (p3.x**2+p3.y**2)**0.5)*180/math.pi
-10.51973489065862
>>> p3 =
Point(2, -5, 6)
>>>
math.atan2(p3.z, (p3.x**2+p3.y**2)**0.5)*180/math.pi
48.091152118196739
>>> [/code]
Another method would
be to calculate the angle between the vector p3 and a reference vector using
the dot product and unit vector methods of the Point objects:
[code]>>>
p3
Point(-2.000000,
5.000000, 6.000000)
>>> p4 =
Point(-2, 5, 0)
>>>
math.acos(p3.uv().dot(p4.uv()))*180/math.pi
48.091152118196732
>>> p5 =
Point(0,0,1)
>>>
math.acos(p3.uv().dot(p5.uv()))*180/math.pi
41.908847881803261
>>> [/code]
Point object p4
represents horizontal, and p5 represents vertical.
Another method would
be to use the Point object methods uv() (unit vector) and dot()
(dot product) and a
Plane3D object (obj.N0) normal unit vector to calculate the cosine
of the angle between
the unit vectors:
[code]>>>
>>> p3
Point(-2.000000,
5.000000, 6.000000)
>>> vector
= p3.uv()
>>> a =
Plane3D(Point(0,0,0), Point(1,0,0), Point(0,1,0))
>>> a.N0
Point(0.000000,
0.000000, 1.000000)
>>>
vector.dot(a.N0)
0.74420840753525075
>>> math.acos(vector.dot(a.N0))
0.73144738125491826
>>>
math.acos(vector.dot(a.N0))*180/math.pi
41.908847881803261
>>>
math.acos(vector.dot(a.N0))*180/math.pi+math.atan2(p3.z,
(p3.x**2+p3.y**2)**0.5)*180/math.pi
90.0
>>> [/code]
In the example
above, I created a Plane3D object representing the horizontal plane
(3 points define a
plane - counter-clockwise results in a unit vector going up).
The plane normal
unit vector, represented by the attribute 'N0' is truly vertical.
In this case, the
Plane3D unit vector is pointing down (points are clockwise):
[code]>>> b
= Plane3D(Point(0,0,0), Point(0,1,0), Point(1,0,0))
>>>
math.acos(vector.dot(b.N0))*180/math.pi
138.09115211819676
>>> [/code]
'''
'''
import timeit
from PointPlane3D import
Point
def
magnitude1(self):
'''
Returns the scalar length of the vector
'''
# return (self.x**2 + self.y**2 +
self.z**2)**0.5
return sum([i**2 for i in self])**0.5
def
magnitude2(self):
'''
Returns the scalar length of the vector
'''
return (self[0]**2 + self[1]**2 +
self[2]**2)**0.5
def
magnitude3(self):
return sum((self)**2)**0.5
def add1(p1, p2):
return Point(*[a+b for a,b in zip(p1,p2)])
def add2(p1, p2):
return Point(p1.x+p2[0], p1.y+p2[1], p1.z+p2[2])
def dist1(p1, p2):
return sum((p1-p2)**2)**0.5
def dist2(p1, p2):
p = p1-p2
return (p.x**2 + p.y**2 + p.z**2)**0.5
def dist3(p1, p2):
return ((p1-p2).x**2 + (p1-p2).y**2 +
(p1-p2).z**2)**0.5
if __name__ ==
'__main__':
t = timeit.Timer("magnitude1(Point(123.0,456.0,789.0))",
"from __main__ import magnitude1, Point")
print t.timeit(10000)
t =
timeit.Timer("magnitude2(Point(123.0,456.0,789.0))", "from
__main__ import magnitude2, Point")
print t.timeit(10000)
t =
timeit.Timer("magnitude3(Point(123.0,456.0,789.0))", "from
__main__ import magnitude3, Point")
print t.timeit(10000)
t = timeit.Timer("add1(Point(1,2,3),
Point(4,5,6))", "from __main__ import add1, Point")
print t.timeit(10000)
t = timeit.Timer("add2(Point(1,2,3),
Point(4,5,6))", "from __main__ import add2, Point")
print t.timeit(10000)
t =
timeit.Timer("dist1(Point(1.3,2.5,3.6), Point(4.123,5.567,6.987))",
"from __main__ import dist1, Point")
print t.timeit(10000)
t = timeit.Timer("dist2(Point(1.3,2.5,3.6),
Point(4.123,5.567,6.987))", "from __main__ import dist2, Point")
print t.timeit(10000)
t =
timeit.Timer("dist3(Point(1.3,2.5,3.6), Point(4.123,5.567,6.987))",
"from __main__ import dist3, Point")
print t.timeit(10000)
'''
'''
>>>
0.313378655668
0.29080539566
0.522658885417
0.778659578242
0.705695276913
1.00513277526
0.738010836576
1.31580494169
>>>
0.316868205318
0.30270746701
0.525020638095
0.776005330285
0.706250654763
1.0033546163
0.745588056583
1.30352069886
>>>
'''