/*
**  grav.c 1.1	Solaris 2.3 Unix	940906	SIO/ODF	fmd
**
**  Calculate local gravity for a given latitude and pressure.
**
**  Reference:
**  =========
**
**	Anon (1970) Geodetic Reference System 1967. Bulletin Geodesique,
**	Special Publication of the Association Internationale
**	de Geodesie, Paris
*/
#include <math.h>
#include <odflib.h>

#ifndef FORTRAN
double				/* <- local earth gravity (m/sec**2) */
Gravity(deglat, pdb)
	double	deglat;		/* -> degrees latitude */
	double	pdb;		/* -> pressure (decibars) */
{
	double	slatsq, gmss;

	slatsq  = sin(deglat*0.0174532925);
	slatsq *= slatsq;
	gmss    = ((2.36e-5    * slatsq   +
		    0.0052788) * slatsq   +
		    1.0)       * 9.780318 +
		    2.184e-6   * pdb;
	  /* gravity acceleration m/sec**2 */
	return (gmss);
} /* Gravity() */

#else  /* FORTRAN */

void
grav_(deglat, pdb, gmss)
	float	*deglat, *pdb, *gmss;
{
	*gmss = Gravity(*deglat, *pdb);
}
#endif /* FORTRAN */