#include <stdio.h>
#include <math.h>
#include "IceDraft.h"
/*****************************************************************************
** double sound(double *soundspeed, double pressure, double temperature,
* 	double salinity)
* 
*  This routine returns an error code and computes a sound speed in m/sec given:
*
* 	Pressure (gauge) in Dbars
* 	Temperature in C degrees
* 	Salinity in s in parts/1000 (ppt)
*
* 	Formula from:
* 	Del Grosso V.A., "New Equations for the Speed of Sound in Natural Waters
*	(with Comparisons to other Equations".  The Journal of the Acoustical Society
*	of America, v. 56, pp. 1084-1091. (1974)
*****************************************************************************/

#define SQUARE(x) (x*x)
#define CUBE(x)   (x*x*x)

int ips_sound(double *soundspeed, double pressure,double temperature,double salinity)
{
	double p,s,t,co,ct,cs,cp,cstp;
	int Error = 0;
	
	if(salinity == 0.) 
		return(ERROR_CTD_INVALID_SALINITY);

	// convert pressure in decibars to kg/cm2 
	p = pressure / GRAVITY;
	s = salinity;
	t = temperature;
	co = 1402.392;
	ct = (0.501109398873e1*t)-(0.550946843172e-1*SQUARE(t));
	ct += (0.221535969240e-3*CUBE(t));
	cs = (0.132952290781e1*s)+(0.128955756844e-3*SQUARE(s));
	cp = (0.156059257041*p)+(0.244998688441e-4*SQUARE(p));
	cp -= (0.883392332513e-8*CUBE(p));
	cstp  = (-0.127562783426e-1*t*s)+(0.635191613389e-2*t*p);
	cstp += (0.265484716608e-7*(SQUARE(t))*(SQUARE(p)));
	cstp -= (0.159349479045e-5*t*SQUARE(p));
	cstp += (0.522116437235e-9*t*CUBE(p));
	cstp -= (0.438031096213e-6*(CUBE(t))*p);
	cstp -= (0.161674495909e-8*(SQUARE(s))*(SQUARE(p)));
	cstp += (0.968403156410e-4*(SQUARE(t))*s);
	cstp += (0.485639620015e-5*t*(SQUARE(s))*p);
	cstp -= (0.340597039004e-3*t*s*p);
	*soundspeed = co + ct + cs + cp + cstp;
	
	return(Error);
	
	
	
}