/************* wave.c: Basic Waveguide Code *****************/



#include <stdio.h>

#include <math.h>

#include <fcntl.h>

#include <string.h>



#include <sys/types.h>



#define PREMUL 32768           /* Fixed exponent */

#define MAXN   8192            /* Max # nodes    */

#ifndef PI

#define PI 3.1415926535

#endif



int notcl = 0, nosound = 0;    /* to prevent (time expensive) output */



int e[MAXN], v[MAXN];          /* displacement and velocity arrays   */

int u1 =1024, u2 =1024;        /* combined values of (c/m)*PREMUL*dt */

int n;                         /* number of nodes                    */

int steps;                     /* number of (full wave) iter. steps  */

int fd;                        /* file descriptor for sound file     */

int binoutnode;                /* node to output in sound file       */



init(n)

int n;

{

   int i;

   for (i=0; i<n; i++) {

      e[i] = 0; v[i] = 0;      /* init to zero energy / strain       */

   }

   for (i=1; i<n-1; i++) {

      e[i] = (int) (PREMUL*10000 * 

                (cos(

                   2.0 * PI * (double) i/ (double) (n-2)) 

                 ) 

                 /exp(

                   pow(

                      (8.0 * (double) (i-n/2)

                         /(double) n

                      ),2)

                 ));

      v[i] = 0; /* (int) (PREMUL/2 * cos(2.0 * PI * (double) i/ (double) n) ); */

   }

   return;



   e[0] = 4 * PREMUL;         /* impulse of 10 units                */

   e[32]= -4 * PREMUL;        /* impulse of -10 units               */

   e[64] = 4 * PREMUL;         /* impulse of 10 units                */

   e[96]= -4 * PREMUL;        /* impulse of -10 units               */

   e[128] = 4 * PREMUL;         /* impulse of 10 units                */

   e[160]= -4 * PREMUL;        /* impulse of -10 units               */

   e[192] = 4 * PREMUL;         /* impulse of 10 units                */

   e[224]= -4 * PREMUL;        /* impulse of -10 units               */

}



init_fixend(n)

int n;

{

   int i;

   for (i=0; i<n; i++) {

      e[i] = 0; v[i] = 0;      /* init to zero energy / strain       */

   }

   for (i=0; i<n; i++) {

      e[i] = (int) (PREMUL*30000 * 

                (sin(

                  2.0 * PI * ((double) (i))/ (double) (n-1)) 

                ));

      /* e[i] += (int) (PREMUL*2000 * 

                (sin(

                   16.0 * PI * ((double) (i))/ (double) (n-1)) 

                )); */

      v[i] = 0; /* (int) (PREMUL/2 * cos(2.0 * PI * (double) i/ (double) n) ); */

   }

   return;

}



init_string(n)

int n;

{

   int i;

   for (i=0; i<n; i++)

      e[i] = v[i] = 0;

   e[61] =  2000*PREMUL;

   e[62] =  5000*PREMUL;

   e[63] =  7000*PREMUL;

   e[64] =  8000*PREMUL;

   e[65] =  7000*PREMUL;

   e[66] =  5000*PREMUL;

   e[67] =  2000*PREMUL;

}



inittcl(n)

int n;

{

   int i;

   if (notcl) return;

   printf("if {[winfo exists .c] == 0} {canvas .c; pack .c -expand y -fill both}\n");

   printf("proc graph {n} {global a; for {set i 0} {$i<$n} {set i [expr $i+1]} {.c create line [expr $i] 300 [expr $i] [expr 300-$a($i)/90] -width 1 -fill red}\n.c create line 0 300 [expr  $n] 300 -fill green}\n");

   printf("proc graphb {n} {global a; for {set i 1} {$i<$n} {set i [expr $i+1]} {.c create line [expr $i] [expr 300-$a([expr $i-1])/90] [expr $i] [expr 300-$a($i)/90] -width 1 -fill purple}\n.c create line 0 300 [expr $n] 300 -fill blue}\n");

   printf("set nnodes %d\n",n);

   printf(".c delete all\n");

   for (i=0; i<n; i++)

      printf("set a(%d) %d\n",i,e[i]/PREMUL);

   printf("graphb $nnodes\n update \n after 2000\nset vt 2\n");

}



