by K. Urner, 4D Solutions
Feb  4, 2001: appended Fraction class
Oct  7, 2000: recursive version of mkdomain, rewrite mkfunction, mkderiv
              and canned functions to handle optional args w/ apply
May  4, 2000: added semi-circle
May  3, 2000: overhauled/generalized to accommodate parametric, catenary,
              derivative features -- more consistent use of domain concept
Mar 23, 2000: modified sinewave function to match linked Java applet
Mar 05, 2000: removed any hard-coded pov file references

import math, povray, primes
from coords import Vector

def mkdomain(low, high, interval):
    Builds domain from low to high, stepping by interval
    output = [low]     # the output list
    i=1                # number of intervals from low
    while output[-1]+interval<=high:
        output.append(low + i*interval)
        i=i+1          # increment intervals counter
    return output

def mkdomain1(start,stop,step=1):
    Recursive, uses addition -- would seem more straightforward
    than mkdomain's multiplication but accumulates floating point
    "dryer lint"
    if start>stop: return []	
    else: return [start] + mkdomain(start+step,stop,step)
def mkfunction(domain, rule, args={}):
    accept list of domain values, rule, optional arguments
    return (domain,range) pairs of corresponding function
    output = []
    for x in domain:
    return output

def mkderiv(domain,rule,args={},h=1e-10):
    Approximate derivative function over domain, allowing optional 
    constants to the function
    pairs = []
    for x in domain:  # for each member of domain...
        rvalue = (apply(rule,(x+h,),args)-apply(rule,(x-h,),args))/(2*h)
    return pairs

def mkpara(domain,rule1,rule2,rule3):
    parametric function builder
    accept list of domain values, rules 
    return (domain,range) pairs of corresponding function
    outlist = []
    for t in domain:  # for each member of domain...
        rvalue  = (rule1(t),rule2(t),rule3(t)) # ... compute range
        outlist.append((t,rvalue)) # append tuple
    return outlist

def parabola(x,a=1,d=0):
    return a*(x**2) - d

def catenary(x,k=1):
    return math.cosh(k*x)
def sinewave(x,a=1,f=1,c=0):  
    return a * math.sin(x*f - c)

def semicircle(x,r=1):
    return (r**2 - x**2) ** 0.5

def getlength(myfunc):
    approximates length of curve as sum of segments
    for list of (x,f(x)) tuples
    length = 0
    for i in range(len(myfunc)-1):
        v1 = Vector((myfunc[i]+(0,)))
        v2 = Vector((myfunc[i+1]+(0,)))
        length = length + (v2-v1).length()
    return length

def getlength1(myfunc,deltax):
    approximates length of curve as sum of
    root( 1+(df(x)/dx)^2 ) * deltax for list
    of (x,df(x)/dx) tuples (discrete)
    sum = 0
    for i in myfunc:
        sum = sum + deltax * ((1 + i[1]**2 )**0.5)
    return sum
def graphfunc(myfunc, myfile):
    draws function for list of (x,f(x)) tuples
    for i in range(len(myfunc)-1):
        v1 = Vector((myfunc[i][0],myfunc[i][1],0))
        v2 = Vector((myfunc[i+1][0],myfunc[i+1][1],0))
        myfile.edge(v1,v2)   # draw edge between the two

def graphpoints(myfunc, myfile):
    draws points for list of (x,f(x)) tuples
    for i in myfunc:
        v1 = Vector((i[0],i[1],0))
        myfile.point(v1)      # draw point 
def graphpara(myfunc, myfile):
   draws function for list of (t,(f(t),g(t),h(t))) tuples
   for i in range(len(myfunc)-1):
       v1 = Vector(myfunc[i][1])    # vector of (t,vector) pair
       v2 = Vector(myfunc[i+1][1])  # vector from next pair
       myfile.edge(v1,v2)           # draw edge between the two
def xyzaxes(file,n):

    tmp = file.cylcolor
    posx = Vector((n,0,0))
    file.cylcolor = 'Green'
    posy = Vector((0,n,0))
    file.cylcolor = 'Orange'
    posz = Vector((0,0,n))
    file.cylcolor = 'Cyan'

    # leave color unchanged 
    file.cylcolor = tmp

def example0():
    graphs sine wave from -pi to pi
    domain = mkdomain(-math.pi,math.pi,0.01)   
    function = mkfunction(domain,sinewave)
    myfile = povray.Povray("sinewave.pov",20)    
    graphfunc(function, myfile)

def example1():
    graphs sine wave points from -pi to pi
    set = mkdomain(-math.pi,math.pi,0.01)   
    function = mkfunction(set,sinewave,{'a':3,'f':4})
    myfile = povray.Povray("sine1.pov",28,0)    
    graphpoints(function, myfile)
def example2():
    graphs catenary and derivative
    domain = mkdomain(-2,2,0.01)   
    function = mkfunction(domain,catenary,{'k':2.0/3.0})
    dfunction = mkderiv(domain,catenary,{'k':2.0/3.0},0.01)
    print getlength(function)
    print getlength1(dfunction,0.01)
    myfile = povray.Povray("testgraph.pov",8,0)    
    graphfunc(function, myfile)
    myfile.cylcolor = 'Red'
    graphfunc(dfunction, myfile)

class Fraction:
    p   = 1L
    q   = 1L

    def mkfract(self,n):
        if type(n).__name__ in ['int','float','long int']:
            return Fraction(n)
        else: return n

    def __init__(self,num=1L,den=1L):
        self.p = long(num)
        self.q = long(den)

    def __mul__(self,other):
        fract = self.mkfract(other)
        return Fraction(self.p * fract.p,
                        self.q * fract.q)
    def __div__(self,n):
        # divide self by another fraction
        fract = self.mkfract(n)
        recip = Fraction(fract.q,fract.p)
        return self.__mul__(recip)
    def simplify(self):
        # reduce numerator/denominator to lowest terms
        divisor = primes.gcd(self.p,self.q)
        if abs(divisor) > 1:
            self.p /= divisor
            self.q /= divisor
    def __add__(self,n):
        # add self to another fraction
        fract  = self.mkfract(n)
        common = primes.lcm(self.q,fract.q)
        term1  =  self.p * common/self.q
        term2  = fract.p * common/fract.q
        return Fraction(term1+term2,common)

    def __sub__(self,n):
        # subtract another fraction from self
        return self.__add__(-n)

    def __neg__(self):
        # negate self
        return Fraction(-self.p,self.q)
    __rmul__ = __mul__

    def __repr__(self):
        return str(self.p)+".0/"+str(self.q)+".0"
    def list(self):
        return [self.p,self.q]
    def float(self):
        return (self.p*1.0)/(self.q*1.0)