#define	TRCORR	1.0 /* Use this to disable trcorrection */
#define RSHIFT  0.0
#define	PSHIFT 	0.0
#include        <math.h>

#include        "..\include\globals.h"
      
extern double trcorrect[], trcorrtog;

/* compute 6 radar products aranged as:  */
/* velocity (m/s) */
/* power (dBm) */
/* NCP */
/* spectral width (m/s) */
/* reflectivity (dBZ) */
/* coherent reflectivity (dBZ) */

/* the products fill a floating point array with scientific parameters */
/* in regular units (i.e. m/s, dBm, dBZ, Hz ... ) */

/* note dwell->dataformat describes: 0 regular abp, 1 abpab polypulsepair, or 2 abpab dualprt */

void products(DWELL *dwell, RADAR *radar, float *prods)      
   {
   switch(dwell->header.dataformat)
      {
      case DATA_POLYPP:   polypp(dwell,radar,prods);      break;
      case DATA_DUALPP:   dualprt(dwell,radar,prods);     break;
      case DATA_POL1:     dualpol1(dwell,radar,prods);    break;
      case DATA_SIMPLEPP: simplepp(dwell,radar,prods);    break;  
      case 6:
      case 7:             newsimplepp(dwell,radar,prods); break;
      default:            simplepp(dwell,radar,prods);    break;
      }
   }

void simplepp(DWELL *dwell, RADAR *radar, float *prods)      
   {
   int  i;   
   float        *aptr,*pptr, rng;
   double       a,b,p,cp,pcorrect;
   double       logstuff,noise,velconst;
   double       dbm,widthconst,range,rconst,width,r12;
   double       sqrt();   
   printf("rconst %f, radar->rconst %f\n", rconst, radar->rconst);

   rconst = radar->rconst - 
   		20.0 * log10(radar->xmit_pulsewidth / 
			dwell->header.rcvr_pulsewidth);
   noise = (radar->noise_power > -10.0) ? 0.0 : 0.0;
   velconst = -C / (2.0 * radar->frequency * 2.0 * M_PI * dwell->header.prt);
   pcorrect = radar->data_sys_sat 
	    - 20.0 * log10(0x2000000 * dwell->header.rcvr_pulsewidth / 1.25E-7) 
	    - 10.0 * log10((double)dwell->header.hits)
	    - radar->receiver_gain
	    + PSHIFT;
   widthconst = (C / radar->frequency) / dwell->header.prt / (2.0 * sqrt(2.0) * M_PI);
   
   aptr = (float *)dwell->abp;
   pptr = prods;
   range = 0.0;

   for(i=0; i<dwell->header.gates; i++)
      {
      a = *aptr++;
      b = *aptr++;
      p = *aptr++;
      
      rng = (float)i * .0005 * C * dwell->header.rcvr_pulsewidth;
      
      if(i)     range = 20.0 * log10(rng);

      /* compute floating point, scaled, scientific products */
      *pptr++ = velconst * atan2(b,a);                 /* velocity in m/s */
      /* This one includes OPRA TR recovery correction */
      	*pptr++ = dbm = 0.004 * p + pcorrect + trcorrtog*trcorrect[i];      /* power in dBm */
      *pptr++ = (cp = sqrt(r12 = a*a+b*b))/p;                /* NCP no units */
      if((width = log(fabs((p-noise)/cp))) < 0.0) width = 0.0001;
      *pptr++ = sqrt(width) * widthconst;  /* s.w. in m/s */
      *pptr++ = dbm + rconst + range;                          /* in dBZ */
      *pptr++ = 10.0 * log10(fabs(cp)) + pcorrect + rconst + range;  /* in dBZ */
      pptr += 2;  /* 8 products total */
      }
   }

void polypp(DWELL *dwell, RADAR *radar, float *prods)      
   {
   int  i;   
   float        *aptr,*pptr;
   double       a,b,p,a2,b2,cp,pcorrect,cp2;
   double       logstuff,noise,velconst;
   double       dbm,widthconst,range,rconst,width,r12,r22;
   double       sqrt();   

   rconst = radar->rconst - 20.0 * log10(radar->xmit_pulsewidth / dwell->header.rcvr_pulsewidth);
   noise = (radar->noise_power > -10.0) ? 0.0 : 0.0;
   velconst = C / (2.0 * radar->frequency * 2.0 * M_PI * dwell->header.prt);
   pcorrect = radar->data_sys_sat
	    - 20.0 * log10(0x2000000 * dwell->header.rcvr_pulsewidth / 1.25E-7) 
	    - 10.0 * log10((double)dwell->header.hits)
	    - radar->receiver_gain;
   widthconst = 0.33333333333333333 * (C / radar->frequency) / dwell->header.prt / (2.0 * sqrt(2.0) * M_PI);
   
   aptr = (float *)dwell->abp;
   pptr = prods;
   range = 0.0;

   for(i=0; i<dwell->header.gates; i++)
      {
      a = *aptr++;
      b = *aptr++;
      p = *aptr++;
      a2 = *aptr++;
      b2 = *aptr++;
      
      r12 = a * a + b * b;
      r22 = a2 * a2 + b2 * b2;
      cp = sqrt(r12) * pow(r12/r22,0.1666666666666666);

      if(i)     range = 20.0 * log10(i * 0.0005 * C * dwell->header.rcvr_pulsewidth);

      /* compute floating point, scaled, scientific products */
      *pptr++ = velconst * atan2(-b,a);                 /* velocity in m/s */
      *pptr++ = dbm = 10.0 * log10(fabs(p)) + pcorrect;      /* power in dBm */
      *pptr++ = cp / p;                /* NCP no units */
      if((width = 0.5 * log(r12/r22)) < 0.0) width = 0.0001;
      *pptr++ = sqrt(width) * widthconst;  /* s.w. in m/s */
      *pptr++ = dbm + rconst + range;                          /* in dBZ */
      *pptr++ = 10.0 * log10(cp) + pcorrect + rconst + range;  /* in dBZ */
      pptr += 2;        /* 8 products total */
      }
   }

