#* HISTORY *
#* B.Simon	09-Nov-92	converted to spp from asurv 1.2

# REGRES -- Linear regression routine

# This subprogram calculates linear regression coeffs., variance, 
# and correlation coeffs., based on Philip R. Bevington, "Data Reduction 
# and Error Analysis for the Physical Sciences", McGraw Hill (NY:NY), 
# program 9-1 of P. 172

procedure regres (x, y, ntot, nvar, alpha, sigmaa, array)

double	x[nvar,ntot]		# i: independent variable array
double	y[ARB]			# i: dependent variable array
int	ntot			# i: number of data points
int	nvar			# i: number of independent variables
double	alpha[ARB]		# o: regression coefficients
double	sigmaa[ARB]		# o: std dev of coeeficients
double	array[nvar,nvar]	# o: correlation coefficient matrix
#--
double	det, a1, ymean, sigma, varnce, fnpts, free1, freen, freej
int	i, j, k
pointer	sp, a, r, xmean, sigmax, yfit

string	badmatrix  "%4w Warning : Regression matrix singular, results invalid"

begin
	# Allocate dynamic memory for temporary arrays

	call smark (sp)
	call salloc (a, nvar, TY_DOUBLE)
	call salloc (r, nvar, TY_DOUBLE)
	call salloc (xmean, nvar, TY_DOUBLE)
	call salloc (sigmax, nvar, TY_DOUBLE)
	call salloc (yfit, ntot, TY_DOUBLE)

	# Initialize degrees of freedom

	fnpts = ntot
	free1 = ntot - 1
	freen = ntot - nvar
	freej = nvar

	# Compute means of independent and dependent variables

	ymean = 0.0
	call aclrd (Memd[xmean], nvar)
	do i = 1, ntot {
	    ymean = ymean + y[i]
	    do j = 1, nvar {
		Memd[xmean+j-1] = Memd[xmean+j-1] + x[j,i]
	    }
	}

	ymean = ymean / fnpts
	do j = 1, nvar {
	    Memd[xmean+j-1] = Memd[xmean+j-1] / fnpts
	}

	# Accumulate matrices r and array

	sigma = 0.0
	call aclrd (Memd[r], nvar)
	call aclrd (Memd[sigmax], nvar)
	call aclrd (array, nvar * nvar)

	do i = 1, ntot {
	    sigma = sigma + (y[i] - ymean) ** 2

	    do j = 1,nvar {
		Memd[sigmax+j-1] = Memd[sigmax+j-1] + 
				   (x[j,i] - Memd[xmean+j-1]) ** 2
		Memd[r+j-1] = Memd[r+j-1] + (x[j,i] - Memd[xmean+j-1]) *
			      (y[i]-ymean)

		do k = 1, j {
		    array[j,k] = array[j,k] + (x[j,i] - Memd[xmean+j-1]) * 
				 (x[k,i] - Memd[xmean+k-1])
		}
	    }
	}

	sigma = sqrt (sigma / free1)

	do j = 1, nvar {
	    Memd[sigmax+j-1] = sqrt (Memd[sigmax+j-1] / free1)
	    Memd[r+j-1] = Memd[r+j-1] / (free1 * Memd[sigmax+j-1] * sigma)
	    do k = 1, j {
		array[j,k] = array[j,k] / 
			     (free1 * Memd[sigmax+j-1] * Memd[sigmax+k-1])
		array[k,j] = array[j,k]
	    }
	}

	# Invert symmetric matrix

	call matinv (array, nvar, det)
	if (det == 0.0) {
	    call aclrd (alpha, nvar+1)
	    call eprintf (badmatrix)
	    return
	}

	# Calculate coefficients

	a1 = ymean
	call aclrd (Memd[a], nvar)
	call aclrd (Memd[yfit], ntot)

	do j = 1, nvar {
	    do k = 1, nvar {
		Memd[a+j-1] = Memd[a+j-1] + Memd[r+k-1] * array[j,k]
	    } 

	    Memd[a+j-1] = Memd[a+j-1] * sigma / Memd[sigmax+j-1]
	    a1 = a1 - Memd[a+j-1] * Memd[xmean+j-1]

	    do i = 1, ntot {
		Memd[yfit+i-1] = Memd[yfit+i-1] + Memd[a+j-1] * x[j,i]
	    }
	}

	# Calculate uncertainties

	varnce = 0.0

	do i = 1, ntot {
	    Memd[yfit+i-1] = Memd[yfit+i-1] + a1
	    varnce = varnce + (y[i] - Memd[yfit+i-1]) ** 2
	}

	varnce = varnce / freen

	do j = 1, nvar {
	    sigmaa[j+1] = array[j,j] * varnce / (free1 * Memd[sigmax+j-1] ** 2)
	    sigmaa[j+1] = sqrt (sigmaa[j+1])
	}

	sigmaa[1] = varnce / fnpts

	do j = 1, nvar {
	    do k = 1, nvar {
		sigmaa[1] = sigmaa[1] + varnce * Memd[xmean+j-1] * 
			    Memd[xmean+k-1] * array[j,k] / (free1 * 
			    Memd[sigmax+j-1] * Memd[sigmax+k-1])
	    }
	}

	sigmaa[1] = sqrt (sigmaa[1])

	# Copy regression coefficients to output array

	alpha[1] = a1
	do j = 1, nvar {
	    alpha[j+1] = Memd[a+j-1]
	}

	call sfree (sp)
	return
end