#
#
# Import plotting routines
from numpy import *
from pylab import *
from visual import *
from visual.graph import *
from RandomArray import *
from ODEIntegrator import RK3DIntegrator
from Scientific.Geometry import Vector
from math import sin, cos, pi

def GenOmega(omegamag,theta):
   wx = omegamag*sin(theta)
   wy = 0
   wz = omegamag*cos(theta)
   return wx,wy,wz

def GenerateSphere(pts):
   r = 1.0
   xmin = -1.0;xmax = 1.0
   ymin = -1.0;ymax = 1.0
   zmin = -1.0;zmax = 1.0
   dx = (xmax-xmin)/(pts-1)
   dy = (ymax-ymin)/(pts-1)
   dz = (zmax-zmin)/(pts-1)
#
#  generate tuple for sphere skeleton
   xyplane = []
   yzplane = []
   xzplane = []
   for n in xrange(0,pts):
#     xy-plane
      x = xmin + dx*(n)
      y = sqrt(r**2 - x**2)
      z = 0.0
      xyplane.append((x,y,z))
#     yz-plane
      x = 0.0
      y = ymin + dy*(n)
      z = sqrt(r**2 - y**2)
      yzplane.append((x,y,z))
#     xz-plane
      x = xmin + dx*(n)
      y = 0.0
      z = sqrt(r**2 - x**2)
      xzplane.append((x,y,z))
#  complete sphere skeleton
   for n in xrange(0,pts):
#     xy-plane
      x = xmax - dx*(n)
      y = -sqrt(r**2 - x**2)
      z = 0.0
      xyplane.append((x,y,z))
#     yz-plane
      x = 0.0
      y = ymax - dy*(n)
      z = -sqrt(r**2 - y**2)
      yzplane.append((x,y,z))
#     xz-plane
      x = xmax - dx*(n)
      y = 0.0
      z = -sqrt(r**2 - x**2)
      xzplane.append((x,y,z))

   return xyplane,yzplane,xzplane
#
def Generate(funcX,funcY,funcZ,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter):
   print """
   Right button drag to rotate camera to view scene.
   Middle button drag up or down to zoom in or out
   """
   scene.title = ftitle 
   #scene.center = vector(0,0,25)
   scene.center = vector(0,0,0)
   # If you can cross your eyes to see stereo, uncomment the next two lines.
   #scene.stereo = 'crosseyed'
   #scene.stereodepth = 2
   scene.height = 600
   scene.width = 600
   #scene.background = (1,1,1)#white
   scene.background = (0,0,0)
   scene.forward = (0,1,-0.5)
   scene.lights = [ vector(-0.25, -0.5, -1.0) ,vector(0.25, 0.5, 1.0) ]
   
   # Draw axes
   curve( pos = [ (-2,0,0), (2,0,0) ], color = (1,1,1), radius = 0.02 )
   curve( pos = [ (0,-2,0), (0,2,0) ], color = (1,1,1), radius = 0.02 )
   curve( pos = [ (0,0,-2), (0,0,2) ], color = (1,1,1), radius = 0.02 )
   curve( pos = [ (-3,0,0), (-2.5,0.0,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-3,0,0), (-3.0,0.5,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-3,0,0), (-3.0,0.0,0.5) ], color=color.red, radius = 0.01 )
   label( pos = ( -2.6,0.0,0.0), text = ' X', box = 0)
   label( pos = ( -3.0,0.5,0.0), text = ' Y', box = 0)
   label( pos = ( -3.0,0.0,0.5), text = ' Z', box = 0)
   
   function = curve( color = color.black, radius=0.02 )
   
   y     = vector(x0,y0,z0)
   yold  = vector(x0,y0,z0)
   dydt  = vector(0,0,0)
   #theta = 3.0 * pi / 4.0
   
   #for t in arange(0,iter,dt):
   for t in xrange(0,iter):
     y = RK3DIntegrator(alpha,wx,wy,wz,y[0],y[1],y[2],
                        funcX,
                        funcY,
                        funcZ,dt)
   
     dydt[0] = (y[0]-yold[0])/dt
     dydt[1] = (y[1]-yold[1])/dt
     dydt[2] = (y[2]-yold[2])/dt
     yold = y
     # Draw lines colored by speed
     c = clip( [mag(dydt) * 1.5], 0, 1 )[0]
     #xt = cos(theta) * y[0] - sin(theta) * y[1]
     #yt = sin(theta) * y[0] + cos(theta) * y[1]
     #function.append( pos=(xt,yt,y[2]), color=(c,0.1, 1-c) )
     function.append( pos=(y[0],y[1],y[2]), color=(c,0.1, 1-c) )
#
def GenerateDotSpread(funcX,funcY,funcZ,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter):

   print ""
   print "Please enter the number if initial conditions"
   nICs   = input()
   
   print """
   Right button drag to rotate camera to view scene.
   Middle button drag up or down to zoom in or out
   """
   scene.title = ftitle 
   scene.center = vector(0,0,0)
   # If you can cross your eyes to see stereo, uncomment the next two lines.
   #scene.stereo = 'crosseyed'
   #scene.stereodepth = 2
   scene.height = 600
   scene.width = 600
   scene.background = (0,0,0)#black
   #scene.background = (1,1,1)#white
   scene.forward = (0,1,-0.5)
   scene.lights = [ vector(-0.25, -0.5, -1.0) ,vector(0.25, 0.5, 1.0) ]
   
   # Draw axes
   curve( pos = [ (-1,0,0), (1,0,0) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (0,-1,0), (0,1,0) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (0,0,-1), (0,0,1) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (-2,0,0), (-1.5,0.0,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-2,0,0), (-2.0,0.5,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-2,0,0), (-2.0,0.0,0.5) ], color=color.red, radius = 0.01 )
   label( pos = ( -1.6,0.0,0.0), text = ' X', box = 0)
   label( pos = ( -2.0,0.5,0.0), text = ' Y', box = 0)
   label( pos = ( -2.0,0.0,0.5), text = ' Z', box = 0)
#
## cannot generate transparent droplet shell
#  ball = sphere(pos=(0,0,0),radius = 1,opacity = 0.0)

