/* nslope3.c - Subroutine to calculate nettleton slope statistics.

	This new version does not remove a mean or median, and uses
	weighted standard deviations to calculate the slope.  */

#include <stdio.h>
#include <math.h>

void nslope3_(grav,topo,weight,ns,wsum,gbar,sigg,tbar,sigt,r,prob)
float	grav[], topo[], weight[];
int	*ns;
float	*wsum, *gbar, *sigg, *tbar, *sigt, *r, *prob;
{
	/* Subroutine to calculate nettleton slope statistics.
	 *
	 * Input:
	 * grav   - vector of gravity values
	 * topo   - vector of topography values
	 * weight - vector of weights associated with these values
	 * ns     - number of points in grav, topo, weight
	 *
	 * Output:
	 * gbar   - median gravity (set equal to zero)
	 * sigg   - standard deviation of gravity
	 * tbar   - median topography (set equal to zero)
	 * sigt   - standard deviation of topography
	 * r      - ordinary linear correlation coefficient
	 * prob   - Probability that r this big indicates true correlation
	 * wsum	  - Sum of weights on input.
	 *
	 * Note:  All values are called by adress; this is FORTRAN callable.
	 *
	 * W. H. F. Smith, 25 October, 1993.   */

	int	i;
	double	gsum, tsum, gtsum, rho, w, z, ws;

	*wsum = 0.0;
	gsum = tsum = gtsum = 0.0;
	for (i = 0; i < *ns; i++) {
                ws  = weight[i];
		gsum += (grav[i] * grav[i])*ws;
		tsum += (topo[i] * topo[i])*ws;
		gtsum += (grav[i] * topo[i])*ws;
		*wsum += weight[i];
	}
        if (*wsum < 4.) {
           *gbar=0.;
           *sigg=0.;
           *tbar=0.;
           *sigt=0.;
           *r=0.;
           *prob=0.;
        }
        else {
           rho = gtsum/sqrt(gsum * tsum);
           w = 0.5 * (log(1.0 + rho) - log(1.0 - rho));
           z = w * sqrt((double)(*wsum - 3));
           *prob = (float)erf(0.707106781 * fabs(z));
           *r = (float)rho;
           *gbar = 0.0;
           *tbar = 0.0;
           *sigg = (float)sqrt(gsum/ *wsum);
           *sigt = (float)sqrt(tsum/ *wsum);
        }
	return;
}