MARTIAN HOMEWORK...

WHAT WE DID IN CLASS TODAY

TO DO: IF NOT ALREADY IN YOUR FOLDER

CUT AND PASTE THE SOURCE CODE BELOW TO AN OPEN FILE WINDOW IN IDLE

AND SAVE IT AS

stickworks.py

"""
stickworks.py

Some infrastructure for working with Vectors and Edges, including
an xyplotter generator and axes maker.

By Kirby Urner, Sept 13, 2006

Updated Sept 29, 2006:
make Edge color a class-level attribute
add funky derivative demo
refactor a bit

Updated for Martian Math, 2010
Brought over Qvector subclass from coords.py

Code:
http://www.4dsolutions.net/ocn/python/stickworks.py

For colorized source:
http://www.4dsolutions.net/cgi-bin/py2html.cgi?script=/ocn/python/stickworks.py

Some relevant discussion:
http://mail.python.org/pipermail/edu-sig/2006-September/007145.html
http://mail.python.org/pipermail/edu-sig/2006-September/007149.html
http://mail.python.org/pipermail/edu-sig/2006-September/007150.html
http://mail.python.org/pipermail/edu-sig/2006-September/007312.html
"""

from visual import vector, cylinder, cross, dot, diff_angle, color
import visual
from operator import add, sub, mul, neg

root2   = 2.0**0.5

class Vector (object):

    """
    A wrapper for visual.vector that expresses a cylinder via draw(),
    always pegged to the origin
    """

    radius = 0.03
    
    def __init__(self, xyz, color=(0,0,1)):
        self.v = vector(*xyz)
        self.xyz = xyz
        self.color = color
        self.cyl = None

    def draw(self):
        """define and render the cylinder"""        
        self.cyl = cylinder(pos = (0,0,0), axis = self.v, radius = self.radius, color = self.color)

    def erase(self):
        """toss the cylinder"""        
        if self.cyl:
            self.cyl.visible = 0
        self.cyl = None

    def __repr__(self):
        return 'Vector @ (%s,%s,%s)' % self.xyz

    # some vector ops, including scalar multiplication
    
    def diff_angle(self, other):
        return self.v.diff_angle(other.v)

    def cross(self, other):
        temp = cross(self.v, other.v)
        return Vector((temp.x, temp.y, temp.z))
    
    def dot(self, other):
        return dot(self.v, other.v)

    def __sub__(self, other):
        temp = self.v - other.v
        return Vector((temp.x, temp.y, temp.z))

    def __add__(self, other):
        temp = self.v + other.v
        return Vector((temp.x, temp.y, temp.z))

    def __mul__(self, scalar):
        temp = self.v * scalar
        return Vector((temp.x, temp.y, temp.z))

    __rmul__ = __mul__

    def __neg__(self):
        return Vector((-self.v.x, -self.v.y, -self.v.z))

    def _length(self):
        return pow(self.v.x ** 2 + self.v.y ** 2 + self.v.z ** 2, 0.5)

    length = property(_length)

    def spherical(self):
        """Return (r,phi,theta) spherical coords based on current (x,y,z)"""
        r = self.length()

        if self.xyz[0]==0:
            if   self.xyz[1]==0: theta =   0.0
            elif self.xyz[1]< 0: theta = 270.0
            else:                theta =  90.0
            
        else:            
            theta = math.atan(self.xyz[1]/self.xyz[0]) * rad2deg
            if   self.xyz[0]<0 and self.xyz[1]==0:  theta =    180
            elif self.xyz[0]<0:                     theta =    180 + theta
            elif self.xyz[0]>0 and self.xyz[1]<0:   theta =    360 + theta

        if r==0: phi=0.0
        else: phi = math.acos(self.xyz[2]/r) * rad2deg
        
        return (r,phi,theta)

    def quadray(self):
        """return (a,b,c,d) quadray based on current (x,y,z)"""
        x=self.xyz[0]
        y=self.xyz[1]
        z=self.xyz[2]
        a = (2./root2) * ((x>=0)*x   + (y>=0)*y   + (z>=0)*z)
        b = (2./root2) * ((x<0)*(-x) + (y<0)*(-y) + (z>=0)*z)
        c = (2./root2) * ((x<0)*(-x) + (y>=0)*y   + (z<0)*(-z))
        d = (2./root2) * ((x>=0)*x   + (y<0)*(-y) + (z<0)*(-z))
        return self.norm((a,b,c,d))

    def norm(self,plist):
        """Normalize such that 4-tuple all non-negative members."""
        return tuple(map(sub,plist,[min(plist)]*4)) 
    
    def norm0(self):
        """Normalize such that sum of 4-tuple members = 0"""
        q = self.quadray()
        return tuple(map(sub,q,[reduce(add,q)/4.0]*4))
    