#  draw a cube around the droplet
   xy,yz,xz = GenerateSphere(101)
         
   curve(pos=xy,color=color.blue,radius=0.02)
   curve(pos=yz,color=color.blue,radius=0.02)
   curve(pos=xz,color=color.blue,radius=0.02)
   
   # The number of iterations to throw away
   nTransients = 200
   # The number of time steps to integrate over
   nIterates = 1000
   
   xState,yState,zState = RK3DIntegrator(alpha,wx,wy,wz,x0,y0,z0,\
                                  funcX,funcY,funcZ,dt)
   
   # The number of time steps to integrate over
   nIterates = iter
   # Number of initial conditions in ensemble
   #nICs = 500
   # Radius of ensemble
   eRadius = 0.03
   State = [ ]
   for i in xrange(0,nICs):
      State.append(sphere(pos=(xState+normal(0.,eRadius),\
         yState+normal(0.,eRadius),zState+normal(0.,eRadius)),
         color=(1,0.2,0.2),radius=0.01))
   for n in xrange(0,nIterates):
      for i in xrange(0,nICs):
         State[i].pos[0],State[i].pos[1],State[i].pos[2] =\
            RK3DIntegrator(alpha,wx,wy,wz,State[i].pos[0],State[i].pos[1],State[i].pos[2],\
                           funcX,funcY,funcZ,dt)
#      rate(250)

def GenerateCurve(funcX,funcY,funcZ,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter,opt):
   if opt == 0:
      show_curve=False
   else:
      show_curve=True

   print """
   Right button drag to rotate camera to view scene.
   Middle button drag up or down to zoom in or out
   """
#
   if show_curve:
      graph = gdisplay(title=ftitle)
      #xzcurve1 = gcurve(color=color.red)
      xzcurve1 = gdots(color=color.white)
#

   scene.title = ftitle 
   #scene.center = vector(0,0,25)
   scene.center = vector(0,0,0)
   # If you can cross your eyes to see stereo, uncomment the next two lines.
   #scene.stereo = 'crosseyed'
   #scene.stereodepth = 2
   scene.height = 600
   scene.width = 600
   #scene.background = (1,1,1)#white
   scene.background = (0,0,0)#black
   scene.forward = (0,1,-0.5)
   scene.lights = [ vector(-0.25, -0.5, -1.0) ,vector(0.25, 0.5, 1.0) ]
   
   # Draw axes
   curve( pos = [ (-2,0,0), (2,0,0) ], color = (1,1,1), radius = 0.02 )
   curve( pos = [ (0,-2,0), (0,2,0) ], color = (1,1,1), radius = 0.02 )
   curve( pos = [ (0,0,-2), (0,0,2) ], color = (1,1,1), radius = 0.02 )
   curve( pos = [ (-3,0,0), (-2.5,0.0,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-3,0,0), (-3.0,0.5,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-3,0,0), (-3.0,0.0,0.5) ], color=color.red, radius = 0.01 )
   label( pos = ( -2.6,0.0,0.0), text = ' X', box = 0)
   label( pos = ( -3.0,0.5,0.0), text = ' Y', box = 0)
   label( pos = ( -3.0,0.0,0.5), text = ' Z', box = 0)
   
   function = curve( color = color.black, radius=0.02 )
   
   y1     = vector(x0,y0,z0)
   y1old  = vector(x0,y0,z0)
   dydt1  = vector(0,0,0)
   y2     = vector(-x0,y0,z0)
   y3     = vector(x0,y0,-z0)
   y4     = vector(-x0,y0,-z0)
   
   #for t in arange(0,iter,dt):
   hist = 0
   for t in xrange(0,iter):
     hist += 1
     y1 = RK3DIntegrator(alpha,wx,wy,wz,y1[0],y1[1],y1[2],
                         funcX,
                         funcY,
                         funcZ,dt)
     dydt1[0] = (y1[0]-y1old[0])/dt
     dydt1[1] = (y1[1]-y1old[1])/dt
     dydt1[2] = (y1[2]-y1old[2])/dt
     # Draw lines colored by speed
     c = clip( [mag(dydt1) * 1.5], 0, 1 )[0]
     function.append( pos=(y1[0],y1[1],y1[2]), color=(c,0.1, 1-c) )
#     if hist < 100:
#        function.append( pos=(y1[0],y1[1],y1[2]), color=(c,0.1, 1-c) )
#     else:
#        hist = 0
#        function = curve( color = color.black, radius=0.02 )
#        function.append( pos=(y1[0],y1[1],y1[2]), color=(c,0.1, 1-c) )
#
     if show_curve:
        test = y1[1]*y1old[1]
        if test <= 0.0:
           xzcurve1.plot(pos=(y1[0],y1[2]))
#
     y1old = y1
#
def GenerateAutoPoincareMap(funcX,funcY,funcZ,ftitle,dt,iter):
   print """
   Generating Poincare Maps for the xz-Plane
   """
   nplot = 2
   a     = 1.0
   x0    = 0.3
   y0    = 0.5
   z0    = 0.3
