/*
**  bvfreq.c 1.1	Solaris 2.3 Unix	940907	SIO/ODF	fmd
**
**  Calculate Brunt-Vaisala frequency (buoyancy frequency)
**  and stability (N**2).
**
**  Reference:
**  =========
**
**      Lynne Talley (SIO) from Bob Millard, Breck Owens &
**	N.P. Fofonoff (WHOI).
**
**      This implementation uses a Gaussian filter to smooth the
**      calculated density gradient.
**
**  Check Values:
**  ===== ======
**
**    for nObs=3,
**           s=35.0, 35.0, 35.0,
**           t=5.0, 4.5, 4.0,
**           p=1000, 1001, 1002:
**
**    BruntVaisalaFreq() = 14.5784 cph,
**    *pPav              = 1001 db, and
**    *pE                = 0.000647402 rad/sec.
**
*/
#include <math.h>
#include <odflib.h>

double				/* <- buoyancy freq (cycles/hour) */
BruntVaisalaFreq(nObs, s, sInc, t, tInc, p, pInc, pPav, pE)
	int	nObs;		/* -> number of points in s, t, p series */
	double	*s;		/* -> salinity series (PSUs) */
	int	sInc;		/* -> byte increment for salinity series */
	double	*t;		/* -> temperature series (deg C) */
	int	tInc;		/* -> byte increment for temperature series */
	double	*p;		/* -> pressure series (decibars) */
	int	pInc;		/* -> byte increment for pressure series */
	double	*pPav;		/* <- Mean pressure for interval (decibars) */
	double *pE;		/* <- stability (N**2) (radians/second) */
{
	int	count;
	static	double	*wts;
	static	int	wtCount;
	double	E, *pS, *pT, *pP, *pW,
		bvFreq, cX, cY, cXX, cXY, dP, dVdP, spVol, svAn, sigma;

	if (nObs<3) {
	    if (wts)
		(void)free((char *)wts);
	    *pPav = *pE = quiet_nan();
	    return (*pE);
	}
    /*
    **  Generate a Gaussian filter series, if we need to.
    */
	if (!wts || wtCount != nObs) {
	    if (wts)
		(void)free((char *)wts);
	    wts     = GaussianWeights(nObs);
	    wtCount = nObs;
	}
	if (!wts) {
	    *pPav = *pE = quiet_nan();
	    return (*pE);
	}
    /*
    **  Calculate the mean (filtered) pressure.
    */
	cX  = cY = cXX = cXY = 0.0;
	for (count = nObs, pW = wts, pP = p; count;
	    count--, cX += *pW++ * *pP, pP = (double *)((char *)pP+pInc));
	*pPav = cX;
    /*
    **  Calculate the density gradient.
    */
	for (count = nObs, pS = s, pT = t, pP = p, pW = wts; count;
	    count--, pW++,
	    pS = (double *)((char *)pS+sInc),
	    pT = (double *)((char *)pT+tInc),
	    pP = (double *)((char *)pP+pInc)) {
	    svAn = IESSpVolAn(*pS, PotentialTemp(*pS, *pT, *pP, cX),
			      cX, &sigma);
	    dP   = (*pP - cX);
	    cXY += *pW * svAn * dP;
	    cY  += *pW * svAn;
	    cXX += *pW * dP * dP;
	}
	dVdP  = cXY/cXX;
	spVol = (1.0/(sigma+1000.0))-svAn+cY;
    /*
    **  Stability (N**2) radians/second.
    */
	E = -0.0096168423*dVdP/(spVol*spVol);
	if (pE)
	    *pE   = E;
    /*
    **  Buoyancy frequency (N) cycles/hour.
    */
	bvFreq= 572.9578*sqrt(fabs(E));
	if (E < 0.0) bvFreq = -bvFreq;
	return (bvFreq);

} /* BruntVaisalaFreq() */