class Qray(Vector):
    """Subclass of Vector that takes quadray coordinate args"""
    
    def __init__(self, arg, color=(0,0,1),*flag):
        """Initialize a vector at an (a,b,c,d) tuple (= arg).
        
        NOTE: in accompanying essay, xyz units = sphere diameter
        i.e. Vector((1,0,0)).length() is 1 D, therefore quadray
        inputs must be scaled by 1/2 to fit this context, i.e.
        tetra edge defined by 2 basis quadrays = 1 D."""

        if len(arg)==3:  arg = Vector(arg).quadray() # if 3-tuple passed
            
        self.coords = self.norm(arg)

        a,b,c,d     =  self.coords
        self.xyz    = ((0.5/root2) * (a - b - c + d),
                       (0.5/root2) * (a - b + c - d),
                       (0.5/root2) * (a + b - c - d))
        self.v = vector(self.xyz)
        self.color = color
        self.cyl = None

    def __repr__(self):
        return "Qvector " + str(self.coords)
                    
    def dot(self,v1):
        """Return the dot product of self with another vector.

        return a scalar"""
        scalar = 0
        return 0.5*reduce(add,map(mul,self.norm0(),v1.norm0()))

    def cross(self,v1):
        """Return the cross product of self with another vector.
        return a Qvector"""
        
        A=Qray((1,0,0,0))
        B=Qray((0,1,0,0))
        C=Qray((0,0,1,0))
        D=Qray((0,0,0,1))
        a1,b1,c1,d1 = v1.quadray()
        a2,b2,c2,d2 = self.quadray()
        k= (2.0**0.5)/4.0
        thesum =   (A*c1*d2 - A*d1*c2 - A*b1*d2 + A*b1*c2
               + A*b2*d1 - A*b2*c1 - B*c1*d2 + B*d1*c2 
               + b1*C*d2 - b1*D*c2 - b2*C*d1 + b2*D*c1 
               + a1*B*d2 - a1*B*c2 - a1*C*d2 + a1*D*c2
               + a1*b2*C - a1*b2*D - a2*B*d1 + a2*B*c1 
               + a2*C*d1 - a2*D*c1 - a2*b1*C + a2*b1*D)
        return k*thesum
    
    def quadray(self):
        return self.coords
    
class Edge (object):

    """
    Edges are defined by two Vectors (above) and express as cylinder via draw().
    """

    radius = 0.03
    color = (1,0,0)

    def __init__(self, v0, v1, color=None):
        if not color==None:
            self.color = color
        self.v0 = v0
        self.v1 = v1
        self.cyl = None        

    def draw(self):
        """define and render the cylinder"""
        temp = (self.v1 - self.v0).xyz
        self.cyl = cylinder(pos = self.v0.xyz, axis = vector(*temp), 
                              radius = self.radius, color = self.color)

    def erase(self):
        """toss the cylinder"""
        if self.cyl:
            self.cyl.visible = 0
        self.cyl = None
        
    def __repr__(self):
        return 'Edge from %s to %s' % (self.v0, self.v1)


def getedges(faces):
        """
        Extract edges from the faces list.
        """ 
        edges = set()
        for f in faces:
            pairs = zip(f , f[1:]+(f[0],))
            for p in pairs:
                edges.add(tuple(sorted(p)))
        return list(edges)                

def xyplotter(domain, f):
    """
    domain should be an initialized generator, ready for next() triggering.
    f is any function of x.  Consecutive Vectors trace connected edges.
    """
    x0 = domain.next()    
    y0  = f(x0)        
    while True:
        x1 = domain.next()
        y1 =  f(x1)
        e = Edge( Vector((x0, y0, 0)), Vector((x1, y1, 0)) )
        e.draw()
        yield None
        x0, y0 = x1, y1
        
def axes(x=0,y=0,z=0):
    """
    Draw some axes on the VPython canvas
    """
    v0 = Vector((x,0,0))
    v0.draw()
    v0 = Vector((-x,0,0))
    v0.draw()

    v0 = Vector((0,y,0))
    v0.draw()
    v0 = Vector((0,-y,0))
    v0.draw()

    v0 = Vector((0,0,z))
    v0.draw()
    v0 = Vector((0,0,-z))
    v0.draw()

    
def dgen(start, step):
    """
    generic domain generator
    """
    while True:
       yield start
       start += step

       
def testme():
    """
    >>> from stickworks import testme
    Visual 2005-01-08
    >>> testme()

    See:
    http://www.4dsolutions.net/ocn/graphics/cosines.png
    """
    
    from math import cos

    def f(x):  return cos(x)

    d = dgen(-5, 0.1)
    axes(-5,1,0)
    graph = xyplotter(d, f)

    for i in xrange(100):
        graph.next()

def testmemore():
    """
    See:
    http://www.4dsolutions.net/ocn/graphics/pycalculus.png
    """

    def snakeywakey(x):
        """
        Polynomial with x-axis crossings at 3,2,-3,-7, with scaler
        to keep y-values under control (from a plotting point of view)
        """
        return 0.01 * (x-3)*(x-2)*(x+3)*(x+7)

    def deriv(f, h=1e-5):
        """
        Generic df(x)/dx approximator (discrete h)
        """
        def funk(x):
            return (f(x+h)-f(x))/h
        return funk
    

    d1 = dgen(-8, 0.1)
    d2 = dgen(-8, 0.1)
    d3 = dgen(-8, 0.1)
    
    axes(-8,5,3)

    deriv_snakeywakey = deriv(snakeywakey)
    second_deriv = deriv(deriv_snakeywakey)
    
    graph1 = xyplotter(d1, snakeywakey)
    graph2 = xyplotter(d2, deriv_snakeywakey)
    graph3 = xyplotter(d3, second_deriv)    

    Edge.color = (1,0,0)  # make snakeywakey red
    
    for i in xrange(130):
        graph1.next()

    Edge.color = (0,1,0)  # make derivative green
    
    for i in xrange(130):
        graph2.next()

    Edge.color = (0,1,1)  # make 2nd derivative cyan

    for i in xrange(130):
        graph3.next()

if __name__ == '__main__':
    # testme()
    testmemore()

testmemore( )


PREV | NEXT | TOC | HOME