/* This program displays two waveforms to see where they */
/* line up.  RWB, Jan. 16, 1996 */

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <graph.h>
#include <string.h>
#include <process.h>
#include <conio.h>
#include <math.h>
#include <errno.h>

#define LWAVE 16384     /* Length of the digitized waveform. */
#define LSEQ 3100       /* Length of one cycle of the m-sequence waveform. */
#define NCCOFF 2000      /* Offset of cross-correlation function. */

void incsamp( char *filename ); /* Increments the sample-number part of the file name. */
void incrun( char *filename );  /* Increments the run-number part of the file name. */
void init( void );      /* Initialization. */
void getwaves( void );
void showaves( void );
void nextfile( void );
void checkey( void );


char xa[LWAVE],xb[LWAVE];
unsigned char uxa[LWAVE],uxb[LWAVE];
int incra=0,incrb;      /* 0 means increment sample number.     /*
			/* 1 means increment run number.        */
			/* 2 means don't increment.             */
FILE *fpin;
char infilea[80],infileb[80];
char *infile1 = "j:\\darmat\\data\\naut\\0797\\09\\wav70000.011";
int ichana = 1; /* Channel A waveform (1 or 2). */
int ichanb = 2; /* Channel B waveform (1 or 2). */
int iquit=0, inewfile=0, iprint=0;
int ninc = 1; /* Amount by which to increment the file number. */
long int icor,lcor,n1,ndisp;
int ifactor;
int yoffa, yoffb;  /* Offsets of two waveforms (down from top). */

int main( int argc, char *argv[] )
{
   struct rccoord oldpos,newpos;
   if( argc <= 1 ) strcpy( infilea, infile1 );
   else strcpy( infilea, argv[1] );

   init();

   iquit = 0;
   while( !iquit ) {
      getwaves();
      showaves();
      nextfile();
   }
   return 0;
}


/* Initialization. */
void init( void )
{
   char instring[80];
   lcor = LWAVE - LSEQ + 1;     /* Maximum length of the cross-correlation array. */
   n1=10;  /* number of waveform points per screen point */
   icor=NCCOFF;  /* Starting value for Ch B offset. */
   ndisp=0;    /* Display shift. */
   yoffa=200;
   yoffb=200;
   ifactor = 20;        /* Multiplicative scale factor. */

	/***** Choosing of input file names. *****/

   printf( "\nChannel A: input file name or return to use %s:\n", infilea );
   gets( instring );
   if( strlen( instring ) > 0 ) strcpy( infilea, instring );

   printf( "\nChannel A waveform (1 or 2): %i", ichana );
   gets( instring );
   if( strlen( instring )>0 ) sscanf( instring, "%i", &ichana );

   printf( "\nIncrement sample number (0), run number (1), or don't increment (2): %i", incra );
   gets( instring );
   if( strlen( instring )>0 ) sscanf( instring, "%i", &incra );

   printf( "\Channel B waveform:  return to use same, + to use following, or enter file name: \n" );
   gets( instring );

   strcpy( infileb, infilea );
   incrb = incra;
   if( strcmp( instring, "+" ) == 0 ) {
      if( incrb == 2 ) {
	 incrb = 0;
	 printf( "\nIncrement sample number (0) or run number (1): %i", incrb );
	 gets( instring );
	 if( strlen( instring )>0 ) sscanf( instring, "%i", &incrb );
      }
      if( incrb==1 ) incrun( infileb );
      else incsamp( infileb );
   }
   else if( strlen( instring )>0 ) {
      strcpy( infileb, instring );
      if( incrb == 2 ) {
	 incrb = 0;
	 printf( "\nIncrement sample number (0) or run number (1): %i", incrb );
	 gets( instring );
	 if( strlen( instring )>0 ) sscanf( instring, "%i", &incrb );
      }
   }

   printf( "\nChannel B waveform (1 or 2): %i", ichanb );
   gets( instring );
   if( strlen( instring )>0 ) sscanf( instring, "%i", &ichanb );

   if( iprint>=0 ) {
      printf( "\nCh A file name:  %s", infilea );
      printf( "\nCh B file name:  %s", infileb );
      printf( "\nichana: %i; ichanb: %i; incra: %i; incrb: %i",
	 ichana, ichanb, incra, incrb );
      getch();
   }
}


