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

#include        <signal.h>

#include        "\djgpp\allegro\allegro.h"

#include        "..\include\globals.h"

extern double trcorrtog;

void histogram(DWELL *dwell, RADAR *radar, float *prods, PARAMDAT *dat, int param);
void proddisplay(DWELL *dwell, RADAR *radar, float *prods, PARAMDAT *dat, int param);
void displayfft(DWELL *dwell, float array[],int size,double filter);
void iqphase(DWELL *dwell, float array[],int size,double filter);
void idisplayts(DWELL *dwell, float array[],int size,double filter);
void qdisplayts(DWELL *dwell, float array[],int size,double filter);
void drawgraphscale(char *xbuf,double xlo, double xhi, char *ybuf, double ylo, double yhi);
void uncursor(double,double);
void cursor(double,double);
void    helpindex(void);


#define STEP    1638.4

#define XTEXT   460
#define YTEXT   24
#define XTITLE  180
#define YTITLE  448
#define THEIGHT 12

float           *array;
double          TSSCALE;

int     PARMNUM;  /* number of bytes per gate */

static  PARAMDAT        Dat[16];
static  int   Lastparm = -1;

int     getseres(float *,int,float *,int);

double DACVOLTAGE = -4.2;

main(int argc,char *argv[])
   {
   FILE *fp;
   DWELL        *dwell;
   CONFIG       *config;
   PIRAQ        *piraq;
   RADAR        *radar;
   DISPLAY      *disp;
   float        *prods;
   double       isum,qsum,v,wp;
   int  i,n,a,lasta,*g,I,Q,p,line,parameter=0;
   float *temp;   
   int  gate0,oldgates,oldmode,Click=0;
   char c,fname[80],buf[80];

   a = lasta = p = gate0 = 0;
   TSSCALE = 1.0;

   /* get either the number of samples, or the filename */
   strcpy(fname,"config");
   n = 256;
   if(argc > 1)
     {
     if(isdigit(argv[1][0]))
	n = atoi(argv[1]);
     else   
	strcpy(fname,argv[1]);   
     }
   
   config = (CONFIG *)malloc(sizeof(CONFIG));
   if(!config)  {printf("not enough mem for config structure\n"); exit(0);}
   
   /* read in the config file (always do this first) */
   readconfig(fname,config);

   config->timeseries = 1;  /* turn on the time series */
   
   /* allocate array for timeseries */
   array = (float *)malloc(4 * 2000 * sizeof(float));
   if(!array)  {printf("not enough mem for time series\n"); exit(0);}
   for(i=0; i<8000; i++)        array[i] = 0.0;   

   /* allocate a piraq structure */
   piraq = (PIRAQ *)malloc(sizeof(PIRAQ));
   if(!piraq)  {printf("not enough mem for piraq structure\n"); exit(0);}
   
   /* allocate a product array */
   prods = (float *)malloc(16 * 1000 * sizeof(float));
   if(!prods)  {printf("not enough mem for products array\n"); exit(0);}
   
   /* allocate a radar structure */
   disp = (DISPLAY *)malloc(sizeof(DISPLAY));
   if(!radar)  {printf("not enough mem for radar structure\n"); exit(0);}
   
   /* allocate a radar structure */
   radar = (RADAR *)malloc(sizeof(RADAR));
   if(!radar)  {printf("not enough mem for radar structure\n"); exit(0);}
   
   /* allocate a dwell structure */
   dwell = (DWELL *)malloc(sizeof(DWELL));
   if(!dwell)  {printf("not enough mem for dwell structure\n"); exit(0);}
   
   /* initialize all the pointers to various piraq values */
   setpiraqptr(config,piraq);   
   
   /* do an fft and get the malloc out of the way */
   for(i=0; i<2*n; i++) array[i] = 1.0;
   fft(array,n);  /* allocate and calculate the fft lookup table for zeta */
   
   /* load the dsp code into dual port ram */
   initsystem(config);   /* map dual port ram and enable it. turn timer,dsp off */
   cleardpram(config,piraq);         /* clear dual port ram */
   load(config,piraq);  /* load code into dpram */
   
   readradar("",radar);
   
   /* set values in dpram to command the piraq to the config state */ 
   setdpram(config,piraq);

   v = -5.8; /* readv(config);*/   /* read the last dac voltage */

   /* read in the display configuration file */
   readdisplay(fname,disp);
   
   signal(SIGFPE,SIG_IGN);      /* try to ignore floating point exceptions */

   timerset(config);
   if(!start(config))     exit(-1);

   /* make sure there is data in the PIRAQ */
   if(!wait(config,piraq))  exit(-1);
   if(!wait(config,piraq))  exit(-1);

   filldwell(dwell,piraq,config);

   initparmdat(Dat,radar,dwell,disp);

   for(i=0; i<15; i++) /* make the scale and offset suitable to ascope */
      {
      Dat[i].scale = 2.0 / (Dat[i].high - Dat[i].low);
      Dat[i].offset = -(Dat[i].high + Dat[i].low) / (Dat[i].high - Dat[i].low);
      }

   switch(dwell->header.dataformat) 
      {
      case DATA_SIMPLEPP:   PARMNUM = 3 * sizeof(float); break;
      case   DATA_DUALPP:   PARMNUM = 6 * sizeof(float); break;
      case     DATA_POL1:   PARMNUM = 6 * sizeof(float); break;
      case   DATA_POLYPP:   PARMNUM = 5 * sizeof(float); break;
      case 6:               PARMNUM = 3 * sizeof(short); break;    
      }

   glibinit();  /* do another malloc for glib */   

   setcolor(15);
   
   helpindex();

   while(1)
      {
      if(!wait(config,piraq))  break;

      /* if not doing expanded gate time series */
      if(a < 6)        filldwell(dwell,piraq,config);
      
      c = kbhit() ? toupper(getch()) : 1;       /* get key (1 is benign) */
      if(c == 0)   /* if function key */
	 {
	 c = getch() - 0x3b;    /* which function key (0 - 9) */
	 if(c < 6)  parameter = c;
	 }

      if(isdigit(c)) 
	 if(c - '1' < 9)    a = c - '1';
	 
      if(c == 'X') trcorrtog = trcorrtog == 1? 0: 1; /* Toggle TR correction */ 
      
      if(c == 'Q' || c == 27)     break;  /* get the hell out */
      
      /* print info */
      line = YTEXT;
      sprintf(buf,"TSeries Gate= %4d",config->tsgate);   grputs(XTEXT,line+=THEIGHT,buf);
      sprintf(buf," Dac Voltage= %5.2f",v);   grputs(XTEXT,line+=THEIGHT,buf);
      sprintf(buf,"  Trig Delay= %4d",config->delay);   grputs(XTEXT,line+=THEIGHT,buf);
      sprintf(buf,"  Full Scale= +-%5f",2048.0/TSSCALE);   grputs(XTEXT,line+=THEIGHT,buf);
      sprintf(buf,"Clutter Filt= %s",config->clutterfilter?"ON ":"OFF");   grputs(XTEXT,line+=THEIGHT,buf);
      sprintf(buf,"Phase Correc= %s",config->pcorrect?"ON ":"OFF");   grputs(XTEXT,line+=THEIGHT,buf);
      sprintf(buf,"Noise Correc= %6.1f",radar->noise_power);   grputs(XTEXT,line+=THEIGHT,buf);

      if(a > 5 && !gate0)  /* if 0mode requested but not executed */
	 {
	 gate0 = 1;
	 stop(config);
	 oldgates = config->gates;
	 oldmode = config->gate0mode;
	 config->gates = 1;
	 config->gate0mode = 1;
	 setdpram(config,piraq);
	 timerset(config);
	 if(!start(config))     {tmode();  exit(-1);}
	 continue;
	 }

      if(a < 6 && gate0)   /* if 0mode off requested but not executed */
	 {
	 gate0 = 0;
	 stop(config);
	 config->gates = oldgates;
	 config->gate0mode = oldmode;
	 setdpram(config,piraq);
	 timerset(config);
	 if(!start(config))     {tmode();  exit(-1);}
	 continue;
	 }

      /* time series display scale */
      if(c == 'U')     TSSCALE *= 2.0;
      if(c == 'D')     TSSCALE *= 0.5;

      if(c == 'N')    radar->noise_power += 0.1;
      if(c == 'M')    radar->noise_power -= 0.1;
     
      /* afclock */
      if(c == 'A' && config->pcorrect)        
	 {
	 tmode();
	 afclock(config,piraq); 
	 gmode();
	 }

      /* dac output voltage */
      if(c == ']')     {dac(v+=0.05); p = 1;}
      if(c == '[')     {dac(v-=0.05); p = 1;}
      
      if((mouse_b & 1) && !Click)
	{      
	p = 1;
	config->tsgate = dwell->header.gates * (mouse_x - 5) / 400;
	setdpram(config,piraq);
	}
      Click = mouse_b & 1;

      if(c == '=')      /* increment time series gate */
	 {
	 p = 1;
	 stop(config);
	 config->tsgate += 1;
	 setdpram(config,piraq);
	 timerset(config);
	 if(!start(config))     {tmode();  exit(-1);}
	 continue;
	 }

      if(c == '-')     /* decrement time series gate */
	 {
	 p = 1;
	 stop(config);
	 config->tsgate -= 1;
	 setdpram(config,piraq);
	 timerset(config);
	 if(!start(config))     {tmode();  exit(-1);}
	 continue;
	 }

      if(c == '.')     /* increment delay */
	 {
	 p = 1;
	 stop(config);
	 config->delay += 1;
	 setdpram(config,piraq);
	 timerset(config);
	 if(!start(config))     {tmode();  exit(-1);}
	 continue;
	 }

      if(c == ',')     /* decrement delay */
	 {
	 p = 1;
	 stop(config);
	 config->delay -= 1;
	 setdpram(config,piraq);
	 timerset(config);
	 if(!start(config))     {tmode();  exit(-1);}
	 continue;
	 }

      if(c == 'C')     /* toggle clutter filter */
	 {
	 config->clutterfilter = 1 - config->clutterfilter;
	 setdpram(config,piraq);
	 }

      if(c == 'O')     /* toggle phase correction */
	 {
	 config->pcorrect = 1 - config->pcorrect;
	 setdpram(config,piraq);
	 }
      
       if(a == 6)     /* 16 MHz time series */
	  {
	  g = piraq->tsptr[*(piraq->flag) & 1];
	  for(i=0; i<config->pulsewidth * 2; i++)
	     {
	     array[i] = TSSCALE * toieee(*g--) / (double)0x2000000;
	     }
	  avedisp(array,config->pulsewidth * 2,1.0);
	  grputs(XTITLE,YTITLE,"Gate 0 time series");
	  continue;  /* go back to begining */
	  }
	  
       if(a == 7)     /* 16 MHz fft */
	  {
	  g = piraq->tsptr[*(piraq->flag) & 1];
	  for(i=0; i<config->pulsewidth * 4; i+=2)
	     {
	     array[i] = toieee(*g--) * 2.0;
	     array[i+1] = 0.0;
	     }
	  hann(array,config->pulsewidth*2);
	  fft(array,config->pulsewidth*2);
   
	  wp = 1.0 - 0.4 * log10((double)0x2000000);
   
	  for(i=0; i<config->pulsewidth * 2; i++)     
	     {
	     I = i*2;
	     Q = I + 1;
	     if(array[I] == 0.0 && array[Q] == 0.0) {array[i] = 0.0; continue;}
	     array[i] = 0.2 * log10(array[I]*array[I]+array[Q]*array[Q]) + wp; 
	     }

	  avedisp(array,config->pulsewidth,0.1);
	  grputs(XTITLE,YTITLE,"Gate 0 FFT        ");
	  continue;  /* go back to begining */
	  }
	    
      /* call getseries which removes the time series from dwell */
      /* it builds an array of length n from multiple smaller arrays */
      /* of size m.  Any leftover is used on the next call */
      if(getseries(array,2*n,(char *)(dwell->abp)+PARMNUM*dwell->header.gates,2*dwell->header.hits)) 
	 {
//         if(config->gate0mode)  afc(config,piraq);  

	if(lasta != a) 
	   {
	   lasta = a;
	   Lastparm = -1;
	   if(a) {gclear(); helpindex();}
	   }


	 switch(a) /* a */
	    {
	    case 0:     /* show the current product */
			products(dwell,radar,prods); 
			proddisplay(dwell,radar,prods,Dat,parameter);
			grputs(XTITLE,YTITLE,"Parameter vs Range ");
			break;
	    case 1:     /* show the current product */
			products(dwell,radar,prods); 
			histogram(dwell,radar,prods,Dat,parameter);
			grputs(XTITLE,YTITLE,"Parameter Histogram");
			break;
	    case 2:     idisplayts(dwell,array,n,1.0);
			grputs(XTITLE,YTITLE,"I vs time          ");
			break;
	    case 3:     qdisplayts(dwell,array,n,1.0);
			grputs(XTITLE,YTITLE,"Q vs time          ");
			break;
	    case 4:     iqphase(dwell,array,n,1.0);
			grputs(XTITLE,YTITLE,"I vs Q             ");
			break;
	    case 5:     displayfft(dwell,array,n,.10);
			grputs(XTITLE,YTITLE,"Power vs Frequency ");
			break;
	    }
	 }
      }

   Prf(1000,1,0);
   tmode(); 
   savev();
   }

