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


int _stklen=100000;

#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       */
FILE *fp;

int el[MAXN], ew[512], flag;


/*
 * this example uses DPMI to initialize Mode 13h video mode and 
 * draws three lines on the screen in different colors
 */
#include <conio.h>
#include <dpmi.h>

#define PIXEL(x,y) (videobase + (x) + (y) *320)
/* this will point to the text area of video memory */
char * videobase;

void GetVidBase(void)
{
  ULONG base;
  SELECTOR temp;
  /* get the data segment selector */
  asm mov [temp],ds
  /* now get its base address */
  dpmi_get_sel_base(&base,temp);
  /* video base is address zero + 0xa0000 for VGA graphics area */
  videobase = (char *)(-base + 0xa0000);

	/* notice we do NOT change the DS selector limit.  Our segments have max
   * possible length so you can wrap around.
   */
}
void ModeSet(int mode)
{
  DPMI_REGS regs;
  regs.b.ah = 0x0;
  regs.b.al = mode;
  regs.h.flags = 0x72;
  regs.h.ss = 0;
  regs.h.sp = 0;
  dpmi_simulate_real_interrupt(0x10,&regs);
}




init(n)
int n;
{
   int i;
   double fl, h;

   for (i=0; i<n; i++) {
      e[i] = 0; v[i] = 0;      /* init to zero energy / strain       */
   }

   h = 3*PI;
   for (i=1; i<n-1; i++) {
      fl = (double) i;
	      fl = fl - (double) n/2.0 ;
      fl = fl /(double) n;
/*      fl = pow( (8.0 * fl ),2.0); */
      fl=1.5*8*fl*fl;
      e[i] = (int) (PREMUL*5000 *  (
          sin(h/2.0+ 1.0  * PI * (double) i/ (double) (n-2)) + 
   (1/3 )*sin(h/2.0+ 3.0  * PI * (double) i/ (double) (n-2)) +
   (1/5 )*sin(h/2.0+ 5.0  * PI * (double) i/ (double) (n-2)) +
   (1/7 )*sin(h/2.0+ 7.0  * PI * (double) i/ (double) (n-2)) +
   (1/9 )*sin(h/2.0+ 8.0  * PI * (double) i/ (double) (n-2)) +
   (1/11)*sin(h/2.0+ 11.0 * PI * (double) i/ (double) (n-2)) +
   (1/13)*sin(h/2.0+ 13.0 * PI * (double) i/ (double) (n-2)) +
   (1/2 )*sin(h/2.0+ 2.0  * PI * (double) i/ (double) (n-2)) +
   (1/4 )*sin(h/2.0+ 4.0  * PI * (double) i/ (double) (n-2)) +
   (1/6 )*sin(h/2.0+ 6.0  * PI * (double) i/ (double) (n-2)) +
   (1/8 )*sin(h/2.0+ 8.0  * PI * (double) i/ (double) (n-2)) +
   (1/12)*sin(h/2.0+ 12.0 * PI * (double) i/ (double) (n-2)) +
   (1/14)*sin(h/2.0+ 14.0 * PI * (double) i/ (double) (n-2)) 
           ) / exp( fl )
      );
      v[i] = 0; /* (int) (PREMUL/2 * cos(2.0 * PI * (double) i/ (double) n) ); */
   }

   return 0;

   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               */

   return 0;
}



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*300 * (sin( 2.0 * PI * ((double) (i))/ (double) (n-1)) ));


      v[i] = 0; 
   }

   return 0;

}



init_string(n)

int n;

{

   int i;

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

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

   e[1] =  200*PREMUL;

   e[2] =  500*PREMUL;

   e[3] =  700*PREMUL;

   e[4] =  800*PREMUL;

   e[5] =  700*PREMUL;

   e[6] =  500*PREMUL;

   e[7] =  200*PREMUL;

   return 0;
}



inittcl(n)
int n;
{
   int i;

   if (notcl) return 0;

   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");

   return 0;
}



int initbin(n,s)                     /* could add a .wav header here   */
int n;
char *s;
{

   if (nosound) return 0;

   fp = (FILE *) fopen(s,"wb");
   if (fp == NULL) return(-1);

/*   fd = open(s,_F_RDWR|_F_BIN); */

   return 0;
}



