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

[Shawn D. <sh...@ty...>]

 

Poul wrote...

> Yes, I know, but even the original example 'faces_heightfield' is 

> too complicated for me to see how I can add these colors. Can 

> someone give me just a little hint?

 

 

Below is some added color.  Search for “s.dube”

 

-Shawn Dube

 

 

 

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)

 

# function added by s.dube

def PosArrayToHueColorArray( aPos, index ):

   # aPos is an array of positions

   # index chooses from which axis (x,y,z) to base color

   aPos= aPos.copy()

   v= hsplit(aPos,3)[index]

   v -= v.min()

   v *= 5.0 / v.max()

   # red

   r= v+1.0

   r %= 6.0

   r -= 3.0

   r= abs(r)

   r -= 1.0

   r= r.clip(0.0,1.0)

   # green

   g = v+3.0

   g %= 6.0

   g -= 3.0

   g= abs(g)

   g -= 1.0

   g= g.clip(0.0,1.0)

   # blue

   b= v

   b += 5.0

   b %= 6.0

   b -= 3.0

   b= abs(b)

   b -= 1.0

   b= b.clip(0.0,1.0)

   # wrap-up

   return hstack((r,g,b))

 

 

class Mesh (Model):

   def __init__(self, xvalues, yvalues, zvalues):

       global g

       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.color= PosArrayToHueColorArray( self.model.pos, 1 ) # added by s.dube

       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 )