## PointPlane3D.py

## Version 1.00 (May 30, 2007)

## Copyright (c) 2007 Bruce Vaughan, BV Detailing & Design, Inc.

## All rights reserved.

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

            '''

                           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

>>>

'''