int outputtcl(n)

int n;

{

   int i;

   if (notcl) return 0;

   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");

   return 0;
}



int outputbin(i)

int i;

{

   char *cp;

   int t;

/*      printf("%d ",e[i]); */
      t = e[i] / (20*PREMUL); t+= 128;
/*      printf("%d\n",t);  */
   if (fp  != NULL)
      fwrite(&t,1,1,fp);

   return 0;
   if (fd > 0) {

      t = e[i]/PREMUL;

      cp = (char *) &t;

/*      write(fd,cp,1);

      write(fd,cp+1,1);
*/
   }

   return 0;
}



int 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;

   return 0;
}



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

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

{                              /* (more optimization possible here)  */

   /* register int i; */
   int i, it;
   double fl;

   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];

   }

   
   it = e[1]; 
   fl = it; fl *= 0.95; it = (int) fl;
   e[1] = it;
   it = e[n-2]; 
   fl = it; fl *= 0.93; it = (int) fl;
   e[n-2] = it;

   it = e[2]; 
   fl = it; fl *= 0.98; it = (int) fl;
   e[2] = it;
   it = e[n-3]; 
   fl = it; fl *= 0.97; it = (int) fl;
   e[n-3] = it;

   it = e[3]; 
   fl = it; fl *= 0.995; it = (int) fl;
   e[3] = it;
   it = e[n-4]; 
   fl = it; fl *= 0.99; 

   fl =fl / (80*PREMUL); fl = fl / 1.0;
   if (fl > 1.0) { fl=sqrt(sqrt(fl)); flag=1; 
/*      printf("******"); */
   } else flag = 0;
   fl =fl * (80*PREMUL ); fl = fl * 1.0;
   it = (int) fl;
   e[n-4] = it;


/*
      e[1] *= 9;
      e[1] /= 10;
      e[n-2] *= 9;
      e[n-2] /= 10;
*/

   return 0;
}





void 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");

}

output_vid(n)
{
   int i,j;
   for (i=0; i< n; i++)
      *PIXEL(i,128-el[i]) = BLACK;

   for (i=0; i< n; i++)
      el[i]=(int) e[i]/(80*PREMUL);

   if (flag == 1) j = BLUE; else j = YELLOW;
   for (i=0; i< n; i++)
      *PIXEL(i,128-el[i]) = j;


   for (i=0; i< 255; i++) {
      *PIXEL(i,128-ew[i]) = BLACK;
      *PIXEL(i,128-(ew[i]=ew[i+1]) ) = RED;
   }

   ew[255]=(int) e[binoutnode]/(80*PREMUL);

}

int simulate(n,tn)
int n, tn;
{
   int t, i;
   for (t=0; t<tn; t++) {
      /* propagate(n); */
      propagate_fixend(n);
      output_vid(n);
      output(t,n);
      outputbin(binoutnode);
      i=0;
      if (kbhit()) {
         if (getch() != ' ') { while (!kbhit()); getch(); } 
         else while (!kbhit());
      }
   }
   return 0;
}



void closeup()

{

   if (fp != NULL) fclose(fp);
/*   if (fd > 0) close(fd); */

   return;
}


init_vid()
{
  int i;
  /* get a base pointer to video mem */
  GetVidBase();

  /* 320x200x16 graphics mode */
  ModeSet(0x13);
  for (i=0; i< MAXN; i++) el[i]=0;
  for (i=0; i< 256; i++) ew[i]=0;

/*  for (i=0; i< 8*256; i++) *(videobase+i*128)=(i);
  for (i=0; i< 8*256; i++) *(videobase+1024+i+16)=1<<((i/8) % 8);
  for (i=0; i< 32; i++) *(videobase+18000-4096+i/2+300)=(1<<(i%8));
key_wait(); 
exit(0); */
}

key_wait() {
  /* wait for a key to be hit */
  while (!kbhit());
  getch();

  /* text mode */
  ModeSet(3);
}

void 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);
      }
   } 
   return ;
}


int main(argc, argv)

int argc;

char *argv[];

{

   notcl = 1;

   process_args(argc,argv);

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

   inittcl(n);

   simulate(n,steps);

   key_wait();
   closeup();

   return 0;

}