/***************************************************************

      deriv.c

         Routines for calculating differences and derivatives.

         Header: vector.h
         Library: vector_x.lib

*/
#include "common.h"
#include "ioserv.h"  /* check_error() */
#include "vector.h"

/***************************************************************

   dif1 and avg_adjacent1

      Functions for taking first differences.  The input and
      output arrays can be the same.  The difference x(i+1)-x(i)
      is placed in element i of the output array by dif1.
      avg_adjacent1 similarly puts the average of the i+1th and
      the ith elements of the source into the ith element of
      the destination.

*/

void dif1(float array[], float dif[], int n)
{
   int i;
   float *x0, *x1, *d;

   x0 = array;
   x1 = array + 1;
   d = dif;
   for (i=0; i<n-1; i++, x0++, x1++, d++)
   {
      if (good_float(*x0) && good_float(*x1))
         *d = *x1 - *x0;
      else  *d = BADFLOAT;
   }
   *d = BADFLOAT;    /* the last element */
   return;
}                       /* dif1 */

void avg_adjacent1(float array[], float array_aa[], int n)
{
   int i;
   float *x0, *x1, *a;

   x0 = array;
   x1 = array + 1;
   a = array_aa;
   for (i=0; i<n-1; i++, x0++, x1++, a++)
   {
      if (good_float(*x0) && good_float(*x1))
         *a = (*x1 + *x0) * 0.5;
      else  *a = BADFLOAT;
   }
   *a = BADFLOAT;    /* the last element */
   return;
}                       /* avg_adjacent */

void dif2(float array[], float dif[], int n)
{
   dif1(array, dif, n);
   dif1(dif, dif, n);
   return;
}                       /* dif2 */

void avg_adjacent2(float array[], float array_aa[], int n)
{
   avg_adjacent1(array, array_aa, n);
   avg_adjacent1(array_aa, array_aa, n);
   return;
}                       /* avg_adjacent2 */

/***************************************************************

   deriv1, deriv2

      These functions calculate the first and second
      derivatives, and again, the source and destination
      arrays can be the same or different.  The functions
      calculate dy/dx, and modify the x array accordingly.

*/

void deriv1(float x[], float y[], float xa[], float yd[], int n)
{
   int i;
   float *dx, *pdx, *pdy;

   dx = (float *)calloc(n, sizeof(float));
   check_error((dx==NULL), "Out of memory in deriv.");

   dif1(x, dx, n);
   dif1(y, yd, n);
   for (i=0, pdx=dx, pdy=yd; i<(n-1); i++, pdx++, pdy++)
               /* bug fix: n to n-1, and extra check */
   {
      if (good_float(*pdy) && good_float(*pdx) && (*pdx) != 0.0)
         *pdy = (*pdy)/(*pdx);
      else
         *pdy = BADFLOAT;
   }
   *pdy = BADFLOAT;  /* the last element */
   avg_adjacent1(x, xa, n);

   free(dx);
   return;
}                       /* deriv1 */

void deriv2(float x[], float y[], float xa[], float yd[], int n)
{
   deriv1(x, y, xa, yd, n);
   deriv1(xa, yd, xa, yd, n);
   return;
}                       /* deriv2 */
#if 0
#include "prtarray.h"  /* print_array_*() */

int main(void)
{
   float data[20];
   float depth[20];
   make_uniform_grid(1, 20, 20, depth);
   print_array_float(stdout, depth, 20, 4, "%10.f ");
   make_sinusoid(data, 20, 1, 3, 0);
   print_array_float(stdout, data, 20, 4, " %10f");
   dif1(data, data, 20);
   print_array_float(stdout, data, 20, 4, " %10f");
   dif2(data, data, 20);
   print_array_float(stdout, data, 20, 4, " %10f");
   avg_adjacent1(depth, depth, 20);
   print_array_float(stdout, depth, 20, 4, "%10.f ");
   avg_adjacent2(depth, depth, 20);
   print_array_float(stdout, depth, 20, 4, "%10.f ");
   return 0;
}
#endif