void    helpindex(void)   
   {
   int line;

   line = 12 * THEIGHT;
   grputs(XTEXT,line+=THEIGHT,  "         COMMANDS");
   grputs(XTEXT,line+=THEIGHT*2," 1  parameters 1 - N");
   grputs(XTEXT,line+=THEIGHT,  " 2  expectation 1 - N");
   grputs(XTEXT,line+=THEIGHT,  " 3  I time series");
   grputs(XTEXT,line+=THEIGHT,  " 4  Q time series");
   grputs(XTEXT,line+=THEIGHT,  " 5  I vs Q");
   grputs(XTEXT,line+=THEIGHT,  " 6  FFT of time series");
   grputs(XTEXT,line+=THEIGHT,  " 7  Gate 0 time series");
   grputs(XTEXT,line+=THEIGHT,  " 8  FFT of gate 0");
   
   grputs(XTEXT,line+=THEIGHT*2," N  Inc H noise power");
   grputs(XTEXT,line+=THEIGHT  ," M  Dec H noise power");
   
   grputs(XTEXT,line+=THEIGHT*2," ]  Inc DAC Voltage");
   grputs(XTEXT,line+=THEIGHT,  " [  Dec DAC Voltage");
   grputs(XTEXT,line+=THEIGHT,  " =  Inc T Series Gate");
   grputs(XTEXT,line+=THEIGHT,  " -  Dec T Series Gate");
   grputs(XTEXT,line+=THEIGHT,  " .  Inc Trigger Delay");
   grputs(XTEXT,line+=THEIGHT,  " ,  Dec Trigger Delay");
   grputs(XTEXT,line+=THEIGHT,  " U  Inc Display Gain");
   grputs(XTEXT,line+=THEIGHT,  " D  Dec Display Gain");
   
   grputs(XTEXT,line+=THEIGHT*2," C  Toggle Filter");
   
   grputs(XTEXT,line+=THEIGHT*2," Q  Quit");
   }

