#include <math.h>
#include "IceDraft.h"
/*****************************************************************************
** int state(double *value,int ip, double pressure, double temperature, 
*		double salinity)
* 
*		this routine returns an error code and determines values for the state 
*		of seawater
*
*		ip = 1   density  - in kg/m**3
*		ip = 0   specific volume - m**3
*		ip = -1  secant bulk modulus
*		input :
*			pressure (gauge) in dbars
*			temperature in degrees celcius
*			salinity in    0/00
*
*		Formula from:
*		Millero, F. J., C. T. Chen, A. Bradshaw and K. Schleicher, A new high 
*		pressure equation of state for seawater, Deep-Sea Res., 27, 255-264, 1980.
*****************************************************************************/
int state(double *value, int ip,double pressure,double temperature,double salinity)
{
	int Error = 0;
	double bw,aw,kw,b,a,kst0,kstp,pw,pst01a,pst01b,pst0,vst0,vstp,pstp;
	
	// pressure in bars
	pressure = pressure/10.;

	bw = 8.50935e-5 + temperature * (-6.12293e-6 + temperature * 5.2787e-8);

	aw = 3.239908 + temperature * (1.43713e-3 + temperature * (1.16092e-4 - temperature * 5.77905e-7));

	kw = 19652.21 + temperature * (148.4206 + temperature * (-2.327105 + temperature * (1.360477e-2 + temperature * (-5.155288e-5))));

	b = bw + salinity * (-9.9348e-7 + temperature * (2.0816e-8 + temperature * 9.1697e-10));

	a = aw + salinity * (2.2838e-3 + temperature * (-1.0981e-5 - temperature * 1.6078e-6)) + 1.91075e-4 * pow(salinity,1.5);

	// Check for error
	if(salinity < 0)
		return ERROR_CTD_INVALID_SALINITY;
	
	kst0 = kw + salinity * (54.6746 + temperature * (-.603459 + temperature * (1.09987e-2 + temperature * (-6.167e-5)))) + (7.944e-2 + temperature * (1.6483e-2 + temperature * (-5.3009e-4))) * pow(salinity,1.5);

	kstp = kst0 + a * pressure + b * pressure * pressure;
	if(ip < 0 )
	{
		*value = kstp;
		return(Error);
	}
	
	pw = 999.842594 + temperature * (6.793952e-2 + temperature * (-9.09529e-3 + temperature * (1.001685e-4 + temperature * (-1.120083e-6 + temperature * 6.536332e-9))));

	pst01a = salinity * (8.24493e-1 + temperature * (-4.0899e-3 + temperature * (7.6438e-5 + temperature * (-8.2467e-7 + temperature * 5.3875e-9))));
	pst01b = (-5.72466e-3 + temperature * (1.0227e-4 + temperature * (-1.6546e-6))) *  pow(salinity,1.5) + 4.8314e-4 * salinity * salinity;
 
 	pst0 = pw + pst01a + pst01b;

	// Check for error
	if(pst0 == 0)
		return ERROR_CTD_PST0_EQ0_IN_STATE;
	if(kstp == 0)
		return ERROR_CTD_KSTP_EQ0_IN_STATE;
		
	vst0 = 1. / pst0;
	vstp = vst0 * (1. - pressure / kstp);
	pstp = pst0 / (1. - pressure / kstp);

	if(ip > 0 )
	{
		*value = pstp;
	}else{
		*value = vstp;
	}
		
		
	return Error;
		
}