/*******************************************************************
                    VECTORS.C

     Functions for working with one-dimensional arrays.
     Modified 10/17/87
********************************************************************/
#include "common.h"
#include "vector.h"

#ifndef M_PI
#define M_PI		3.14159265358979323846
#endif

void scale_elements(float scales[], float array[], float out_array[], int n )
{
     int i;
     float aa, ss;
     float bf = BADFLOAT * .999999;

     for (i=0; i<n; i++)
     {
          aa = array[i];
          ss = scales[i];
          if ( (aa>=bf) || (ss>=bf) )  out_array[i] = BADFLOAT;
          else  out_array[i] = aa * ss;
     }
}


/*********************************************
   |---------|---------|---------|---------|
  start                          |        end
  ar[0]     ar[1]    ......     ar[n-1]
**********************************************/

void make_uniform_grid(double start, double end, int n, float array[])
{
     int i;
     double incr;
     double grid_val;

     incr = (end - start) / n;
     array[0] = grid_val = start;
     for (i = 1; i < n; i++)
     {
          grid_val += incr;
          array[i] = grid_val;
     }
     return;
}

void make_sinusoid(float array[], int n, double amp, double cycles, double degrees)
{
     int i;
     double omega,
            phi;
     omega = 2 * M_PI * cycles / n;
     phi = degrees * M_PI / 180.0;
     for (i=0; i<n; i++)  array[i] = amp * cos( omega*i - phi );
     return;
}

/* Trapezoidal rule integration of a vector f specified on a grid z.
     All arrays have n elements.  integral[0] is always zero.
     There is no checking for BADFLOAT or other problems.
     The output (integral) array can be the same as the input (f);
     in other words, the integration can be done in place.
*/
void integrate_tr(float f[], float z[], float integral[], int n)
{
     int i;
     float fm, fi;
     fi = f[0];
     integral[0] = 0.0;
     for (i=1; i<n; i++)
     {
          fm = fi;
          fi = f[i];
          integral[i] = integral[i-1] + 0.5 * (fi + fm) * (z[i] - z[i-1]);
     }
     return;
}

/* Trapezoidal rule integration of a vector f specified on a grid z,
     used to calculate an average value.
     All arrays have n elements.
     There is no checking for BADFLOAT or other problems.
*/

double average_tr(float f[], float z[], int n)
{
     double integral = 0.0;
     int i;
     for (i=1; i<n; i++)
          integral += 0.5 * (f[i] + f[i-1]) * (z[i] - z[i-1]);
     return ( integral / (z[n-1] - z[0]) );
}

double scalar_product(float a[], float b[], int n)
{
     int i;
     double xx = 0.0;
     for (i=0; i<n; i++)  xx += a[i] * b[i];
     return (xx);
}

double mode_projection(float a[], float b[], int n)
{
     int i;
     double xx;

     xx = 0.5 * (a[0] * b[0] + a[n-1] * b[n-1]);
     for (i=1; i<n-1; i++)  xx += a[i] * b[i];
     return (xx);
}

void stats(float a[], int n, double *mean_ptr, double *variance_ptr)
{
     double s = 0.0, 
            ss = 0.0,
            ai;
     int i;

     for (i=0; i<n; i++)
     {
          ai = a[i];
          s += ai;
          ss += ai * ai;
     }
     *mean_ptr = s / n;
     *variance_ptr = ss / n - *mean_ptr * *mean_ptr;
     return;
}