#declare X_RES = 35;  // number of points in x direction
#declare Y_RES = 35;  // number of points in y direction
#declare C_SQUARED = 0.3;  // c^2 constant in wave equation
#declare DT = 0.1;   // the time step


// declare the arrays for the wave height and velocity
// and two for temporary values
#declare waveHeight    = array[X_RES][Y_RES];
#declare waveVelocity  = array[X_RES][Y_RES];
#declare waveHeightT   = array[X_RES][Y_RES];
#declare waveVelocityT = array[X_RES][Y_RES];


// initialise the arrays to zero
#local X=0;
#while(X<X_RES)

  #local Y=0;
  #while(Y<Y_RES)

    #declare waveHeight[X][Y]   = 0;
    #declare waveVelocity[X][Y] = 0;

  #local Y=Y+1;
  #end

#local X=X+1;
#end

#declare waveHeightT   = waveHeight;
#declare waveVelocityT = waveVelocity;












// perform N iterations of the wave equation
#local ITERATION=0;
#while(ITERATION<=300)

  #local X=1;
  #while(X<X_RES-1)

    #local Y=1;
    #while(Y<Y_RES-1)

      // first calculate the second derivatives in the X and Y direction
      #local XD = waveHeight[X+1][Y] - 2*waveHeight[X][Y] + waveHeight[X-1][Y];
      #local YD = waveHeight[X][Y+1] - 2*waveHeight[X][Y] + waveHeight[X][Y-1];
      
      // now calculate the acceleration
      #local A = C_SQUARED * (XD + YD);
      
      // and add it to the velocity using simple numerical integration
      #local waveVelocityT[X][Y] = waveVelocityT[X][Y] + A * DT;
      
      // provide some damping
      #local waveVelocityT[X][Y] = waveVelocityT[X][Y] * 0.995;
            
      // add the velocity to the position
      #local waveHeightT[X][Y] = waveHeightT[X][Y] + waveVelocity[X][Y] * DT;

    
    #local Y=Y+1;
    #end

  #local X=X+1;
  #end

  #local waveHeight = waveHeightT;
  #local waveVelocity = waveVelocityT;

  // make the center point move up and down to create the wave
  #declare waveHeight[X_RES/2][Y_RES/2] = 0.25*sin(ITERATION/500*4.5*pi);


#local ITERATION=ITERATION+1;
#end










// plot the array as points joined by cylinders
light_source{ 100 color rgb 1 }
camera {
   perspective
   location <1,2,.5>
   up z
   right x*4/3
   angle 40
   sky z
   look_at <.5,.5,0>
}

union{

    #local X=1;
    #while(X<X_RES)
    
      #local Y=1;
      #while(Y<Y_RES)

         cylinder{ <X/X_RES,Y/Y_RES,waveHeight[X][Y]> 
                   <(X-1)/X_RES,Y/Y_RES,waveHeight[X-1][Y]> 
                   0.1/X_RES }    
    
         cylinder{ <X/X_RES,Y/Y_RES,waveHeight[X][Y]> 
                   <X/X_RES,(Y-1)/Y_RES,waveHeight[X][Y-1]> 
                   0.1/X_RES }    
    
    
      #local Y=Y+1;
      #end
    
    #local X=X+1;
    #end


    pigment{
      gradient z
      color_map{
       [0.0 color rgb <0,1,1>]
       [0.5 color rgb <1,1,0>]
       [1.0 color rgb <1,0,0>]
       }
       translate -z/2
       scale 0.5
    
    
       }
    finish{ specular 0.5 roughness 0.01 }
}