/* 
Name  :radfills
Author: Carlos Lopez-Carrillo
Date  : December-19-1996.
Last Modification before New version: --
Name to Archive : radfills1.c
Status: Exp
Description: the idea is to run over the grid and substitute every
value on the grid of the mentioned fields by a weighted average. 
The points to make the interpolation must be in the same z-plane 
and inside or over the radius of influence. Also there must be a 
minimun number of points into that radius to make a sensible average. 
the value for the radius of influenc(radius),the minimun value of good
points into that radius(minp) and the list of fields to be interpolated  
must be provided as an input parameters.
Default values are:
                  radius = sqrt(dx*dx + dy*dy) 
this will take into account the value of the point plus the eignt 
surrounding it -- provided they lie inside the actual grid. 
                  minp  = 1;
this allow a good point to survive even if there no other good points 
arround. 

Note: As we go on the grid, we calculate the average values for every
mentioned field based on the originals values of the input fields, the
average value is then saved on a parallel array. When we are done with
every point on the grid, we made another run over the grid to update
the fields before they are written to the output.

MOD-Feb-14-1997:This is to add two fields to the output: 
The Detreined-Volume-Flux Density (dvfd) and the Detrained-Mass-Flux Density
The reason is to get the detreined fluxes corrected by a
column-by-column divergence correction on them. the idea goes as folows:

Let D be the Uncorrected divergence from the smoothed velocity fields u,v. 
and Let Dc be the Corrected divergence on each column such that Dc = D + C
where C is the correction term. To find the correction term, we apply the 
condition that the Integral of Dc times the density, R, from zero to Hmax 
must be ZERO.  Then (in Latex notation) we would write,

       C = - \Int_{0}^{Hmax} R*D dz  /  \Int_{0}^{Hmax} R dz

Note: D = \part{u} / \part{x} + \part{v} / \part{y}

Note: that the assumed names for the horizontal velocities are
"u" and "v" this fields must exist in the input file. This program will
check how many of this velocities are on the command line, if this
number is 2 then the filed for the detreined volume flx density will
be created. If the names in the input file are different from those ,
the smoothing is still going on if the actual names are provided in the
command line. However the field (dvfd) will not be created. The names 
assumed are easily changed on the defiene statement below. I dont like
this but for the time being is gonig to be like that.
End-of-MOD-Feb-14-1997

MOD-Jul-14-1997: I find a trivial bug. The bug was that the statements
inside an if were note properly surronded by braces. This was not a problem 
until I decide to smooth fiels which were not the two components of the horizontal
velocity so the fields "dvfd" and "dmfd" were not added. The if
causing the troble was that before requesting memory for "dvfd" and
"dmfd". Since those field were not added the request of memory driven
by the if was the cause of the troble. Bug fixed!
End-of-MOD-Jul-14-1997

Compilation line:
gcc -c -O radfills.c && gcc -o radfills radfills.o -lm -lcdf
*/

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

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

#define I3(ix,iy,iz) (ix)*ny*nz + (iy)*nz + (iz)  /* access to 3-d arrays */
#define I2(ix,iy) (ix)*ny + (iy)                  /* access to 2-d arrays */

#define MAXFIELDS 25

/* names used for the grid steps*/
#define DX "dx"
#define DY "dy"
#define DZ "dz"

/*names asummed for the velocity fields*/
#define Vu "u"
#define Vv "v"

/* candis stuff */
static char hbin[HBMAX][LINE];
static char cline[LINE],pline[LINE];
static float bad,badlim;
static float *sbout,*vbout;
static long nsbin,nvbin,nsbout,nvbout;
static struct field *fp;

/*global variables */
static int nx,ny,nz;                     /* array limits */
static long num;                         /* total ponints on the grid */
static float dx,dy,dz;                   /* grid steps*/
static float *fields[MAXFIELDS];         /* fields to be interpolated */
static float *ofields[MAXFIELDS];        /* holder for ouput values of the fields*/
static float *dvfd;                      /* detreined volume flux density*/
static float *dmfd;                      /* detrained mass flux density*/
static float *rho;                       /* mass density*/
static int nfields;                      /* # fields to be interpolated */
static float radius;                     /* radius of influence */
static int minp;                         /* min num of good points inside influence circule*/
static int nhvel=0;                      /* Number of horizontal velocities to be smmothed*/
static int mvx;                          /* marck to "u" on ofields*/
static int mvy;                          /* marck to "v" on ofields*/
static float column_mass;
static float *fpaux2;