void dualprt(DWELL *dwell, RADAR *radar, float *prods)      
   {
   int  i;   
   float        *aptr,*pptr,*bptr;
   double       a,b,p,a2,b2,cp,pcorrect,biga,bigb;
   double       logstuff,noise,velconst;
   double       dbm,widthconst,range,rconst,width,r12;
   double       sqrt();   

   rconst = radar->rconst - 20.0 * log10(radar->xmit_pulsewidth / dwell->header.rcvr_pulsewidth);
   noise = (radar->noise_power > -10.0) ? 0.0 : 0.0;
   velconst = C / (2.0 * radar->frequency * 2.0 * M_PI * fabs(dwell->header.prt - dwell->header.prt2));
   pcorrect = radar->data_sys_sat
	    - 20.0 * log10(0x2000000 * dwell->header.rcvr_pulsewidth / 1.25E-7) 
	    - 10.0 * log10((double)dwell->header.hits)
	    - radar->receiver_gain;
   widthconst = (C / radar->frequency) / dwell->header.prt / (2.0 * sqrt(2.0) * M_PI);
   
   aptr = (float *)dwell->abp;
   bptr = aptr + 3 * dwell->header.gates;
   pptr = prods;
   range = 0.0;

   for(i=0; i<dwell->header.gates; i++)
      {
      a = *aptr++;
      b = *aptr++;
      p = *aptr++;
      a2 = *bptr++;
      b2 = *bptr++;
      p += *bptr++;
      p /= 2.0;      

      if(i)     range = 20.0 * log10(i * 0.0005 * C * dwell->header.rcvr_pulsewidth);

      /* compute floating point, scaled, scientific products */
      biga = a * a2 + b * b2;
      bigb = a2 * b - a * b2;
      
      *pptr++ = velconst * atan2(bigb,biga);                 /* velocity in m/s */
      *pptr++ = dbm = 10.0 * log10(fabs(p)) + pcorrect;      /* power in dBm */
      *pptr++ = (cp = sqrt(r12 = a*a+b*b))/p;                /* NCP no units */
      if((width = log(fabs((p-noise)/cp))) < 0.0) width = 0.0001;
      *pptr++ = sqrt(width) * widthconst;  /* s.w. in m/s */
      *pptr++ = dbm + rconst + range;                          /* in dBZ */
      *pptr++ = 10.0 * log10(fabs(cp)) + pcorrect + rconst + range;  /* in dBZ */
      pptr += 2;        /* 8 products total */
      }
   }

