/*******************************************************************
                    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

#if PROTOTYPE_ALLOWED
void scale_elements(float scales[], float array[], float out_array[], int n )
#else
void scale_elements( scales, array, out_array, n )
float scales[];
float array[];
float out_array[];   /* can be identical to array[] */
int n;
#endif
{
     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]
**********************************************/

#if PROTOTYPE_ALLOWED
void make_uniform_grid(double start, double end, int n, float array[])
#else
void make_uniform_grid( start, end, n, array )
double start;
double end;
int n;
float array[];
#endif
{
     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;
}

#if PROTOTYPE_ALLOWED
void make_sinusoid(float array[], int n, double amp, double cycles, double degrees)
#else
void make_sinusoid(array, n, amp, cycles, degrees)
float array[];
int   n;             /* number of elements; array length is n units */
double amp,           /* amplitude; should by positive */
       cycles,        /* frequency in cycles per array length */
       degrees;       /* phase in degrees of maximum past zero */
#endif
{
     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.
*/
#if PROTOTYPE_ALLOWED
void integrate_tr(float f[], float z[], float integral[], int n)
#else
void integrate_tr( f, z, integral, n)
float     f[], z[], integral[];
int       n;
#endif
{
     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.
*/

#if PROTOTYPE_ALLOWED
double average_tr(float f[], float z[], int n)
#else
double average_tr( f, z, n)
float     f[], z[];
int       n;
#endif
{
     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]) );
}

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

#if PROTOTYPE_ALLOWED
double mode_projection(float a[], float b[], int n)
#else
double mode_projection(a, b, n)
float a[], b[];
int   n;
#endif
{
     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);
}

#if PROTOTYPE_ALLOWED
void stats(float a[], int n, double *mean_ptr, double *variance_ptr)
#else
void stats(a, n, mean_ptr, variance_ptr)
float a[];
int n;
double *mean_ptr, *variance_ptr;
#endif
{
     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;
}