/* functions headers */
void fun_usage(char *s);
void fun_stop(char *s);
void fun_make_header(int argc,char **argv);
void fun_alloc(int argc,char **argv);
void fun_make_s_slice();
void fun_cgeta(int argc,char **argv);
void fun_initia(float *radius,int *minp,long *num,float **ofields);
void fun_grid_winterp(float **ofields);
void fun_update_fields(float **fields);
void fun_ini_dvmfd();
void fun_get_dvfd();
float fun_deriv(float *f, long xm,long xc,long xp, float dx,int N);
float fun_vertint(float *f ,float dl, int N);
float fun_vertint2(float *f ,float dl, int N);
float fun_column_mass(float *f );
/*------------------------------------------------------------------*/

main(int argc,char **argv){
int i,j,k;              /* looping variable*/
long cnt = 0;           /* v-slice counter */
float faux;

    /* check and get the command line Parameters */
    fun_cgeta(argc,argv);

    /* read inheader and make outheader */
    fun_make_header(argc,argv);

    /* allocate memory (get buffer) for static & variable  fields */
    fun_alloc(argc,argv);

    /* get s slice and make the newone*/
    fun_make_s_slice();

    /* initialization */
    fun_initia(&radius,&minp,&num,ofields);

    while(getslice(stdin,nvbin,vbout) != EOF) { /* read_slice_while */

         fun_grid_winterp(ofields);

         fun_update_fields(fields);

         if (nhvel == 2) fun_get_dvfd();
         
         putslice(stdout,nvbout,vbout);
         cnt++;

    } /* read_slice_while */

    exit(0);

}/*Eofmain*/
/*-------------------------------------------------------------------- 
       F U N C T I O N S    I M P L E M E N T A T I O N S
--------------------------------------------------------------------*/

void fun_usage(char *s){
char *USAGE = "radfills [-r radius] [-m p] field1 field2 ... < stdin > stdout\n";

    if (s != " ") fprintf(stderr,"radfills@cegeta:%s\n",s);
    fprintf(stderr,"%s",USAGE);
    exit(1);
} /*End-of-fun_usage*/

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

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

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

void fun_cgeta(int argc,char **argv) {
char c;
char *STDMESSAGE=" ";
int cop = 0; /*counter options */

/* getopt stuff */
extern char *optarg;
extern int optind;


  if (argc < 2) fun_usage(STDMESSAGE);

  radius = -1.0;
  minp = -1;
  
  nfields = 0;

  while ((c = getopt(argc,argv,"hr:m:")) != -1) {
   switch (c) {
       case 'r' : if(sscanf(optarg,"%f",&radius) != 1)
                    fun_usage("incorrect format for -r option");
                 cop +=2;
       break;
       case 'm' : if(sscanf(optarg,"%i",&minp) != 1)
                    fun_usage("incorrect format for -m option");
                 cop +=2;
       break;
       case 'h' : /* get help */
                 fun_usage(STDMESSAGE);
       break;
       default  : if ( (argc - 1 - cop) < MAXFIELDS )
                    fun_usage("too many fields");
       break;
    }
  }/*while-getopt*/

  nfields = argc - 1 - cop; /* Number of fields to be interpolated */

}/*End-of-cegeta*/


