#include <math.h>
#include "IceDraft.h"
/*****************************************************************************
** int pss78(double *salinity, double pressure, double temperature, double conductivity)
*
*  this routine returns an error code and computes a salinity value given: 
* 
* 	pressure (gauge) in dbars
* 	temperature in degrees celsius
* 	conductivity ratio  c(s,t,p)/c(35,15,0)
* 
* 	formula from:
*	1978 PRACTICAL SALINITY SCALE EQUATIONS, from IEEE Journal of Oceanic Engineering, 
*	Vol. OE-5, No. 1, January 1980, page 14.
*****************************************************************************/

int pss78(double *salinity, double pressure, double temperature, double conductivity)
{
	double alpha,r0,rt,rs;
	int Error = 0;

	alpha = pressure * ( 2.07e-5 + pressure * (-6.37e-10 + pressure * (3.989e-15)))/(1. + temperature * (3.426e-2 + temperature * 4.464e-4 - conductivity * 3.107e-3) + conductivity * 4.215e-1);

	// Check for error
	if(alpha == -1)
		return ERROR_CTD_ALPHA_EQNEG1_IN_PSS78;
	
	r0 = conductivity / (1. + alpha);

	rt = r0/(.6766097 + temperature * (2.00564e-2 + temperature * (1.104259e-4 + temperature * ( -6.9698e-7 + temperature * 1.0031e-9))));

	// Check for error
	if(rt < 0)
		return ERROR_CTD_RT_LESS0_IN_PSS78;
		
	rs = sqrt(rt);
	
	*salinity = (.008 + rs * (-.1692 + rs * (25.3851 + rs * (14.0941 + rs * (-7.0261 + rs * 2.7081))))) + (temperature - 15.) * (.0005 + rs * (-.0056 + rs * (-.0066 + rs * (-.0375 + rs * (.0636 - rs * .0144))))) / (1. + .0162 * (temperature - 15.));
	
	if(*salinity < 0)
		return ERROR_CTD_INVALID_SALINITY;
		
	return(Error);
	
}