#   Param = [(a,0.2938900,wy,0.404508000,x0,y0,z0,dt,iter),
#            (a,0.5877850,wy,0.809016990,x0,y0,z0,dt,iter),
#            (a,0.6641973,wy,0.914189203,x0,y0,z0,dt,iter),
#            (a,0.6877100,wy,0.946549800,x0,y0,z0,dt,iter),
#            (a,0.7347000,wy,1.011270000,x0,y0,z0,dt,iter),
#            (a,0.8228990,wy,1.132600000,x0,y0,z0,dt,iter),
#            (a,1.0286000,wy,1.415770000,x0,y0,z0,dt,iter),
#            (a,1.4694600,wy,2.022540000,x0,y0,z0,dt,iter),
#            (a,2.3411400,wy,3.236067900,x0,y0,z0,dt,iter),
#            (a,4.7022820,wy,6.472135900,x0,y0,z0,dt,iter),
#            (a,0.0003142,wy,0.099999990,x0,y0,z0,dt,iter),
#            (a,0.0031410,wy,0.099999990,x0,y0,z0,dt,iter),
#            (a,0.0309000,wy,0.095150000,x0,y0,z0,dt,iter),
#            (a,0.0707106,wy,0.070710600,x0,y0,z0,dt,iter),
#            (a,0.0951056,wy,0.030901600,x0,y0,z0,dt,iter),
#            (a,0.0062831,wy,1.999999900,x0,y0,z0,dt,iter),
#            (a,0.0125662,wy,1.999961000,x0,y0,z0,dt,iter),
#            (a,0.0420942,wy,1.999556900,x0,y0,z0,dt,iter),
#            (a,0.0628215,wy,1.999013100,x0,y0,z0,dt,iter),
#            (a,0.1255810,wy,1.996053400,x0,y0,z0,dt,iter),
#            (a,1.1755700,wy,1.618033988,x0,y0,z0,dt,iter),
#            (a,1.9021130,wy,0.618033988,x0,y0,z0,dt,iter),
#            (a,1.9960530,wy,0.125581030,x0,y0,z0,dt,iter),
#            (a,2.0000000,wy,3.000000000,x0,y0,z0,dt,iter),
#            (a,2.0000000,wy,3.000000000,0.0,y0,z0,dt,iter)]
   Param = [(a,0.10,0.2*pi,x0,y0,z0,dt,iter),
            (a,0.50,0.2*pi,x0,y0,z0,dt,iter),
            (a,1.00,0.2*pi,x0,y0,z0,dt,iter),
            (a,1.13,0.2*pi,x0,y0,z0,dt,iter),
            (a,1.17,0.2*pi,x0,y0,z0,dt,iter),
            (a,1.25,0.2*pi,x0,y0,z0,dt,iter),
            (a,1.40,0.2*pi,x0,y0,z0,dt,iter),
            (a,1.75,0.2*pi,x0,y0,z0,dt,iter),
            (a,2.50,0.2*pi,x0,y0,z0,dt,iter),
            (a,4.00,0.2*pi,x0,y0,z0,dt,iter),
            (a,8.00,0.2*pi,x0,y0,z0,dt,iter),
            (a,0.10,0.001*pi,x0,y0,z0,dt,iter),
            (a,0.10,0.100*pi,x0,y0,z0,dt,iter),
            (a,0.10,0.170*pi,x0,y0,z0,dt,iter),
            (a,0.10,0.250*pi,x0,y0,z0,dt,iter),
            (a,0.10,0.270*pi,x0,y0,z0,dt,iter),
            (a,0.10,0.290*pi,x0,y0,z0,dt,iter),
            (a,0.10,0.333*pi,x0,y0,z0,dt,iter),
            (a,0.10,0.400*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.001*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.002*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.0067*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.010*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.020*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.060*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.100*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.200*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.300*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.350*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.400*pi,x0,y0,z0,dt,iter),
            (a,2.00,0.480*pi,x0,y0,z0,dt,iter)]
   for n in Param:
      print "alpha, omega x, omega y, omega z, x, y, z, dt, timesteps"
      wx,wy,wz = GenOmega(n[1],n[2])
      print n[0],wx  ,wy  ,wz  ,n[3],n[4],n[5],n[6],n[7]
      PoincareMap(funcX,funcY,funcZ,ftitle,n[0],wx  ,wy  ,wz  \
                                               ,n[3],n[4],n[5]\
                                               ,n[6],n[7],nplot   )
      print "... Done"

# define the function
def PoincareMap(funcX,funcY,funcZ,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter,nplot):
   
   x    = Vector(x0,y0,z0)
   xold = Vector(x0,y0,z0)
   xz1plot = [ ]   # The iterates
   xz2plot = [ ]
   xy1plot = [ ]   # The iterates
   xy2plot = [ ]
   yz1plot = [ ]   # The iterates
   yz2plot = [ ]
   
   for t in xrange(0,iter):
     xold = x
     x = RK3DIntegrator(alpha,wx,wy,wz,x[0],x[1],x[2],
                         funcX,
                         funcY,
                         funcZ,dt)
#    xz-plane
     test = x[1]*xold[1]
     if test <= 0.0:
        xz1plot.append(x[0])
        xz2plot.append(x[2])

##    xy-plane
#     test = x[2]*xold[2]
#     if test <= 0.0:
#        xy1plot.append(x[0])
#        xy2plot.append(x[1])
#
##    yz-plane
#     test = x[0]*xold[0]
#     if test <= 0.0:
#        yz1plot.append(x[1])
#        yz2plot.append(x[2])

  # xz-plane
  # make background black
#   subplot(311,axisbg=(0,0,0))
   #subplot(111,axisbg=(0,0,0))
   subplot(111,axisbg=(1,1,1))
  # call plot function
   #plot(xz1plot,xz2plot, 'w.')
   plot(xz1plot,xz2plot, 'k.')
#
#  # xy-plane
#  # make background black
#   subplot(312,axisbg=(0,0,0))
#  # call plot function
#   plot(xy1plot,xy2plot, 'w.')
#
#  # yz-plane
#  # make background black
#   subplot(313,axisbg=(0,0,0))
#  # call plot function
#   plot(yz1plot,yz2plot, 'w.')

  # set figure xz limits
   xlim(xmin=-1.0)
   xlim(xmax= 1.0)
   ylim(ymin=-1.0)
   ylim(ymax= 1.0)

  # Setup the plot
   figure(1,(8,6.8))
   TitleString = ftitle + ':' + '(wx,wy,wz) ' + str(wx) + ',' + str(wy) + ',' + str(wz)
   title(TitleString)
   xlabel('x' )   # set x-axis label
   ylabel('z')    # set z-axis label

   if nplot == 1:
      # Display plot in window
      show()
   elif nplot == 2:
      # save picture to file
      filename = 'ChaoticDrop_'+str(alpha)+'_'+str(wx)+'_'+str(wy)+'_'+str(wz)+'.ps'
      savefig(filename, format='ps', dpi = 800)

   # clear figure
   cla()
 
def GenDotSpreadPCurve(funcX,funcY,funcZ,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter):
   show_curve=True

   print ""
   print "Please enter the number if initial conditions"
   nICs   = input()
   
   print """
   Right button drag to rotate camera to view scene.
   Middle button drag up or down to zoom in or out
   """
#
   if show_curve:
      graph = gdisplay(title=ftitle)
      xzcurve1 = gdots(color=color.white)