/* display the current parameter */
double  Lastx,Lasty;
void proddisplay(DWELL *dwell, RADAR *radar, float *prods, PARAMDAT *dat, int param)
   {
   int          i;
   char         buf[80];

   if(param != Lastparm)
      {
      gclear();
      helpindex();
      drawgraphscale("Range (km)",0.0,dwell->header.gates * radar->xmit_pulsewidth * .15 / 1e-6,
		     dat[param].label,dat[param].low,dat[param].high);
      Lastparm = param; /* set lastparm to current param */
      }
   else
      uncursor(Lastx,Lasty);
   
   prods += param;
   
   Lastx = dwell->header.tsgate / (double)(dwell->header.gates - 1);
   Lasty = prods[16*dwell->header.tsgate] * dat[param].scale + dat[param].offset;

   sprintf(buf,"Cursor: %8.3f %s",prods[16*dwell->header.tsgate],dat[param].units);   
   grputs(XTEXT,10 * THEIGHT,buf);
   sprintf(buf,"  gate: %5.1f Km",0.001*C*dwell->header.rcvr_pulsewidth*dwell->header.tsgate*0.5);
   grputs(XTEXT,11 * THEIGHT,buf);
   
   for(i=0; i<dwell->header.gates; i++,prods+=16)     
      array[i] = *prods * dat[param].scale + dat[param].offset;  /* scale to +/-1 */
   
   avedisp(array,dwell->header.gates,1.0); /* domain = +/- 1 */
   
   cursor(Lastx,Lasty);
   
   }

