#include <stdio.h> 
#include <stdlib.h>
#include <string.h>
#include <math.h> 
#define Nx 300    /* the number of intervals  */
#define Ny 300    /* rem more points than intervals!  */
#define Nframes 240 // number of final frame 
#define Ndt 60      // number of subintervals per frame

int main() {
    float u0[Nx+1][Ny+1],u1[Nx+1][Ny+1],u2[Nx+1][Ny+1],v0[Nx+1][Ny+1];
    int i,j,result,iter,array[Nx+1][Ny+1];
    int irows=Nx+1;
    int icols=Ny+1;
    int iframe=0;
    char namebuf[13];
    float dx=1.0; 
    float c=2;  /* node spacing and wave speed (about 1 node per frame) */
    float dt,lambda,phase_wave;
    float bdx,bdy;
    int ndt = Ndt;
    float center=(float)Nx/2.0;
    float br=center-2.0;
    dt=1/(float)Ndt;      //this makes interframe time interval = 1   
    lambda = (c*c*dt*dt/(dx*dx));
     
//initialize, try to keep "disturbance height" on order of +/- unity    
    //initial velocity and position
    for(i=0 ; i<=Nx ; i++) {
        for(j=0 ; j<=Ny ; j++) {
            u0[i][j]=0.0;
            u1[i][j]=0.0;
            v0[i][j]=0.0;    
        }
    }
/*
    //softened initial disturbance
     for(i=-4 ; i<=4 ; i++) {
        for(j=-4 ; j<=4 ; j++) {
        u0[Nx/3+i][Ny/4+j]=exp(-(float)(i*i+j*j)/2.0);  //singular disturbance
        }
    }
*/
 
  
    //step to next position
        //u1 = u0 +v0*dt +(lambda/2)*(D^2) u (lambda*dt/6)*(D^2)v
            //algorithm does not include edges (which are usually the boundaries)
    for(i=1 ; i<Nx ; i++) {
        for(j=1 ; j<Ny ; j++) {
            u1[i][j]=
                u0[i][j]+v0[i][j]*dt
                +(.5*lambda)*(u0[i+1][j]+u0[i-1][j]+u0[i][j+1]+u0[i][j-1]-4.0*u0[i][j])
                +(lambda*dt/6.0)*(v0[i+1][j]+v0[i-1][j]+v0[i][j+1]+v0[i][j-1]-4.0*v0[i][j]);    
        }
    }

            //enforce boundary conditions
                //clamped circle? boundaries
            for(i=0 ; i<=Nx ; i++) {
                bdx=(float)i-center;
                bdy=br*br-bdx*bdx;
                if(bdy<0){
                    bdy=0.;
                    }
                bdy=pow(bdy,.5);
                for(j=0 ; j<=center-bdy ; j++) {
                    u2[i][j]=0;
                    }
                for(j=center+bdy; j<=Ny ; j++) {
                    u2[i][j]=0;
                    }
                }
        
        
    //output this frame
        for(i=0 ; i<irows ; i++) {
            for(j=0 ; j<icols ; j++) {
                (int)array[i][j]=(int)(20000.*u0[i][j]+65535./2.0);
            }
        }
        sprintf( namebuf, "wdm%.4i.tga", iframe );
        result=tgaout(irows, icols, namebuf, &array[0][0]);
        printf("Writing  frame %i\n",iframe);

// cycle through remaining frames
    for(iframe=1 ; iframe<=Nframes ; iframe++) {

    //evolve to next frame
        for(iter=0 ; iter<ndt ; iter++) {
        
            //step to next position
                //u2 = 2*u1 -u0 +(lambda)*(D^2) u
                    //algorithm does not include edges (which are usually the boundaries)
            for(i=1 ; i<Nx ; i++) {
                for(j=1 ; j<Ny ; j++) {
                    u2[i][j]=2*u1[i][j]-u0[i][j]+lambda*(u1[i+1][j]+u1[i-1][j]+u1[i][j+1]+u1[i][j-1]-4.0*u1[i][j]);    
                }
            }
            
            //enforce boundary conditions
                //clamped diamond boundaries
            for(i=0 ; i<=Nx ; i++) {
                bdx=(float)i-center;
                bdy=br*br-bdx*bdx;
                if(bdy<0){
                    bdy=0.;
                    }
                bdy=pow(bdy,.5);
                for(j=0 ; j<=center-bdy ; j++) {
                    u2[i][j]=0;
                    }
                for(j=center+bdy; j<=Ny ; j++) {
                    u2[i][j]=0;
                    }
                }
        
            //add a few oscillations from a little force oscillator
            phase_wave=.5*(iframe+(float)iter/(float)ndt);
            if(phase_wave<=3*3.14159){
                u2[Nx/2][Ny/4]=.15*sin(phase_wave);
            }
                    //exchange for next iteration
            for(i=0 ; i<=Nx ; i++) {
                for(j=0 ; j<=Ny ; j++) {
                    u0[i][j]=u1[i][j];    
                    u1[i][j]=u2[i][j];    
                }
            }
                
    }    
    
    
    
     //output this frame
        for(i=0 ; i<irows ; i++) {
            for(j=0 ; j<icols ; j++) {
                (int)array[i][j]=(20000.*u0[i][j]+65535./2.0);
            }
        }
        sprintf( namebuf, "wdm%.4i.tga", iframe );
        result=tgaout(irows, icols, namebuf, &array[0][0]);
        printf("Writing  frame %i\n",iframe);

    }
    






    
    return 0;
}




/*
 * 
 * 
 * adapted by mrg3@psu.edu from
 * 
 * asc2tga.cpp
 * converts ascii matrix to tga heightfield
 * Copyright (C) 1999 by Ortwin Glueck
 * Released under the GNU General Public License (GPL)
 * This is free software. No warranty.
 */

#include <stdio.h> 
#include <stdlib.h>

unsigned char tgaheader[]={0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,24,32}; 

int tgaout(int irows, int icols, char *fname, int *data)
{
 FILE *out;
 unsigned int value,i,j;
 unsigned short rows,cols;
 unsigned char BGR[3];
 unsigned short *target;

 out=fopen(fname,"wb");
 if (out==NULL)
 {
  printf("Could not open output file\n");
  exit(1);
 }
 
 target=(unsigned short*)&tgaheader[12];
 (unsigned short)rows=irows;
 (unsigned short)cols=icols;
 target[0]=rows;
 target[1]=cols;
 fwrite(tgaheader,18,1,out);
 
 BGR[0]=0; //blue is not used
 for (i=0;i<rows;i++)
  for (j=0;j<cols;j++)
  {
   value = *(data+i*cols+j);
   if(value<0){
        value=0;
   }
   if(value>65535){
        value=65535;
    }

   //printf("Writing %u value\n",value);
   *((unsigned short*)&BGR[1])=value;
   fwrite(BGR,3,1,out);
  }
    fclose(out);
    return 0;

}