#
   scene.title = ftitle 
   scene.center = vector(0,0,0)
   # If you can cross your eyes to see stereo, uncomment the next two lines.
   #scene.stereo = 'crosseyed'
   #scene.stereodepth = 2
   scene.height = 600
   scene.width = 600
   scene.background = (0,0,0)#black
   #scene.background = (1,1,1)#white
   scene.forward = (0,1,-0.5)
   scene.lights = [ vector(-0.25, -0.5, -1.0) ,vector(0.25, 0.5, 1.0) ]
   
   # Draw axes
   curve( pos = [ (-1,0,0), (1,0,0) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (0,-1,0), (0,1,0) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (0,0,-1), (0,0,1) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (-2,0,0), (-1.5,0.0,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-2,0,0), (-2.0,0.5,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-2,0,0), (-2.0,0.0,0.5) ], color=color.red, radius = 0.01 )
   label( pos = ( -1.6,0.0,0.0), text = ' X', box = 0)
   label( pos = ( -2.0,0.5,0.0), text = ' Y', box = 0)
   label( pos = ( -2.0,0.0,0.5), text = ' Z', box = 0)
#
## cannot generate transparent droplet shell
#  ball = sphere(pos=(0,0,0),radius = 1,opacity = 0.0)

#  draw a cube around the droplet
   xy,yz,xz = GenerateSphere(101)
         
   curve(pos=xy,color=color.blue,radius=0.02)
   curve(pos=yz,color=color.blue,radius=0.02)
   curve(pos=xz,color=color.blue,radius=0.02)
   
   # The number of iterations to throw away
   nTransients = 200
   # The number of time steps to integrate over
   nIterates = 1000
   
   for n in xrange(0,nTransients):
      xState,yState,zState = RK3DIntegrator(alpha,wx,wy,wz,x0,y0,z0,\
                                     funcX,funcY,funcZ,dt)
   
   # The number of time steps to integrate over
   nIterates = iter
   # Number of initial conditions in ensemble
   #nICs = 500
   # Radius of ensemble
   eRadius = 0.03
   State = [ ]
   y     = [x0,y0,z0]
   yold  = [x0,y0,z0]
   for i in xrange(0,nICs):
      State.append(sphere(pos=(xState+normal(0.,eRadius),\
         yState+normal(0.,eRadius),zState+normal(0.,eRadius)),
         color=(1,0.2,0.2),radius=0.01))
   for n in xrange(0,nIterates):
      for i in xrange(0,nICs):
         y = RK3DIntegrator(alpha,wx,wy,wz,State[i].pos[0],State[i].pos[1],State[i].pos[2],\
                            funcX,funcY,funcZ,dt)
         State[i].pos[0] = y[0]
         State[i].pos[1] = y[1]
         State[i].pos[2] = y[2]
#
         if show_curve:
            test = y[1]*yold[1]
            #if test <= 0.0 and i == 0:
            if test <= 0.0:
               xzcurve1.plot(pos=(y[0],y[2]))
         yold = y
#      rate(250)
 
def GenPCurve(funcX,funcY,funcZ,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter):
   show_curve=True

   #print ""
   #print "Please enter the number if initial conditions"
   #nICs   = input()
   nICs = 1

   print """
   Middle button drag up or down to zoom in or out
   """
#
   if show_curve:
      graph = gdisplay(title=ftitle)
      xzcurve1 = gdots(color=color.white)
#
   # The number of iterations to throw away
   nTransients = 200
   # The number of time steps to integrate over
   nIterates = 1000
   
   for n in xrange(0,nTransients):
      xState,yState,zState = RK3DIntegrator(alpha,wx,wy,wz,x0,y0,z0,\
                                     funcX,funcY,funcZ,dt)
   
   # The number of time steps to integrate over
   nIterates = iter
   # Number of initial conditions in ensemble
   #nICs = 500
   # Radius of ensemble
   snICs = nICs/4
   eRadius = 0.1
   dR = eRadius/nICs
   State = [ ]
   y     = [x0,y0,z0]
   yold  = [x0,y0,z0]
#   for i in xrange(0,snICs):
#      State.append(( xState+dR*i, yState+dR*i, zState+dR*i))
#   for i in xrange(0,snICs):
#      State.append((-xState+dR*i, yState+dR*i, zState+dR*i))
#   for i in xrange(0,snICs):
#      State.append(( xState+dR*i,-yState+dR*i, zState+dR*i))
#   for i in xrange(0,snICs):
#      State.append(( xState+dR*i, yState+dR*i,-zState+dR*i))
   for i in xrange(0,nICs):
      State.append(( xState+dR*i, yState+dR*i, zState+dR*i))

   for n in xrange(0,nIterates):
      for i in xrange(nICs):
         y = State[i]
         y = RK3DIntegrator(alpha,wx,wy,wz,y[0],y[1],y[2],\
                            funcX,funcY,funcZ,dt)
         State[i] = y
#
         if show_curve:
            test = y[1]*yold[1]
            if test <= 0.0:
               xzcurve1.plot(pos=(y[0],y[2]))
         yold = y
#      rate(250)

def DxDtDot(alpha,wx,wy,wz,x,y,z,dx,dy,dz):
        dfdx = 0.5*(9.0*x*x + (5.0-2.0*alpha)*y*y\
                 + (7.0+2.0*alpha)*z*z - 3.0)/(1+alpha)
        dfdy = (5.0-2.0*alpha)*x*y/(1+alpha) - 0.5*wz
        dfdz = (7.0+2.0*alpha)*x*z/(1+alpha) + 0.5*wy
        Df = dfdx*dx + dfdy*dy + dfdz*dz
	return Df

def DyDtDot(alpha,wx,wy,wz,x,y,z,dx,dy,dz):
        dgdx = 0.5*(10.0*x*y*alpha - 4.0*x*y)/(1+alpha)\
                 + 0.5*wz
        dgdy = 0.5*(9.0*y*y + (5.0*alpha-2.0)*x*x\
                 + (7.0*alpha+2.0)*z*z - 3.0*alpha)/(1+alpha)
        dgdz = 0.5*(14.0*alpha+4.0)*y*z/(1+alpha) - 0.5*wx
        Dg = dgdx*dx + dgdy*dy + dgdz*dz
	return Dg

def DzDtDot(alpha,wx,wy,wz,x,y,z,dx,dy,dz):
        dhdx = 0.5*(-10.0*x*z - 4.0*x*z/(1+alpha)-wy)
        dhdy = 0.5*(-10.0*y*z - 4.0*alpha*y*z/(1+alpha)+wx)
        dhdz = 0.5*(-9.0*z*z + (-7.0-5.0*alpha)*x*x/(1+alpha)\
                    +(-5.0-7.0*alpha)*y*y/(1+alpha) + 3.0)
        Dh = dhdx*dx + dhdy*dy + dhdz*dz
	return Dh

