PointPlane3D

## PointPlane3D.py

## Version 1.00 (May 30, 2007)

## NOT FOR SALE. The software is provided "as is" without any warranty.

################################################################################

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]]:

'''

/ \\

e,e0 G. \\

/ \\

p1.-----.-----. p2

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.

/ <---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 a.DistToPlane(Point(1000,1000,1000))

print a

print repr(a)

print a.data1()

print a.data2()

print a.data3()

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 c.DistToPlane(Point(1000,1000,1000))

print c

print repr(c)

print c.data1()

print c.data2()

print c.data3()

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 b.DistToPlane(Point(1000,1000,1000))

print b

print repr(b)

print b.data1()

print b.data2()

print b.data3()

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 '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

>>> m.on_segment2

>>> 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

>>> 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

>>>

'''