/*-----------------------------------------------------------------------
        C A N D I S-- F O R M A T    M O D U L E

Name  :cdfformat -- modified for radfills 
Author:Carlos Lopez-Carrillo
Date  :December-19-1996
Status: Exp
Description: take care of the candis format.i.e. get the input header
build the new one, and handel the memory allocation of the input aoutput fields.

This module contains: 

void fun_make_header(int argc,char **argv);
void fun_alloc(int argc,char **argv);
void fun_make_s_slice();

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

void fun_make_header(int argc,char **argv) {
int i; 
char auxstr[LINE];
char delta[LINE];       /* buffer to hold grid step sizes */
char comment[HBMAX][LINE];
int dime[MAXFIELDS];

    /* get header of input file */
    gethdr(stdin,hbin,HBMAX);
    if (getfmt(hbin,HBMAX) != 'f') {
	fun_stop("make_header-E01:format expected");
    }

    /* get the number of elements in the input file */
    nsbin  = elemcnt(hbin,HBMAX,'s');
    nvbin  = elemcnt(hbin,HBMAX,'v');

    /* get badlim parameter*/
    if (seekpar(hbin,HBMAX,"badlim",pline) == OK) badlim = atof(pline);
    else fun_stop("badlim does not appear in input file");

    /* get bad parameter*/
    if (seekpar(hbin,HBMAX,"bad",pline) == OK) bad = atof(pline);
    else fun_stop("bad does not appear in input file");

   /* check for existence of the input field */
   for(i = 1 ; i <= nfields; i++) { /*loop existence*/
      if((fp = seekfld(hbin,HBMAX,'v',argv[argc-i])) == NULL) {
        fprintf(stderr,"radfills: there is no %s field\n",argv[argc-i]);
        exit(1);
      }
      dime[i] = fp->dim;
      if ( (dime[i] != 2) && (dime[i] != 3) )
         fun_stop("make_header:dimension is nethier 2 nor 3");

      /*check how many horizontal velocities are 
        to be smoohted and mark them*/
       if (strcmp(fp->fname,Vu) == 0) {
           mvx=nfields-i; 
          nhvel++;
       }
       if (strcmp(fp->fname,Vv) == 0) {
          mvy=nfields-i;
          nhvel++;
       }
   } /*loop existence*/

   /* check dimension compatibility among the fields */
   for(i = 1 ; i < nfields; i++) /* dimension compatibility loop*/
      if(dime[i] != dime[i+1]) 
        fun_stop("make_header:fail in dimension compatibility");


   /* since fp shold hold the info from the last filed on the 
      loop existence  and we already pass the dimension compa- 
      tibility loop the we can do the following -- ok! */
   nx=ny=nz=1;
   if ( dime[nfields] >=1 ) nx = fp->dsize1; 
   if ( dime[nfields] >=2 ) ny = fp->dsize2;
   if ( dime[nfields] >=3 ) nz = fp->dsize3;

   /* ok! */

    /* get deltas */
    if(seekpar(hbin,HBMAX,DX,delta) == FAIL) 
      fun_stop("make_header:can not find DX");
    dx = atof(delta);

    if(seekpar(hbin,HBMAX,DY,delta) == FAIL) 
      fun_stop("make_header:can not find DY");
    dy = atof(delta);

    if(dime[nfields] >=3) {
       if(seekpar(hbin,HBMAX,DZ,delta) == FAIL) 
       fun_stop("make_header:can not find DZ"); 
       dz = atof(delta);
    }


   /* make header of output file */
   /*1) comments */

   strcpy(auxstr,"radfills ");
   for(i = 1; i <= (argc - nfields - 1); i++) {
      strcat(auxstr,argv[i]);
      strcat(auxstr," ");
   }
   for(i = (argc - nfields); i <= (argc-1); i++){
     strcat(auxstr,argv[i]);
     strcat(auxstr," ");
   }

   sprintf(cline,"%s\n",auxstr);
   addcline(hbin,HBMAX,cline);

   /* add new fields  */
   /* if both horizontal velocities are to be smoothed then crate a
      field for the dtreined volume flux from the smoothed velocity 
      fields with a column divergence correction on it.*/
   if (nhvel == 2){
      addfld(hbin,HBMAX,'v',"dvfd",10.,0.,'s',3,"x",nx,"y",ny,"z",nz);
      addfld(hbin,HBMAX,'v',"dmfd",10.,0.,'s',3,"x",nx,"y",ny,"z",nz);
   }

    /* write output header */
    puthdr(stdout,hbin,HBMAX);

} /*EOfun_make_header*/
/*--------------------------------------------------------------------*/

void fun_alloc(int argc,char **argv){
int ic;      /*loop var*/
int fc = 0;  /* field counter*/


    
    /* allocate buffers */
    nsbout = elemcnt(hbin,HBMAX,'s');
    nvbout = elemcnt(hbin,HBMAX,'v');
    sbout  = getbuff(nsbout);
    vbout  = getbuff(nvbout);

    /* get pointers for fields to be emploied */
    for(ic = (argc - nfields); ic <= (argc-1); ic++){
       fields[fc] = getptr(hbin,HBMAX,vbout,'v','d',argv[ic]);
       fc++;
    }
    /* get pointers for new fields*/
    if (nhvel == 2) { /* MOD-Jul-14-1997 */
       dvfd = getptr(hbin,HBMAX,vbout,'v','d',"dvfd");
       dmfd = getptr(hbin,HBMAX,vbout,'v','d',"dmfd");
       rho =  getptr(hbin,HBMAX,sbout,'s','d',"rho");
    }


}/*EOfun_alloc*/
/*--------------------------------------------------------------------*/

void fun_make_s_slice(){
int i;

    /* get input static slice */
    getslice(stdin,nsbin,sbout);

    /* write output static slice */
    putslice(stdout,nsbout,sbout);

}/*EOfun_make_s_slice*/