/* display the current parameter */
void histogram(DWELL *dwell, RADAR *radar, float *prods, PARAMDAT *dat, int param)
   {
   int          i,value;
   char         buf[80];

   if(param != Lastparm)
      {
      gclear();
      helpindex();
      drawgraphscale(dat[param].label,dat[param].low,dat[param].high,
		     "Expectation",0.0,(double)dwell->header.gates);
      Lastparm = param; /* set lastparm to current param */
      }
   
   prods += param;
   
   /* do histogram assuming the range is split into 100 bins */
   
   /* clear bins */
   for(i=0; i<100; i++) array[i] = -1.0;
   
   for(i=0; i<dwell->header.gates; i++,prods+=16)     
      {
      value = 50.0 * (1.0 + *prods * dat[param].scale + dat[param].offset);
      if(value >= 0 && value < 100)
	array[value] += 2.0 / (double)dwell->header.gates;
      }

   avedisp(array,100,1.0); /* domain = +/- 1 */
   
   cursor(Lastx,Lasty);
   
   }

/* display fft (dB) of a complex array with 1.0 full scale amplitude */
void displayfft(DWELL *dwell, float array[],int size,double filter)
   {
   double       wp;
   int          i,I,Q;

   wp = 1.0 - 0.4 * log10((double)0x1000000 * dwell->header.rcvr_pulsewidth / 1.25e-7);
   
   hann(array,size);
   fft(array,size);
   
   for(i=0; i<size; i++)     
      {
      I = i * 2;
      Q = I + 1;
      if(array[I] == 0.0 && array[Q] == 0.0) {array[i] = 0.0; continue;}
      array[i] = 0.2 * log10(array[I]*array[I]+array[Q]*array[Q]) + wp; 
      }

   avedisp(array,size,filter);
   }