initbin(n,s)                     /* could add a .wav header here   */

int n;

char *s;

{

   if (nosound) return;

   fd = open(s,O_WRONLY|O_CREAT|O_TRUNC|O_BINARY,S_IWUSR);

}



outputtcl(n)

int n;

{

   int i;

   if (notcl) return;

   for (i=0; i<n; i++)

      printf("set a(%d) %d\n",i,e[i]/PREMUL);

   printf(".c delete all\n");

   printf("graph $nnodes\n update\n after $vt\n");

}



outputbin(i)

int i;

{

   char *cp;

   int t;

   if (fd > 0) {

      t = e[i]/PREMUL;

      cp = (char *) &t;

      write(fd,cp,1);

      write(fd,cp+1,1);

   }

}



propagate(n)                   /* Compute new velocity and           */

int n;                         /* displacements once for all nodes.  */

{                              /* (more optimization possible here)  */

   register int i;

   for (i=0; i<n; i++)

      v[i] += (u1 * (e[(i-1+n)%n]+e[(i+1)%n]-2*e[i]) ) / PREMUL;

   for (i=0; i<n; i++)

      e[i] += (u2 * v[i] ) / PREMUL;

}



propagate_fixend(n)            /* Compute new velocity and           */

int n;                         /* displacements once for all nodes.  */

{                              /* (more optimization possible here)  */

   register int i;

   for (i=1; i<n-1; i++) {

      /* v[i] += (u1 * (e[i-1]+e[i+1]-2*e[i]) )/PREMUL; */

      v[i] += e[i-1];

      v[i] += e[i+1];

      v[i] -= e[i];

      v[i] -= e[i];

   }

   for (i=1; i<n-1; i++) {

      /* e[i] += (u2 * v[i] ) / PREMUL; */

      e[i] += v[i];

   }

      e[1] *= 5;

      e[1] /= 10;

      e[n-2] *= 5;

      e[n-2] /= 10;

}



simulate(n,tn)

int n, tn;

{

   int t;

   for (t=0; t<tn; t++) {

      /* propagate(n); */

      propagate_fixend(n);

      output(t,n);

      outputbin(binoutnode);

   }

}



output(t,n)

int t,n;

{

   int i;

   outputtcl(n);

   return;





   printf("t=%3d|",t);

   for (i=0; i<n; i++)

      printf(" %3d",e[i]);

   printf("\n");

}



closeup()

{

   if (fd > 0) close(fd);

}





process_args(argc,argv)

int argc;

char *argv[];

{

   int i;

   if ((argc < 3) || ( ((argc-3)%2) != 0 )) {

      printf("Usage %s <#nodes> <#steps> [<option> <value>]...\n",

         argv[0]);

      printf("option = binfile  <file.raw>  | tcl <1|0>        | \n");

      printf("         initstate <file.asc> | outnode <node #>   \n");

      exit(-1);

   }

   sscanf(argv[1],"%d",&n);

   sscanf(argv[2],"%d",&steps);

   for (i=3; i<argc; i++) {

      if (strcmp(argv[i],"binfile") == 0) {

           initbin(n,argv[i+1]); i++;

      } else 

      if (strcmp(argv[i],"tcl") == 0) {

         if (argv[i+1][0] == '1') notcl=0; else notcl=1; i++;

      } else

      if (strcmp(argv[i],"initstate") == 0) {

         printf("Warning: initstate not yet implemented\n"); i++;

      } else

      if (strcmp(argv[i],"outnode") == 0) {

         sscanf(argv[i+1],"%d",&binoutnode); i++;

      } else {

         printf("Error: unknown option %s\n",argv[i]);

         exit(-1);

      }

   } 

}



main(argc, argv)

int argc;

char *argv[];

{

   notcl = 1;

   process_args(argc,argv);

   /* init_fixend(n); */ init_string(n);  /* init(n); */

   inittcl(n);

   simulate(n,steps);

   closeup();

   return 0;

}