/*--------------------------------------------------------------------*/
void fun_initia(float *radius,int *minp,long *num,float **ofields){

/*Note: it is not necessary to initialise the values of the ofields
since fun_grid_winterp take care of that in doing so it avoid to go
an extra time over the grid*/

int fc;
float *auxp;

    if (*radius < 0.) *radius = sqrt(dx*dx + dy*dy);
    if (*minp < 0 )   *minp = 1;

    *num = nx*ny*nz;
    for(fc=0; fc < nfields ; fc++){ /*get buffers for holders*/
       auxp = getbuff( (*num) );
       ofields[fc]=auxp;
    }

}/*End-of-initia*/
/*--------------------------------------------------------------------*/
void fun_grid_winterp(float **ofields){

/* -- the default grid is defined as a three dimensional with 
first index in the x direction second in the y-direction and  
third in the z-direction. 

for a 3-d field, the neiborhood is taken to be all point into the radius
of the sqare defined by the four closet neighbours to the point on the z
plane.  otherwise it is defined by the input parameter however the
neighbours are still on the z-plane.

for a 2-d field, we search the same neiborhood and in the same way
as in the 3-d case.
Note: the input fields are not actually modified by this routine, 
instead parallels arrays to the original ones are filled with the
results of the interpolation.  -- */

/* globals: 
float radius 
int avepoints
float dx,dy
float **fields
float badlim,minp 
int nx,ny,nz;
long num
*/

int kc,jc,ic;
int mfield;
int gcx,gcy;
long n;              /*position of the field value we are lookin at */
long p;              /*position of the field values to make w-average*/
float dfp;           /*distance to the point*/
float term;          /*value of the field*/
float weight;        /*geometric weight for the value of the field*/
float wsum;          /*weigthed sum of the terms*/
float deno;          /*sum of the weights to make the normalization*/
int avepoints;       /*actual number of poit to make average*/
int nrx,nry;         /*define a rectangular domain containing the radius*/


  for(kc = 0; kc < nz ; kc++){ /*kc-grid loop*/
     for(jc = 0; jc < ny ; jc++){ /*jc-grid loop*/
        for(ic = 0; ic < nx ; ic++){ /*i-grid loop*/
            
           n   = I3(ic,jc,kc); /*point to be replaced if criteria is satisfied */
               
           for(mfield = 0 ; mfield  < nfields; mfield++){ /*fields loop */
               
              wsum = deno = 0.0;
              avepoints = 0;
              nrx = floor(radius/dx);
              nry = floor(radius/dy);
              for(gcx = -nrx; gcx <= nrx; gcx++){ /*winter-grid loopx*/
                 for(gcy = -nry; gcy <= nry; gcy++){ /*winter-grid loopy*/
                    
                    dfp = sqrt((gcx*dx)*(gcx*dx) + (gcy*dy)*(gcy*dy));
                    if( dfp <= radius) { /*check-on-range if*/
                       
                      p = I3(ic+gcx,jc+gcy,kc);
                      
                      if( p > 0 && p < num){ /* limit-to-actual-grid if*/
                        term = fields[mfield][p];  
                        if(fabs(term) < badlim) { /* check-for-good-points if*/
                          weight  = (radius - dfp)/(radius + dfp);
                          wsum  += term*weight;
                          deno  +=weight;
                          avepoints++; 
                        } /* check-for-good-points if*/
                      } /* limit-to-actual-grid if*/ 
                    } /*check-on-range if*/ 
                 } /*winter-grid loopy*/
              } /*winter-grid loopx*/
                 
              /* normalize field & criteria to keep bad data */ 
              if( (deno < 1.) || (avepoints <  minp) ) 
                ofields[mfield][n] = bad;    /* this avoid an initialization loop */   
              else 
                ofields[mfield][n] = wsum/deno;
               
           } /*fields loop */
           
        } /*ic-grid loop*/
     } /*jc-grid loop*/
  } /*kc-grid loop*/

}/*End-of-grid_winterp*/

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

void fun_update_fields(float **fields){
long lc;
int fc;
  for(lc = 0; lc < num; lc++){
     for(fc = 0; fc < nfields; fc++){
        fields[fc][lc] = ofields[fc][lc];
     }
  }   
}/*End-of-update_fields*/

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