/* display iq in polar coordinates 1.0 full scale amplitude */
void iqphase(DWELL *dwell, float array[],int size,double filter)
   {
   double       wp;
   int          i,I,Q;

   wp = 1.0 / ((double)0x1000000 * dwell->header.rcvr_pulsewidth / 1.25e-7);
   
   for(i=0; i<size; i++)     
      {
      I = i * 2;
      Q = I + 1;
      array[I] *= TSSCALE * wp; 
      array[Q] *= TSSCALE * wp; 
      }
//pdisp(array,size);
   polardisp(array,size,filter);
   }

/* display ts of a complex array with 1.0 full scale amplitude */
void idisplayts(DWELL *dwell, float array[],int size,double filter)
   {
   double       wp;
   int          i,I,Q;

   wp = 1.0 / ((double)0x1000000 * dwell->header.rcvr_pulsewidth / 1.25e-7);
   
   for(i=0; i<size; i++)     
      {
      I = i * 2;
      Q = I + 1;
      array[i] = TSSCALE * wp * array[I]; 
      }
//disp(array,size);
   avedisp(array,size,filter);
   }

/* display ts of a complex array with 1.0 full scale amplitude */
void qdisplayts(DWELL *dwell, float array[],int size,double filter)
   {
   double       wp;
   int          i,I,Q;

   wp = 1.0 / ((double)0x1000000 * dwell->header.rcvr_pulsewidth / 1.25e-7);
   
   for(i=0; i<size; i++)     
      {
      I = i * 2;
      Q = I + 1;
      array[i] = TSSCALE * wp * array[Q]; 
      }
//disp(array,size);
   avedisp(array,size,filter);
   }

void displayp(DWELL *dwell, RADAR *radar, float *prods, double filter)
   {
   double       wp,scale;
   short        *ptr;
   int          i;

   scale =  1.0 - .02 * (radar->data_sys_sat - radar->receiver_gain);
   
   prods += 1;
    
   for(i=0; i<dwell->header.gates; i++,prods+=16)     
      array[i] = 0.02 * *prods + scale; 
   
   avedisp(array,dwell->header.gates,filter); /* range = +/- 1 */
   }