/*
**  soundvel.c	1.1	Solaris 2.3 Unix	940906	SIO/ODF	fmd
**
**  Sound velocity in sea water calculations:
**
**  double SoundVelocity(int method, double s, double t, double p)
**	calculate the velocity of sound in seawater.
**  double SoundVelCandM(double s, double t, double p)
**	Chen & Millero sound velocity.
**  double SoundVelDelGrosso(double s, double t, double p)
**	Del Grosso sound velocity.
**  double SoundVelWilson(double s, double t, double p)
**	Wilson sound velocity.
**
**  Plus soon-to-be obsolete Fortran interface routines.
**
**  References:
**  ==========
**
**	Chen & Millero, 1977,JASA,62,1129-1135.
**	Del Grosso, NLR ii, Navel Research Lab 1974.
**	Tolstoy & Clay (Wilson), Ocean Acoustics 1966, pg. 6.
*/
#include <math.h>
#include <odflib.h>

#ifndef FORTRAN

/*
**  Calculate the velocity of sound in seawater.
*/
double					/* <- sound velocity (m/sec) */
SoundVelocity(method, s, t, p)
	int	method;			/* -> calculation method:
						0 = Chen & Millero,
						1 = Del Grosso,
						2 = Wilson */
	double	s;			/* -> salinity (psu) */
	double	t;			/* -> temperature (degrees C) */
	double	p;			/* -> pressure (decibars) */
{
	switch (method) {
	case 1:
	    return (SoundVelDelGrosso(s, t, p));
	case 2:
	    return (SoundVelWilson(s, t, p));
	default:
	    return (SoundVelCandM(s, t, p));
	}
} /* SoundVelocity() */

double
SoundVelCandM(s, t, p)
	double	s;
	double	t;
	double	p;
{
	double	a, a0, a1, a2, a3, b, b0, b1, c, c0, c1, c2, c3, d, sr;

	p *= 0.1;			/* pressure in bars */
	sr = sqrt(fabs(s));		/* s**2 term */
	d  = 1.727e-3 - 7.9836e-6*p;
    /*
    **  s**(3/2) term
    */
	b1 = 7.3637e-5 +
	     1.7945e-7*t;
	b0 = -1.922e-2 -
	      4.42e-5*t;
	b  = b0 + b1*p;
    /*
    **  s**1 term
    */
	a3 = (-3.389e-13*t    +
	       6.649e-12)*t   +
	       1.10e-10;
	a2 = ((7.988e-12*t    -
	       1.6002e-10)*t  +
	       9.1041e-9)*t   -
	       3.9064e-7;
	a1 = (((-2.0122e-10*t +
		 1.0507e-8)*t -
		 6.4885e-8)*t -
		 1.258e-5)*t  +
		 9.4742e-5;
	a0 = (((-3.21e-8*t    +
		 2.006e-6)*t  +
		 7.164e-5)*t  -
		 1.262e-2)*t  +
		 1.389;
	a  = ((a3*p + a2)*p + a1)*p + a0;
    /*
    **  s**0 term
    */
	c3 = (-2.3643e-12*t   +
	       3.8504e-10)*t  -
	       9.7729e-9;
	c2 = (((1.0405e-12*t  -
		2.5335e-10)*t +
		2.5974e-8)*t  -
		1.7107e-6)*t  +
		3.1260e-5;
	c1 = (((-6.1185e-10*t +
		 1.3621e-7)*t -
		 8.1788e-6)*t +
		 6.8982e-4)*t +
		 0.153563;
	c0 = ((((3.1464e-9*t  -
		 1.47800e-6)*t+
		 3.3420e-4)*t -
		 5.80852e-2)*t+
		 5.03711)*t   +
		 1.402388e3;
	c  = ((c3*p + c2)*p + c1)*p + c0;
	return (c + (a + b*sr + d*s)*s);
} /* SoundVelCandM() */

double
SoundVelDelGrosso(s, t, p)
	double	s;
	double	t;
	double	p;
{
	double	c0, dcstp, dltacp, dltacs, dltact;

	p      *= 0.1019716;	/* pressure in kg/cm**2 gauge */
	c0      = 1402.392;

	dltact  = t*( 0.501109398873e1  +
		  t*(-0.550946843172e-1 +
	          t*  0.221535969240e-3));

	dltacs  = s*(0.132952290781e1  +
		  s* 0.128955756844e-3);

	dltacp  = p*( 0.156059257041e0 +
		  p*( 0.24499868841e-4 +
	          p*(-0.883392332513e-8)));
	dcstp   =  t*(-0.127562783426e-1*s      + 
		   p*( 0.635191613389e-2        +
		   p*( 0.265484716608e-7*t      +
		      -0.159349479045e-5        +
		       0.522116437235e-9*p)     +
		     (-0.438031096213e-6)*t*t)) +
	          s*((-0.161674495909e-8)*s*p*p +
		   t*( 0.968403156410e-4*t      +
		   p*( 0.485639620015e-5*s      +
		     (-0.340597039004e-3))));
	return (c0 + dltact + dltacs + dltacp + dcstp);
} /* SoundVelDelGrosso() */

double
SoundVelWilson(s, t, p)
	double	s;
	double	t;
	double	p;
{
	double	c0, vp, vs, vstp, vt;

	p    = (p+10.13)*0.1019716;
	s   -= 35.0;
	c0   = 1449.14;
	vt   = t*( 4.5721    +
	       t*(-4.4532e-2 +
	       t*(-2.6045e-4 +
	       t*  7.9851e-6)));
	vp   = p*( 0.160272   +
	       p*( 1.0268e-5  +
	       p*( 3.5216e-9  +
	       p*(-3.3603e-12))));
	vs   = s*(1.39799     +
	       s* 1.69202e-3);
	vstp = s*(t*(-1.1244e-2   +
		  p*  3.158e-8    +
		  t*( 7.7711e-7   +
		  p*  1.579e-9))  +
		  p*( 7.7016e-5   +
		  p*(-1.2943e-7)))+
		p*t*(-1.8607e-4   +
		  t*( 7.4812e-6   +
		  t*  4.5283e-8)  +
		  p*(-2.5294e-7   +
		  t*  1.8563e-9   +
		  p*(-1.9646e-10)));
	return (c0 + vt + vp + vs + vstp);
} /* SoundVelWilson() */

#else  /* FORTRAN */

double
svel_(select, p, t, s)
	int	*select;
	float	*p, *t, *s;
{
	return (SoundVelocity(*select, *s, *t, *p));
}
#endif /* FORTRAN */