/* rayexa -- 
  gcc -c -O rayexa.c  && gcc -O -o rayexa rayexa.o -lm -lcd
 */

#include <stdio.h>
#include <math.h>
#include "cdfhdr.h"

#define GATES 512        /*Number of gates in ray*/
#define MAXCHAR 100      /*Max characters on the calling line*/
#define D2R 0.01745329278  /* pi/180 conversion factor (degrees to radians)*/

/* assumed names of radial velocity and reflectivity fields */
#define VR "vr"
#define DBZ "dbz"

/* assumed names defining positioning */
#define RANGE "range"   /* range */
#define ALT "alt"	/* altitude */
#define EL "el"		/* elevation */

/* candis stuff */
static char hbuff[HBMAX][LINE];
static char cline[LINE];
static float bad,badlim;
static float *sbuff,*vbuff;
static long nsbuff,nvbuff;
static struct field *fp;

/* names for pointers to input fields */
static float *range;		/* --> RANGE (index field) */
static float *dbz;		/* --> DBZ */
static float *el;		/* --> EL */
static float *alt;		/* --> ALT */
static float in_vr[512];        /* not required radial velocity --> VR */

/* functions */
void fun_make_header();
void fun_alloc();
void fun_make_s_slice();
void fun_stop(char *s);
void fun_usage(char *s);
void fun_cgeta(int argc,char **argv);
void fun_dbzthresh();
void fun_minrange();
void fun_seaclutter();
void fun_cut_surf_layer();


/* line arguments -- treated as global for simplicity*/
static char argu[MAXCHAR];
static float dbz_thresh;   /* reflectivity threshold (dbz) */
static float min_range;    /* minimum range for which we accept gates */
static float cl_h;         /*clutter layer height (km) */
static int flag_c = 0;     /* Flags for different arguments options */
static int flag_r = 0;
static int flag_t = 0;
static int flag_s = 0;

main(int argc,char **argv)
{
  float refl;           /* reflectivity factor */
  
  /* check and get the command line Parameters */
  fun_cgeta(argc,argv);
  
  /* Candis header */
  fun_make_header();
  
  /* allocate memory for static & variable fields */
  fun_alloc();
  
  /* Candis static slice */
  fun_make_s_slice();
  
  /* loop through input variable slices */
  while (getslice(stdin,nvbuff,vbuff) != EOF) {

      if (flag_t) fun_dbzthresh();
      if (flag_r) fun_minrange();
      if (flag_s) fun_seaclutter();
      if (flag_c) fun_cut_surf_layer(); /*has to be done after*/


      
    putslice(stdout,nvbuff,vbuff);
  } /* read_slice_while */
  
    exit(0);
} /* EOmain */

/*--------------------------------------------------------------------*/
/*refl = 10.*log10(dbz[iloop]);*/
void fun_dbzthresh(){
float aux;
int loop;
      for (loop = 0; loop < GATES; loop++) { /* loop on range */
	  aux = dbz[loop];
          if (aux < dbz_thresh)
             dbz[loop] = bad;
      } /* range_loop */
}
/*--------------------------------------------------------------------*/
void fun_minrange(){
float aux,faux;
int iaux,loop;

      /*find last gate to be supressed*/
      for (loop = 0; loop < GATES; loop++) { /* loop on range */
	  if (range[loop] > min_range) {
	      iaux = loop;
	      break;
          }
      }

      /*supress dbz-values from first to iaux gate.*/
      for (loop = 0; loop < iaux; loop++) {
             dbz[loop] = bad;
      } /* range_loop */
}
/*--------------------------------------------------------------------*/
void fun_cut_surf_layer(){
/* Cleans up below the flight level (i.e. downward pointing rays.) 
   It has to be done after the sea clutter is removed,
   from reflectivity data collected above flight level.
*/
int loop;
int n_c = -1;
   for (loop = 0; loop < GATES; loop++) {
       if ( (*el < 0) && 
            (*alt - (range[loop] * sin((fabs(*el) + 1.0) * D2R )) < cl_h) ) {  
	  n_c = loop - 1; 
	  break; 
       }
   }
   if (n_c > 0 ) { /*clean up everything beyound n_c gate.*/
      for (loop = n_c; loop < GATES; loop++) {
          dbz[loop] = bad;
      }
   }
}
/*--------------------------------------------------------------------*/
   /* Attempt to despecle a ray. Not done yet!!
   This is relevant only if unfolding velocities is required.
   I will implement this in cdfrayvr where, tsegments of good
   data is already done. Hence, the only thing required is to
   get rid of good segments of length one.
   */