/* finish this later */
void dualpol1(DWELL *dwell, RADAR *radar, float *prods)      
   {
   int  i;   
   float        *aptr,*pptr,*bptr;
   double       a,b,p,a2,b2,p2,cp1,cp2,v,dp,pcorrect,biga,bigb;
   double       v1a,v1b,v2a,v2b,cp;
   double       logstuff,noise,velconst;
   double       dbm,widthconst,range,rconst,width,r12;
   double       sqrt();   

   rconst = radar->rconst - 20.0 * log10(radar->xmit_pulsewidth / dwell->header.rcvr_pulsewidth);
   noise = (radar->noise_power > -10.0) ? 0.0 : 0.0;
   velconst = C / (2.0 * radar->frequency * 2.0 * M_PI * fabs(dwell->header.prt - dwell->header.prt2));
   pcorrect = radar->data_sys_sat
	    - 20.0 * log10(0x2000000 * dwell->header.rcvr_pulsewidth / 1.25E-7) 
	    - 10.0 * log10((double)dwell->header.hits)
	    - radar->receiver_gain;
   widthconst = (C / radar->frequency) / dwell->header.prt / (2.0 * sqrt(2.0) * M_PI);
   
   aptr = (float *)dwell->abp;
   bptr = aptr + 3 * dwell->header.gates;
   pptr = prods;
   range = 0.0;

   for(i=0; i<dwell->header.gates; i++)
      {
      a = *aptr++;
      b = *aptr++;
      p = *aptr++;
      a2 = *bptr++;
      b2 = *bptr++;
      p2 = *bptr++;

      if(i)     range = 20.0 * log10(i * 0.0005 * C * dwell->header.rcvr_pulsewidth);

      /* compute floating point, scaled, scientific products */
      v1a = atan2(b,a) + radar->phaseoffset;
      if(v1a > M_PI) v1a -= 2.0 * M_PI; if(v1a < -M_PI)  v1a += 2.0 * M_PI;
      v2a = atan2(b2,a2) - radar->phaseoffset;
      if(v2a > M_PI) v2a -= 2.0 * M_PI; if(v2a < -M_PI)  v2a += 2.0 * M_PI;
      
      v1b += (v1a < 0.0) ? 2.0 * M_PI : -2.0 * M_PI;
      v2b += (v2a < 0.0) ? 2.0 * M_PI : -2.0 * M_PI;

      dp = 0.5 * (v1a - v2a);   /* or maybe v1b - v2b */
      if(dp < (160.0 * M_PI / 180.0) && /* if this doesn't work, */
	 dp > (-20.0 * M_PI / 180.0))   /* then v1b - v2b won't either */
	 {                               
	 v = 0.5 * (v1a + v2a);
	 if(v > M_PI || v < -M_PI)
	    v = 0.5 * (v1b + v2b);
	 }
      else   /* use the cross terms */   
	 {
	 dp = 0.5 * (v1a - v2b);
	 v = 0.5 * (v1a + v2b);
	 if(v > M_PI || v < -M_PI)
	    v = 0.5 * (v1b + v2a);
	 }

      cp1 = sqrt(a * a + b * b); 
      cp2 = sqrt(a2 * a2 + b2 * b2);
      cp = cp1 + cp2;

      *pptr++ = velconst * v;                 /* velocity in m/s */
      *pptr++ = dbm = 10.0 * log10(fabs(p + p2)) + pcorrect;      /* power in dBm */
      *pptr++ = cp / (p + p2);                /* NCP no units */
      if((width = log(fabs((p-noise)/cp))) < 0.0) width = 0.0001;
      *pptr++ = sqrt(width) * widthconst;  /* s.w. in m/s */
      *pptr++ = dbm + rconst + range;                          /* in dBZ */
      *pptr++ = 10.0 * log10(fabs(cp)) + pcorrect + rconst + range;  /* in dBZ */
      *pptr++ = 10.0 * log10(p2/ p);  /* Zdr */
      *pptr++ = dp * 180.0 / M_PI;                    /* PHI dp */
      }
   }

void newsimplepp(DWELL *dwell, RADAR *radar, float *pptr)      
   {
   int  i;   
   short        *aptr;
   double       rng,cp,v,p,pcorrect;
   double       velconst,dbm,widthconst,range,rconst,scale2db,scale2ln;
   double       sqrt();   

   scale2ln = 0.004 * log(10.0) / 10.0;  /* from counts to natural log */
   scale2db = 0.004 / scale2ln;         /* from natural log to 10log10() */
   velconst = C / (2.0 * radar->frequency * 2.0 * fabs(dwell->header.prt) * 32768.0);
   rconst = radar->rconst - 20.0 * log10(radar->xmit_pulsewidth / dwell->header.rcvr_pulsewidth);
   pcorrect = radar->data_sys_sat
	    - 20.0 * log10(0x10000) 
	    - radar->receiver_gain
	    + PSHIFT;
   widthconst = (C / radar->frequency) / dwell->header.prt / (2.0 * sqrt(2.0) * M_PI);
   
   aptr = dwell->abp;
   range = 0.0;
   
   for(i=0; i<dwell->header.gates; i++)
      {
      cp = *aptr++;    /* 0.004 dB / bit */
      v  = *aptr++;    /* nyquist = 65536 = +/- 32768 */
      p  = *aptr++;    /* 0.004 dB / bit */
      
      rng = (float)i * .0005 * C * dwell->header.rcvr_pulsewidth;

      if(i)     range = RSHIFT + 20.0 * log10(rng);

      /* compute floating point, scaled, scientific products */
      *pptr++ = velconst * v;                 /* velocity in m/s */
      /* This one includes OPRA TR recovery correction */
      *pptr++ = dbm = 0.004 * p + pcorrect + trcorrtog*trcorrect[i];      /* power in dBm */
      *pptr++ = exp(scale2ln*(cp - p));                /* NCP no units */
      *pptr++ = sqrt(scale2ln * fabs(p-cp)) * widthconst;  /* s.w. in m/s */
      *pptr++ = dbm + rconst + range;                          /* in dBZ */
      *pptr++ = 0.004 * cp + pcorrect + rconst + range;  /* in dBZ */
      pptr += 10;
      }
   }