/* Increment the run-number part of a file name. */
void incrun( char *filename )
{
   int irun;
   printf( "\nOld file name: %s", filename );
   sscanf( &filename[ strlen(filename)-3 ], "%3d", &irun );
   printf( "\nirun: %i", irun );
   irun+=ninc;
   sprintf( filename,"%3i", irun );
   printf( "\nNew file name: %s", filename );
   return;
}


/* Increment the sample-number part of a file name. */
void incsamp( char *filename )
{
   int isamp;
   char *string1 = "9999";
   sscanf( &filename[ strlen(filename)-8 ], "%4d", &isamp );
   isamp+=ninc;
   sprintf( string1,"%04i", isamp );
   strncpy( &filename[ strlen(filename)-8 ], string1, 4 );
   return;
}


/* Increment the file names. */
void nextfile( void )
{
   if( incra==0 ) incsamp( infilea );
   else if( incra==1 ) incrun( infilea );

   if( incrb==0 ) incsamp( infileb );
   else if( incrb==1 ) incrun( infileb );
}


/* Read in the waveforms from disk files. */
void getwaves( void )
{
   int count,iwave,ibyte;

   /* Channel A input. */
   if( (fpin = fopen( infilea, "rb" )) == NULL ) {
      printf( "\nCould not open file %s", infilea );
      exit( 0 );
   }

   for(iwave = 0; iwave<ichana; iwave++) count = fread( xa, 1, LWAVE, fpin );
   fclose( fpin);

   if( (fpin = fopen( infileb, "rb" )) == NULL ) {
      printf( "\nCould not open file %s", infileb );
      exit( 0 );
   }

   for(iwave = 0; iwave<ichanb; iwave++) count = fread( xb, 1, LWAVE, fpin );
   fclose( fpin);

   if( iprint>=1 ) {
      printf( "\nxa " );
      for( ibyte=0; ibyte<19; ibyte++ ) printf( "%4i", xa[ibyte] );
      printf( "\nxb " );
      for( ibyte=0; ibyte<19; ibyte++ ) printf( "%4i", xb[ibyte] );
      getch();
   }
}

/* Display waveforms, calculate cross correelation function. */
void showaves( void )
{
   struct videoconfig oldvconf, newvconf; /* video config before graphics */
   int defcolor, ntrows;
   long int sum,suma,sumb,i1,ia,i;

   _getvideoconfig( &oldvconf ); /* Save current video configuration. */
   ntrows = _setvideomode( _VRES16COLOR ); /* 640x480, 16 colors (VGA) */
	/* universal color values (all color modes):
	#define _BLACK          0x000000L
	#define _BLUE           0x2a0000L
	#define _GREEN          0x002a00L
	#define _CYAN           0x2a2a00L
	#define _RED            0x00002aL
	#define _MAGENTA        0x2a002aL
	#define _BROWN          0x00152aL
	#define _WHITE          0x2a2a2aL
	#define _GRAY           0x151515L
	#define _LIGHTBLUE      0x3F1515L
	#define _LIGHTGREEN     0x153f15L
	#define _LIGHTCYAN      0x3f3f15L
	#define _LIGHTRED       0x15153fL
	#define _LIGHTMAGENTA   0x3f153fL
	#define _YELLOW         0x153f3fL
	#define _BRIGHTWHITE    0x3f3f3fL  */
   defcolor = _getcolor( );

   inewfile = 0;
   while( !inewfile ) {
      _clearscreen( _GCLEARSCREEN );
      _setcolor( defcolor );
      _moveto( 0,200 );
      _lineto( 739,200 ); /* reference line */

      /* Plot the A waveform. n1 = points per pixel               */
      /* ndisp = offset in waveform points.  */

      _setcolor( 4);
      _moveto( 0, ((int)xa[ndisp+NCCOFF])/2 + yoffa );
      /* Loop over horizontal pixel positions. */
      for( i=0,ia=ndisp; i<640; i++ ) {
	 /* Loop over points sharing one pixel position. */
	 for( i1=0; i1<n1; i1++,ia++ ) {
/*            printf( "\ni: %li; i1: %li; ia: %li; n1: %li; ndisp: %li; yoffa: %i; xa[ia+NCCOFF]: %hi",
	    i, i1, ia, n1, ndisp, yoffa, xa[ia+NCCOFF] ); */
	    _lineto( (int) (i), ((int) xa[ia+NCCOFF])*ifactor/20 + yoffa );
	 }
      }

      /* Plot the B waveform. */
      _setcolor(2);
      _moveto( 0, ((int)xb[ndisp+icor])/2 + yoffb );
      for( i=0,ia=ndisp; i<640; i++ ) {
	 for( i1=0; i1<n1; i1++, ia++ ) {
	    _lineto( (int) (i), ((int) xb[ia+icor])*ifactor/20 + yoffb );
	 }
      }

      /* Calculate the cross-correlation function. */
      sum=0;
      suma=0;
      sumb=0;
      for( i=0; i<LSEQ; i++ ) {
	 sum += (long) ((int) xb[ i + icor ]*(int)xa[i+NCCOFF]);
	 suma+= (long)xa[i+NCCOFF]*(long)xa[i+NCCOFF];
	 sumb+= (long)xb[i+icor]*(long)xb[i+icor];
	 /* if( (i < 0)&&(icor==0) ) printf( "\nxa, xb, sum = %i, %i, %li",
	    (int) xa[i+NCCOFF], (int) xb[i+icor], sum ); */
      }

      _settextposition( 1,1 );
      printf( "\nA: waveform %i from %s", ichana, infilea );
      printf( "   c[%li]: %8li", icor, sum );
      printf( "\nB: waveform %i from %s", ichanb, infileb );
      printf( "   cnorm: %8.5lf", (double)sum/sqrt( (double)suma*(double)sumb ) );
      printf( "\nNow on screen: xa[%li-%li] (red), xb[%li-%li] (green)  rmsA: %3.0f  rmsB: %3.0f",
	 ndisp+NCCOFF,ndisp+NCCOFF+640-1,ndisp+icor,ndisp+icor+640-1,
	 sqrt( suma/(double)LSEQ ), sqrt( sumb/(double)LSEQ ) );
      printf( "\nq=quit; (u,i,j,k,n,m)=shift green and red; (U,I,J,K,N,M)=shift green" );
      printf( "\nz=next WAV file; x=previous one; v=set variables; S,s=shift B; B,b=blowup/down" );
      checkey();
   }
    _setvideomode( oldvconf.mode );
}