# Volume contraction given by
#	 Trace(Jacobian(x,y,z))
def DropODETrJac(alpha,wx,wy,wz,x,y,z):
        dfdx = 0.5*(9.0*x*x + (5.0-2.0*alpha)*y*y\
                 + (7.0+2.0*alpha)*z*z - 3.0)/(1+alpha)
        dgdy = 0.5*(9.0*y*y + (5.0*alpha-2.0)*x*x\
                 + (7.0*alpha+2.0)*z*z - 3.0*alpha)/(1+alpha)
        dhdz = 0.5*(-9.0*z*z + (-7.0-5.0*alpha)*x*x/(1+alpha)\
                    +(-5.0-7.0*alpha)*y*y/(1+alpha) + 3.0)
	return dfdx + dgdy + dhdz

# The fractal dimension from the LCEs (Kaplan-Yorke conjecture)
#   Assume these are ordered: LCE1 >= LCE2 >= LCE3
def FractalDimension3DODE(LCE1,LCE2,LCE3):
	# "Close" to zero ... we're estimating here
	Err = 0.01
	if LCE1 < -Err:	  # Stable fixed point    (-,-,-)
		return 0.0
	elif abs(LCE1) <= Err:
		if LCE2 < -Err:  # Limit cycle	      (0,-,-)
			return 1.0
		else:	        # Torus               (0,0,-)
			return 2.0
	else:	            # Chaotic attractor   (+,0,-)
		return 2.0 + (LCE1+LCE2) / abs(LCE3)

# Tanget space flow (using fourth-order Runge-Kutta integrator)
def TangentFlowRKThreeD(a,wx,wy,wz,x,y,z,df,dg,dh,dx,dy,dz,dt):
	k1x = dt * df(a,wx,wy,wz,x,y,z,dx,dy,dz)
	k1y = dt * dg(a,wx,wy,wz,x,y,z,dx,dy,dz)
	k1z = dt * dh(a,wx,wy,wz,x,y,z,dx,dy,dz)
	k2x = dt * df(a,wx,wy,wz,x,y,z,dx+k1x/2.0,dy+k1y/2.0,dz+k1z/2.0)
	k2y = dt * dg(a,wx,wy,wz,x,y,z,dx+k1x/2.0,dy+k1y/2.0,dz+k1z/2.0)
	k2z = dt * dh(a,wx,wy,wz,x,y,z,dx+k1x/2.0,dy+k1y/2.0,dz+k1z/2.0)
	k3x = dt * df(a,wx,wy,wz,x,y,z,dx+k2x/2.0,dy+k2y/2.0,dz+k2z/2.0)
	k3y = dt * dg(a,wx,wy,wz,x,y,z,dx+k2x/2.0,dy+k2y/2.0,dz+k2z/2.0)
	k3z = dt * dh(a,wx,wy,wz,x,y,z,dx+k2x/2.0,dy+k2y/2.0,dz+k2z/2.0)
	k4x = dt * df(a,wx,wy,wz,x,y,z,dx+k3x,dy+k3y,dz+k3z)
	k4y = dt * dg(a,wx,wy,wz,x,y,z,dx+k3x,dy+k3y,dz+k3z)
	k4z = dt * dh(a,wx,wy,wz,x,y,z,dx+k3x,dy+k3y,dz+k3z)
	dx += ( k1x + 2.0 * k2x + 2.0 * k3x + k4x ) / 6.0
	dy += ( k1y + 2.0 * k2y + 2.0 * k3y + k4y ) / 6.0
	dz += ( k1z + 2.0 * k2z + 2.0 * k3z + k4z ) / 6.0
	return dx,dy,dz
#
def GenerateAutoLCE(funcX,funcY,funcZ,ftitle,dt,iter):
   print """
   Generating LCE for all interesting parameters
   """
   nplot = 2
   a     = 1.0
   x0    = 0.1
   y0    = 0.1
   z0    = 0.1
   wy    = 0.0
   Param = [(a,0.0,wy,0.0,x0,y0,z0,dt,iter),
            (a,0.2938900,wy,0.404508000,x0,y0,z0,dt,iter),
            (a,0.5877850,wy,0.809016990,x0,y0,z0,dt,iter),
            (a,0.6641973,wy,0.914189203,x0,y0,z0,dt,iter),
            (a,0.6877100,wy,0.946549800,x0,y0,z0,dt,iter),
            (a,0.7347000,wy,1.011270000,x0,y0,z0,dt,iter),
            (a,0.8228990,wy,1.132600000,x0,y0,z0,dt,iter),
            (a,1.0286000,wy,1.415770000,x0,y0,z0,dt,iter),
            (a,1.4694600,wy,2.022540000,x0,y0,z0,dt,iter),
            (a,2.3411400,wy,3.236067900,x0,y0,z0,dt,iter),
            (a,4.7022820,wy,6.472135900,x0,y0,z0,dt,iter),
            (a,0.0003142,wy,0.099999990,x0,y0,z0,dt,iter),
            (a,0.0031410,wy,0.099999990,x0,y0,z0,dt,iter),
            (a,0.0309000,wy,0.095150000,x0,y0,z0,dt,iter),
            (a,0.0707106,wy,0.070710600,x0,y0,z0,dt,iter),
            (a,0.0951056,wy,0.030901600,x0,y0,z0,dt,iter),
            (a,0.0062831,wy,1.999999900,x0,y0,z0,dt,iter),
            (a,0.0125662,wy,1.999961000,x0,y0,z0,dt,iter),
            (a,0.0420942,wy,1.999556900,x0,y0,z0,dt,iter),
            (a,0.0628215,wy,1.999013100,x0,y0,z0,dt,iter),
            (a,0.1255810,wy,1.996053400,x0,y0,z0,dt,iter),
            (a,1.1755700,wy,1.618033988,x0,y0,z0,dt,iter),
            (a,1.9021130,wy,0.618033988,x0,y0,z0,dt,iter),
            (a,1.9960530,wy,0.125581030,x0,y0,z0,dt,iter),
            (a,2.0000000,wy,3.000000000,x0,y0,z0,dt,iter),
            (a,2.0000000,wy,3.000000000,0.0,y0,z0,dt,iter)]
   for n in Param:
      print "alpha, omega x, omega y, omega z, x, y, z, dt, timesteps"
      print n[0],n[1],n[2],n[3],n[4],n[5],n[6],n[7],n[8]
      GenerateLCE(funcX,funcY,funcZ,ftitle,n[0],n[1],n[2],n[3]\
                                               ,n[4],n[5],n[6]\
                                               ,n[7],n[8],nplot   )
      print "... Done"


