Re: [Visualpython-users] 3d surface plot animation

Now it works! (please try my further developed example below).
However, it is far too slow.
Is there any possibility to speed it up?

How can I add colors depending on the z-value?


Poul Riis


from visual import *
from time import *

class Model:
    def __init__(self):
        self.frame = frame()
        self.model = faces(frame=self.frame, color=color.blue)
        self.vertices = []

    def FacetedTriangle(self, v1, v2, v3, color=color.red):
        """Add a triangle to the model"""
        for v in (v1,v2,v3):
            self.vertices.append(v)

    def FacetedPolygon(self, *v):
        """Appends a planar polygon of any number of vertices to the
model"""
        for t in range(len(v)-2):
            self.FacetedTriangle( v[0], v[t+1], v[t+2] )
            self.FacetedTriangle( v[0], v[t+2], v[t+1] )

    def DrawNormals(self, scale):
        pos = self.model.pos
        normal = self.model.normal
        for i in range(len(pos)):
            arrow(pos=pos[i], axis=normal[i]*scale)

class Mesh (Model):
    def __init__(self, xvalues, yvalues, zvalues):
        Model.__init__(self)

        points = zeros( xvalues.shape + (3,), float )
        points[...,0] = xvalues
        points[...,1] = yvalues
        points[...,2] = zvalues

        for i in range(zvalues.shape[0]-1):
            for j in range(zvalues.shape[1]-1):
                self.FacetedPolygon( points[i,j], points[i,j+1],
                                     points[i+1,j+1], points[i+1,j] )
                
        self.model.pos = self.vertices
        self.model.make_normals()
        self.model.smooth()
        self.model.make_twosided()
        
dxy=0.02 # Length of x-y-intervals
x = arange(-1,1+dxy,dxy)
y = arange(-1,1+dxy,dxy)
d=0.8 # Distance between two oscillators
dh=d/2 # Half that distance
lambdah=d/5 # Wavelength
k=2*pi/lambdah
T=10 # Oscillator period
A=0.025 # Oscillator amplitude
omega=2*pi/T
m1=12 #Number of wavelengths from oscillator 1 to point of constructive
interference
m2=10 #Number of wavelengths from oscillator 2 to point of constructive
interference
n=m1-m2
m=m2
xci=n*(n+2*m)*lambdah**2/2/d
yci=sqrt((d**2-(n*lambdah)**2)*(((n+2*m)*lambdah)**2-d**2))/2/d
print('Point of constructive interference: (',xci,yci,')')
print(sqrt((xci-dh)**2+yci**2)/lambdah)
print(sqrt((xci+dh)**2+yci**2)/lambdah)
N1=m*20
N2=(m+n)*20
a1points=[]
a2points=[]
for i in range(0,N1+1):
    a1points.append((i*(xci-dh)/N1+dh,2*A*cos(i/N1*m*2*pi),1-i*yci/N1))
for i in range(0,N2+1):
    a2points.append((i*(xci+dh)/N2-dh,2*A*cos(i/N2*(m+n)*2*pi),1-i*yci/N2))
def f(x,y,t):
    return
A*(sin(k*sqrt((x-dh)*(x-dh)+y*y)-omega*t)+sin(k*sqrt((x+dh)*(x+dh)+y*y)-omega*t))

source1=sphere(radius=0.02,pos=(dh,0,1),color=color.yellow)
source2=sphere(radius=0.02,pos=(-dh,0,1),color=color.red)
z = zeros( (len(x),len(y)), float )
x,y = x[:,None]+z, y+z

curve(pos=a1points,color=color.red)
curve(pos=a2points,color=color.yellow)

Ncp=20
ncmax=floor(d/lambdah)
if ncmax*lambdah==d:
    ncmax=ncmax-1
for nc in range(1,ncmax+1):
    bhyp=sqrt(-(nc*lambdah/2)**2+(d/2)**2)
    ahyp=bhyp*nc/sqrt(-nc*nc+(d/lambdah)**2)
    print('ahyp,bhyp:',ahyp,bhyp)
    txmax=acosh(1/ahyp)
    tymax=acosh(2/bhyp)
    tmax=txmax
    if tymax<txmax:
        tmax=tymax
    hyppoints=[]
    for i in range(0,Ncp+1):
        t=i*tmax/Ncp
        hyppoints.append((ahyp*cosh(t),2*A,-bhyp*sinh(t)+1))
    curve(pos=hyppoints,color=color.green)
    hyppoints=[]
    for i in range(0,Ncp+1):
        t=i*tmax/Ncp
        hyppoints.append((-ahyp*cosh(t),2*A,-bhyp*sinh(t)+1))
    curve(pos=hyppoints,color=color.green)
    hyppoints=[]
    hyppoints.append((0,2*A,-1))
    hyppoints.append((0,2*A,1))
    curve(pos=hyppoints,color=color.green)
for i in range(0,100):
    surface=Mesh( x, f(x,y-1,i) , y )
    #sleep(0.01)
    surface.model.visible=False
    del surface
    source1.pos=(-dh,-A*sin(omega*i),1)
    source2.pos=(dh,-A*sin(omega*i),1)
i=100
surface=Mesh( x, f(x,y-1,i) , y )