/* Wait for a command. */
void checkey( void )
{
   char instring[80], ch1;
   ch1 = (char)getch();
   switch( ch1 ) {
      case 'v':
	 printf( "\nicor is the offset of waveform b relative to waveform a." );
	 printf( "\nicor = %li?  ", icor );
	 gets( instring );
	 if( strlen(instring)>0 ) {
	    sscanf( instring, "%li", &icor );
	    printf( "\nNew value for icor = %li", icor );
	 }
	 printf( "\nndisp is a display offset for both waves." );
	 printf( "\nndisp = %li?  ", ndisp );
	 gets( instring );
	 if( strlen( instring ) > 0 ) {
	    sscanf( instring, "%li", &ndisp );
	    printf( "\nNew value for ndisp = %li", ndisp );
	 }
	 printf( "\nn1 is the number of waveform points to plot per pixel (the compression factor)." );
	 printf( "\nn1 = %li?  ", n1 );
	 gets( instring );
	 if( strlen( instring ) > 0 ) {
	    sscanf( instring, "%li", &n1 );
	    printf( "\nNew value for n1 = %li", n1 );
	 }
	 break;
      case 'i':
	 ndisp+=1;
	 break;
      case 'u':
	 ndisp-=1;
	 break;
      case 'k':
	 ndisp+=10;
	 break;
      case 'j':
	 ndisp-=10;
	 break;
      case 'm':
	 ndisp+=100;
	 break;
      case 'n':
	 ndisp-=100;
	 break;
      case 'I':
	 icor+=1;
	 break;
      case 'U':
	 icor-=1;
	 break;
      case 'K':
	 icor+=10;
	 break;
      case 'J':
	 icor-=10;
	 break;
      case 'M':
	 icor+=100;
	 break;
      case 'N':
	 icor-=100;
	 break;
      case 's':
	 yoffb += 10;
	 break;
      case 'S':
	 yoffb -= 10;
	 break;
      case 'b':
	 ifactor = ifactor*2/3;
	 break;
      case 'B':
	 ifactor = ifactor*3/2;
	 break;
      case 'q':
	 inewfile=1;
	 iquit=1;
	 break;
      case 'x':
	 ninc = -1;
	 inewfile=1;
	 break;
      case 'z':
	 ninc = 1;
	 inewfile = 1;
	 break;
      }
      icor=max((long)0,min(icor,LWAVE-LSEQ+1));
      ndisp=max((long)0,min(ndisp,LWAVE-LSEQ+1));
}