def GenerateLCE(DxDt,DyDt,DzDt,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter,nplot):
   # Simulation parameters
   # Integration time step
   # The number of iterations to throw away
   nTransients = 10
   # The number of time steps to integrate over
   nIterates = iter
   
   # The main loop that generates the orbit, storing the states
   xState = x0
   yState = y0
   zState = z0
   # Iterate for some number of transients, but don't use these states
   for n in xrange(0,nTransients):
   	xState,yState,zState = RK3DIntegrator(alpha,wx,wy,wz,xState,yState,zState,DxDt,DyDt,DzDt,dt)
   # Set up array of iterates and store the current state
   x = [xState]
   y = [yState]
   z = [zState]
   for n in xrange(0,nIterates):
   	# at each time step calculate new (x,y,z)(t)
   	xt,yt,zt = RK3DIntegrator(alpha,wx,wy,wz,x[n],y[n],z[n],DxDt,DyDt,DzDt,dt)
   	# and append to lists
   	x.append(xt)
   	y.append(yt)
   	z.append(zt)
   
   # Estimate the LCEs
   # The number of iterations to throw away
   nTransients = 100
   # The number of iterations to over which to estimate
   #  This is really the number of pull-backs
   nIterates = iter
   # The number of iterations per pull-back
   nItsPerPB = 10
   # Initial condition
   xState = x0
   yState = y0
   zState = z0
   # Initial tangent vectors
   e1x = 1.0
   e1y = 0.0
   e1z = 0.0
   e2x = 0.0
   e2y = 1.0
   e2z = 0.0
   e3x = 0.0
   e3y = 0.0
   e3z = 1.0
   # Iterate away transients and let the tangent vectors align
   #	with the global stable and unstable manifolds
   for n in xrange(0,nTransients):
   	for i in xrange(nItsPerPB):
   		xState,yState,zState = RK3DIntegrator(alpha,wx,wy,wz,xState,yState,zState,\
   			DxDt,DyDt,DzDt,dt)
   		# Evolve tangent vector for maximum LCE (LCE1)
   		e1x,e1y,e1z = TangentFlowRKThreeD(alpha,wx,wy,wz,xState,yState,zState,\
   			DxDtDot,DyDtDot,DzDtDot,e1x,e1y,e1z,dt)
   		# Evolve tangent vector for next LCE (LCE2)
   		e2x,e2y,e2z = TangentFlowRKThreeD(alpha,wx,wy,wz,xState,yState,zState,\
   			DxDtDot,DyDtDot,DzDtDot,e2x,e2y,e2z,dt)
   		# Evolve tangent vector for last LCE
   		e3x,e3y,e3z = TangentFlowRKThreeD(alpha,wx,wy,wz,xState,yState,zState,\
   			DxDtDot,DyDtDot,DzDtDot,e3x,e3y,e3z,dt)
   	# Normalize the tangent vector
   	d = sqrt(e1x*e1x + e1y*e1y + e1z*e1z)
   	e1x /= d
   	e1y /= d
   	e1z /= d
   	# Pull-back: Remove any e1 component from e2
   	dote1e2 = e1x * e2x + e1y * e2y + e1z * e2z
   	e2x -= dote1e2 * e1x
   	e2y -= dote1e2 * e1y
   	e2z -= dote1e2 * e1z
   	# Normalize second tangent vector
   	d = sqrt(e2x*e2x + e2y*e2y + e2z*e2z)
   	e2x /= d
   	e2y /= d
   	e2z /= d
   	# Pull-back: Remove any e1 and e2 components from e3
   	dote1e3 = e1x * e3x + e1y * e3y + e1z * e3z
   	dote2e3 = e2x * e3x + e2y * e3y + e2z * e3z
   	e3x -= dote1e3 * e1x + dote2e3 * e2x
   	e3y -= dote1e3 * e1y + dote2e3 * e2y
   	e3z -= dote1e3 * e1z + dote2e3 * e2z
   	# Normalize third tangent vector
   	d = sqrt(e3x*e3x + e3y*e3y + e3z*e3z)
   	e3x /= d
   	e3y /= d
   	e3z /= d
   
   # Okay, now we're ready to begin the estimation
   LCE1 = 0.0
   LCE2 = 0.0
   LCE3 = 0.0
   LCE1p = [LCE1]
   LCE2p = [LCE2]
   LCE3p = [LCE3]
   iterations = [0]
   for n in xrange(0,nIterates):
   	for i in xrange(nItsPerPB):
   		xState,yState,zState = RK3DIntegrator(alpha,wx,wy,wz,xState,yState,zState,\
   			DxDt,DyDt,DzDt,dt)
   		# Evolve tangent vector for maximum LCE (LCE1)
   		e1x,e1y,e1z = TangentFlowRKThreeD(alpha,wx,wy,wz,xState,yState,zState,\
   			DxDtDot,DyDtDot,DzDtDot,e1x,e1y,e1z,dt)
   		# Evolve tangent vector for next LCE (LCE2)
   		e2x,e2y,e2z = TangentFlowRKThreeD(alpha,wx,wy,wz,xState,yState,zState,\
   			DxDtDot,DyDtDot,DzDtDot,e2x,e2y,e2z,dt)
   		# Evolve tangent vector for last LCE
   		e3x,e3y,e3z = TangentFlowRKThreeD(alpha,wx,wy,wz,xState,yState,zState,\
   			DxDtDot,DyDtDot,DzDtDot,e3x,e3y,e3z,dt)
   	# Normalize the tangent vector
   	d = sqrt(e1x*e1x + e1y*e1y + e1z*e1z)
   	e1x /= d
   	e1y /= d
   	e1z /= d
   	# Accumulate the first tangent vector's length change factor
   	LCE1 += log(d)
   	# Pull-back: Remove any e1 component from e2
   	dote1e2 = e1x * e2x + e1y * e2y + e1z * e2z
   	e2x -= dote1e2 * e1x
   	e2y -= dote1e2 * e1y
   	e2z -= dote1e2 * e1z
   	# Normalize second tangent vector
   	d = sqrt(e2x*e2x + e2y*e2y + e2z*e2z)
   	e2x /= d
   	e2y /= d
   	e2z /= d
   	# Accumulate the second tangent vector's length change factor
   	LCE2 += log(d)
   	# Pull-back: Remove any e1 and e2 components from e3
   	dote1e3 = e1x * e3x + e1y * e3y + e1z * e3z
   	dote2e3 = e2x * e3x + e2y * e3y + e2z * e3z
   	e3x -= dote1e3 * e1x + dote2e3 * e2x
   	e3y -= dote1e3 * e1y + dote2e3 * e2y
   	e3z -= dote1e3 * e1z + dote2e3 * e2z
   	# Normalize third tangent vector
   	d = sqrt(e3x*e3x + e3y*e3y + e3z*e3z)
   	e3x /= d
   	e3y /= d
   	e3z /= d
   	# Accumulate the third tangent vector's length change factor
   	LCE3 += log(d)
        # Convert to per-iterate, per-second LCEs and to base-2 logs
        IntegrationTime = dt * float(nItsPerPB) * float(n)
   	# and append to lists
        if IntegrationTime == 0.0:
           IntegrationTime = 0.01
        LCE1p.append(LCE1 / IntegrationTime)
        LCE2p.append(LCE2 / IntegrationTime)
        LCE3p.append(LCE3 / IntegrationTime)
        iterations.append(n)
   
   # Convert to per-iterate, per-second LCEs and to base-2 logs
   IntegrationTime = dt * float(nItsPerPB) * float(nIterates)
   LCE1 = LCE1 / IntegrationTime
   LCE2 = LCE2 / IntegrationTime
   LCE3 = LCE3 / IntegrationTime
   # Setup the parametric plot:
   xlabel('iterations') # set x-axis label
   ylabel('LCE') # set y-axis label
   # Construct plot title
   title('LCE iteration for ' + ftitle)
   #axis('equal')
   axis([0,1000.0,-15.0,1.1])
   # Plot the trajectory in the phase plane
   plot(iterations,LCE1p,'b')
   plot(iterations,LCE2p,'r')
   plot(iterations,LCE3p,'g')
   axhline(0.0,color = 'k')
   axvline(0.0,color = 'k')
   
   # Use command below to save figure
   #savefig('LorenzODELCE', dpi=600)
   
   # Display the plot in a window
   if nplot == 1:
      # Display plot in window
      show()
   elif nplot == 2:
      # save picture to file
      filename = 'ChaoticDropLCEHistory_'+str(alpha)+'_'+str(wx)+'_'+str(wy)+'_'+str(wz)+'.ps'
      savefig(filename, format='ps', dpi = 800)

   # clear figure
   cla()
   # Calculate contraction factor, for comparison.
   Contraction = DropODETrJac(alpha,wx,wy,wz,0,0,0)
   
   # Choose a pair of coordinates from (x,y,z) to show
   # Setup the parametric plot:
   xlabel('x(t)') # set x-axis label
   ylabel('y(t)') # set y-axis label
   # Construct plot title
   LCEString = '(%g,%g,%g)' % (LCE1,LCE2,LCE3)
   PString = '(wx,wy,wz) = (%g,%g,%g)' % (wx,wy,wz)
   CString = 'Contraction = %g, Diff = %g' % (Contraction,abs(LCE1+LCE2+LCE3-Contraction))
   FString   = 'Fractal dimension = %g' % FractalDimension3DODE(LCE1,LCE2,LCE3)
   title('4th order Runge-Kutta Method: '+ ftitle + PString + ':\n LCEs = ' + LCEString + ', ' + CString + '\n ' + FString)
   axis('equal')
   axis([-20.0,20.0,-20.0,20.0])
   # Plot the trajectory in the phase plane
   plot(x,y,'b')
   axhline(0.0,color = 'k')
   axvline(0.0,color = 'k')
   
   # Display the plot in a window
   if nplot == 2:
      # save picture to file
      filename = 'ChaoticDropLCE_'+str(alpha)+'_'+str(wx)+'_'+str(wy)+'_'+str(wz)+'.ps'
      savefig(filename, format='ps', dpi = 800)

   if nplot == 1:
      # Display plot in window
      show()
      return
   elif nplot == 2:
      # save picture to file
      filename = 'ChaoticDropLCE_'+str(alpha)+'_'+str(wx)+'_'+str(wy)+'_'+str(wz)+'.ps'
      savefig(filename, format='ps', dpi = 800)

   # clear figure
   cla()