/*--------------------------------------------------------------------*/
void fun_seaclutter(){ 
/* elimination of specular reflection by discarding gates
   whose range value would put them beyoud the surface regardles where
   the ray is pointing to. This is an attempt of implementing the thechnique
   described by Testud ... complete reference*/

#define N_t 5
#define N_d 10
#define DBZ_INIT -100  /*initial value in the search of max in dbz*/


int loop,n1_d,n2_d,n_max;
float rs,aux,faux,dbz_max;

   /* get nominal ray-distance to surf based on fabs(el)*/
   rs = *alt/fabs(sin((*el) * D2R)); 

   if ( rs > range[N_d] ) {

      /* get interval centered in rs. This interval might be contaminated 
         by sea clutter regardless to where the ray is pointing*/
      for (loop = N_d ; loop < GATES; loop++ ) {
          if ( range[loop] > rs) {
	     n1_d = loop - N_d; 
	     n2_d = loop + N_d; if ( n2_d > (GATES - 1)) n2_d = GATES - 1 ;
	     break;
          }
      }
      /* find max in dbz within the interval centered in rs*/
      dbz_max =  DBZ_INIT;
      n_max   = -1;
      for (loop = n1_d; loop <= n2_d; loop++ ) {
	  aux = dbz[loop];
	  faux = fabs(aux);
	  if ( faux < badlim && aux > dbz_max) {
	        dbz_max = aux;
                n_max   = loop;
	        n1_d = n_max - N_t;
	        n2_d = n_max + N_t;
          }
      }
      if (n_max > 0) { /* discard everything between n1_d and n2_d
	                regardless  where the ray is pointing to*/
         for (loop = n1_d; loop <= n2_d ; loop++ )  { dbz[loop] = bad;}
         if ( *el < 0 ){ /*Further discard everything beyound gate n2_d*/
            for (loop = n2_d; loop < GATES; loop++ ) { dbz[loop] = bad;} 
         }
     } 
   }
}
/*--------------------------------------------------------------------*/
void fun_cgeta(int argc,char **argv) {
  
  /* fun_cgeta stuff */
  char *STDMESSAGE=" "; 
  int i,argtc;		/* total */
  int fixarg = 1;	/* fixed arguments */
  int maxarg = 7;	/* maximum number of arguments */
  
  /* getopt stuff */
  extern char *optarg;
  /* extern int optind = 1; */
  int c;
  int j,k,l;
  
  if (argc > maxarg+1) {fprintf(stderr,"argc = %i\n",argc); exit(1);}
  
  /* Process Arguments */
  while ((c = getopt(argc,argv,"c:r:t:sh")) != -1) {
    switch (c) {
    case 'c' : flag_c = 1.0;
      if (sscanf(optarg," %f",&cl_h) != 1)
	fun_usage(": incorrect format for clutter_layer_height");
      argtc+=2;
      break; 
    case 'r' : flag_r = 1.0;
      if (sscanf(optarg," %f",&min_range) != 1)
	fun_usage(": incorrect format for minimum_range");
      argtc+=2;
      break;
    case 't' : flag_t = 1.0;
      if (sscanf(optarg," %f",&dbz_thresh) != 1)
	fun_usage(": incorrect format for dbz_threshold");
      argtc+=2;
      break;
    case 's' : flag_s = 1.0; /*specular reflection*/
      argtc+=1;
      break;
    case 'h' : fun_usage(STDMESSAGE);
      break;
    default  : fun_usage(STDMESSAGE);
      break;
    } 
  }

  /*get string with all called arguments in it*/
  l = 0;
  for ( k=0 ; k < argc; k++) {
      j=0;
      while ( (c = argv[k][j++]) != '\0' ) { 
         if (l < MAXCHAR) argu[l++] = c;
	 else fun_stop("Not enough room to hold data");
      }
      argu[l++] = ' ';
  } argu[l] = '\n';

  
} /* EOfun_cgeta */

/*--------------------------------------------------------------------*/
void fun_usage(char *s)
{
  char *USAGE = "Usage:cdfraydbz  [-c h(km)][-r min_range][-t dbz_thresh][-s][-h]\n";
  
  if (s != " ") fprintf(stderr,"%s\n",s);
  fprintf(stderr,"%s",USAGE);
  exit(1);
} /* EOfun_usage */

/*--------------------------------------------------------------------*/
void fun_stop(char *s)
{
  fprintf(stderr,"@fun_%s",s);
  exit(1);
}

/*--------------------------------------------------------------------*/
void fun_make_header() {
  
  /* get header of input file */
  gethdr(stdin,hbuff,HBMAX);
  if (getfmt(hbuff,HBMAX) != 'f') {
    fun_stop("make_header:format expected");
  }
  /*get bad and badlim */
  if (seekpar(hbuff,HBMAX,"bad",cline) == OK) bad = atof(cline);
  else {fprintf(stderr,"Could't find bad value\n"); exit(1);}
  if (seekpar(hbuff,HBMAX,"badlim",cline) == OK) badlim = atof(cline);
  else {fprintf(stderr,"Could't find badlim value\n"); exit(1);}
  
  /* make header of output file */
  /* comments */
  /*sprintf(cline,":  %s\n",argu); */
  addcline(hbuff,HBMAX,argu);
  
  /* write output header */
  puthdr(stdout,hbuff,HBMAX);

} /* EOfun_make_header */

/*--------------------------------------------------------------------*/
void fun_alloc() {
  
  /* allocate buffers */
  nsbuff  = elemcnt(hbuff,HBMAX,'s');
  nvbuff  = elemcnt(hbuff,HBMAX,'v');
  sbuff   = getbuff(nsbuff);
  vbuff   = getbuff(nvbuff);
  
  /* get pointers for index fields */
  range = getptr2(hbuff,HBMAX,sbuff,'s','d',RANGE,&fp);
  if (fp->dim != 1 || strcmp(fp->fname,fp->dname1) != 0) {
    fprintf(stderr,"@alloc: %s is not an index field\n",fp->fname);
    exit(1);
  }
  if ( fp->dsize1 != GATES ) { 
      fprintf(stderr,"#of gates is not 512!! \n"); exit(1);
  }

  /* get pointers for fields to be cartesianized */
  dbz    = getptr(hbuff,HBMAX,vbuff,'v','d',DBZ);
  
  /* get pointer  elevation */
  el  = getptr(hbuff,HBMAX,vbuff,'v','d',EL);
  alt = getptr(hbuff,HBMAX,vbuff,'v','d',ALT);
  
} /* EOfun_alloc */

/*--------------------------------------------------------------------*/
void fun_make_s_slice() {
  
  /* get input static slice */
  getslice(stdin,nsbuff,sbuff);
  
  /* write output static slice */
  putslice(stdout,nsbuff,sbuff);
  
} /* EOfun_make_s_slice */

/*--------------------------------------------------------------------*/