void fun_ini_dvmfd() {

/* get uncorrected smoothed fluxes: this function retrive the 
detreined volume density as the horizontal divergence (dvfd) 
from the smoothed velocity fields without any column divergence 
correction and also set the detreined mass flux density (dmfd). */

int i,j,k;
long p,pw,pe,ps,pn;
float dudx,dvdy;
float ro;

    for( k=0; k < nz; k++ ) {
       ro = rho[k];
       for(j=0; j < ny; j++ ) {
          for( i=0; i < nx; i++ ) {

             p=I3(i,j,k);
             pw=I3(i-1,j,k);
             pe=I3(i+1,j,k);
             pn=I3(i,j+1,k);
             ps=I3(i,j-1,k);

             dudx = fun_deriv(fields[mvx],pw,p,pe,dx,nx-1);
             dvdy = fun_deriv(fields[mvy],ps,p,pn,dy,ny-1);

             if ( (fabs(dudx) < badlim) && (fabs(dvdy) < badlim) ){
               dvfd[p] = dudx + dvdy;
               dmfd[p] = dvfd[p] * ro;
             }
             else 
               dvfd[p] = dmfd[p] = bad; 

          }/*i-loop*/
       }/*j-loop*/
    }/*k-loop*/

}/*E-of-fun_get_unsmfl*/

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

float fun_column_mass(float *f ){

int k;
float mass;
float fpaux[nz];

    for(k=0; k < nz; k++){
       if( fabs(f[k]) < badlim)
         fpaux[k] = rho[k];
       else
         fpaux[k] = bad;
    }

    mass = fun_vertint2(fpaux ,dz,nz-1);

    return(mass);

}/*E-of-fun_vertint*/

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

void fun_get_dvfd() {

/*this function generate the corrected version of the detrained volume
flux*/

int i,j,k;
long p;
float cdc;

    /* set the initial value for the detreined volume and mass flux densities.
     -- this is an initialization step and must be done for the whole grid.*/

    fun_ini_dvmfd();

    /*get the corrected detrained-volume-flux density*/
    for( i=0; i < nx; i++ ) {
       for(j=0; j < ny; j++ ) {

          cdc = fun_vertint2(&dmfd[I3(i,j,0)],dz,nz-1);
          column_mass = fun_column_mass(&dmfd[I3(i,j,0)]);

          if (fabs(cdc) < badlim && fabs(column_mass) < badlim ) {
              for( k=0; k < nz; k++ ) {
                 p=I3(i,j,k);
                 if (fabs(dvfd[p]) < badlim ) {
                    dvfd[p] = dvfd[p] - (cdc/column_mass);
                    dmfd[p] = dvfd[p]*rho[k];
                 }
                 else { 
                    dvfd[p] = dmfd[p] = bad; 
                 }
              }/*k-loop*/
         }/*if*/

       }/*j-loop*/
    }/*i-loop*/


}/*E-of-fun_get_dvfd*/

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


float fun_deriv(float *f, long xm,long xc,long xp, float dx,int N){

float deriv;

    /* try centred scheme out of the borders*/
    if ( (xm >= 0) && (xp <= N) && (fabs(f[xp]) < badlim) && (fabs(f[xm]) < badlim) )
          deriv = (f[xp] - f[xm])/(2.0*dx);

    else /* try backward  out of the lower border*/

       if ( (xm >= 0) && (fabs(f[xm]) < badlim) && (fabs(f[xc]) < badlim) )
          deriv = (f[xc] - f[xm])/dx;

       else /* try forward out of the upper border*/ 

           if ( (xp <= N) && (fabs(f[xc]) < badlim) && (fabs(f[xp]) < badlim) )
              deriv = (f[xp] - f[xc])/dx;

           else /* put bad value for part_dx */

             deriv = bad;

    return(deriv); 
}/*E-of-fun_deriv*/

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

float fun_vertint(float *f ,float dl, int N){
/* Trapazoidal Rule */

int i;
float s=0;

    if( fabs(f[0]) < badlim) s = s + 0.5*f[0];
    if( fabs(f[N]) < badlim) s = s + 0.5*f[N]; 
    for(i=1; i < N; i++){
       if( fabs(f[i]) < badlim) s = s + f[i];
    } 
    s=s*dl;

    return(s);
}/*E-of-fun_vertint*/

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

float fun_vertint2(float *f ,float dl, int N){
/* Trapazoidal Rule */

int i;
float s=0;
int count=0;

    for(i=0; i <= N; i++){
       if( fabs(f[i]) < badlim) {
         s = s + f[i];
         count=count+1;
       }
    }
    if ( count > 0 ) {
      s=s*dl;
      if (fabs(s) > badlim)
         fprintf(stderr," WARNING IN INTEGRAL	%f\n");
    }
    else
      s=bad;

    return(s);
}/*E-of-fun_vertint*/

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