#
def GenerateDotHistory(funcX,funcY,funcZ,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter):

   print ""
   print "Please enter the number of points for history"
   nICs   = input()
   
   print """
   Right button drag to rotate camera to view scene.
   Middle button drag up or down to zoom in or out
   """
   scene.title = ftitle 
   scene.center = vector(0,0,0)
   # If you can cross your eyes to see stereo, uncomment the next two lines.
   #scene.stereo = 'crosseyed'
   #scene.stereodepth = 2
   scene.height = 600
   scene.width = 600
   scene.background = (0,0,0)#black
   #scene.background = (1,1,1)#white
   scene.forward = (0,1,-0.5)
   scene.lights = [ vector(-0.25, -0.5, -1.0) ,vector(0.25, 0.5, 1.0) ]
   
   # Draw axes
   curve( pos = [ (-1,0,0), (1,0,0) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (0,-1,0), (0,1,0) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (0,0,-1), (0,0,1) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (-2,0,0), (-1.5,0.0,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-2,0,0), (-2.0,0.5,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-2,0,0), (-2.0,0.0,0.5) ], color=color.red, radius = 0.01 )
   label( pos = ( -1.6,0.0,0.0), text = ' X', box = 0)
   label( pos = ( -2.0,0.5,0.0), text = ' Y', box = 0)
   label( pos = ( -2.0,0.0,0.5), text = ' Z', box = 0)
#
## cannot generate transparent droplet shell
#  ball = sphere(pos=(0,0,0),radius = 1,opacity = 0.0)

#  draw a cube around the droplet
   xy,yz,xz = GenerateSphere(101)
         
   curve(pos=xy,color=color.blue,radius=0.02)
   curve(pos=yz,color=color.blue,radius=0.02)
   curve(pos=xz,color=color.blue,radius=0.02)
   
   # The number of iterations to throw away
   nTransients = 200
   # The number of time steps to integrate over
   nIterates = 1000
   
   xState,yState,zState = RK3DIntegrator(alpha,wx,wy,wz,x0,y0,z0,\
                                  funcX,funcY,funcZ,dt)
   
   # The number of time steps to integrate over
   nIterates = iter
   # Radius of ensemble
   eRadius = 0.03
   State = [ ]
   for i in xrange(0,nICs):
#      State.append(sphere(pos=(xState+normal(0.,eRadius),\
#         yState+normal(0.,eRadius),zState+normal(0.,eRadius)),
#         color=(1,0.2,0.2),radius=0.01))
      State.append(sphere(pos=(xState,yState,zState),color=(1,0.2,0.2),radius=0.01))
   for n in xrange(0,nIterates):
      for i in xrange(1,nICs):
         State[i].pos[0],State[i].pos[1],State[i].pos[2] =\
            RK3DIntegrator(alpha,wx,wy,wz,State[i-1].pos[0],State[i-1].pos[1],State[i-1].pos[2],\
                           funcX,funcY,funcZ,dt)
         if i == nICs-1:
            State[0].pos[0],State[0].pos[1],State[0].pos[2] =\
               RK3DIntegrator(alpha,wx,wy,wz,State[i].pos[0],State[i].pos[1],State[i].pos[2],\
                              funcX,funcY,funcZ,dt)
         
#         State[i].pos[0],State[i].pos[1],State[i].pos[2] =\
#            RK3DIntegrator(alpha,wx,wy,wz,State[i].pos[0],State[i].pos[1],State[i].pos[2],\
#                           funcX,funcY,funcZ,dt)
         rate(250)
#
def GenerateDotHistory2(funcX,funcY,funcZ,ftitle,alpha,wx,wy,wz,x0,y0,z0,dt,iter):

   print ""
   print "Please enter the number of points for history"
   nICs   = input()
   
   print """
   Right button drag to rotate camera to view scene.
   Middle button drag up or down to zoom in or out
   """
   scene.title = ftitle 
   scene.center = vector(0,0,0)
   # If you can cross your eyes to see stereo, uncomment the next two lines.
   #scene.stereo = 'crosseyed'
   #scene.stereodepth = 2
   scene.height = 600
   scene.width = 600
   scene.background = (0,0,0)#black
   #scene.background = (1,1,1)#white
   scene.forward = (0,1,-0.5)
   scene.lights = [ vector(-0.25, -0.5, -1.0) ,vector(0.25, 0.5, 1.0) ]
   
   # Draw axes
   curve( pos = [ (-1,0,0), (1,0,0) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (0,-1,0), (0,1,0) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (0,0,-1), (0,0,1) ], color = (0,0,1), radius = 0.02 )
   curve( pos = [ (-2,0,0), (-1.5,0.0,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-2,0,0), (-2.0,0.5,0.0) ], color=color.red, radius = 0.01 )
   curve( pos = [ (-2,0,0), (-2.0,0.0,0.5) ], color=color.red, radius = 0.01 )
   label( pos = ( -1.6,0.0,0.0), text = ' X', box = 0)
   label( pos = ( -2.0,0.5,0.0), text = ' Y', box = 0)
   label( pos = ( -2.0,0.0,0.5), text = ' Z', box = 0)
#
## cannot generate transparent droplet shell
#  ball = sphere(pos=(0,0,0),radius = 1,opacity = 0.0)

#  draw a cube around the droplet
   xy,yz,xz = GenerateSphere(101)
         
   curve(pos=xy,color=color.blue,radius=0.02)
   curve(pos=yz,color=color.blue,radius=0.02)
   curve(pos=xz,color=color.blue,radius=0.02)
   
   # The number of iterations to throw away
   nTransients = 200
   # The number of time steps to integrate over
   nIterates = 1000
   
   xState1,yState1,zState1 = RK3DIntegrator(alpha,wx,wy,wz,x0,y0,z0,\
                                       funcX,funcY,funcZ,dt)
   xState2,yState2,zState2 = RK3DIntegrator(alpha,wx,wy,wz,x0*1.001,y0*1.001,z0*1.001,\
                                       funcX,funcY,funcZ,dt)
   
   # The number of time steps to integrate over
   nIterates = iter
   # Radius of ensemble
   eRadius = 0.03
   State1 = [ ]
   State2 = [ ]
   for i in xrange(0,nICs):
      State1.append(sphere(pos=(xState1,yState1,zState1),color=color.red  ,radius=0.01))
      State2.append(sphere(pos=(xState2,yState2,zState2),color=color.green,radius=0.01))
#      State2.append(sphere(pos=(xState+normal(0.,eRadius),\
#         yState+normal(0.,eRadius),zState+normal(0.,eRadius)),
#         color=color.green,radius=0.01))
   for n in xrange(0,nIterates):
      for i in xrange(1,nICs):
         State1[i].pos[0],State1[i].pos[1],State1[i].pos[2] =\
            RK3DIntegrator(alpha,wx,wy,wz,State1[i-1].pos[0],State1[i-1].pos[1],State1[i-1].pos[2],\
                           funcX,funcY,funcZ,dt)
         State2[i].pos[0],State2[i].pos[1],State2[i].pos[2] =\
            RK3DIntegrator(alpha,wx,wy,wz,State2[i-1].pos[0],State2[i-1].pos[1],State2[i-1].pos[2],\
                           funcX,funcY,funcZ,dt)
         if i == nICs-1:
            State1[0].pos[0],State1[0].pos[1],State1[0].pos[2] =\
               RK3DIntegrator(alpha,wx,wy,wz,State1[i].pos[0],State1[i].pos[1],State1[i].pos[2],\
                              funcX,funcY,funcZ,dt)
            State2[0].pos[0],State2[0].pos[1],State2[0].pos[2] =\
               RK3DIntegrator(alpha,wx,wy,wz,State2[i].pos[0],State2[i].pos[1],State2[i].pos[2],\
                              funcX,funcY,funcZ,dt)
         
         rate(550)
# end program