c-----------------------------------------------------------------------
c
c  Name:
c      SUBROUTINE LINREG
c
c  Purpose:
c      Compute line of best fit to a dataset defined in X and Y.
c
c  Usage:
c      CALL LINREG( N, X, Y, XBAR, YBAR, XSIG, YSIG, YINT,
c    &   SLOPE, CORR, EVAR, SERR )
c 
c  Input:
c      N       Number of points in X and Y arrays
c      X       X ordinate array
c      Y       Y coordinate array
c
c  Output:
c      XBAR    Mean of X values
c      YBAR    Mean of Y values
c      XSIG    Standard deviation of X values
c      YSIG    Standard deviation of Y values
c      YINT    Intercept on the Y axis of the line of best fit
c      SLOPE   Slope of the line of best fit
c      CORR    Correlation coefficient
c      EVAR    Variance
c      SERR    Standard error
c
c  Revised:
c      Liam.Gumley@ssec.wisc.edu
c      $Id: linreg.f,v 1.3 1999/11/11 20:04:45 gumley Exp $
c
c  Notes:
c      All output parameters are DOUBLE PRECISION.
c
c-----------------------------------------------------------------------

      subroutine linreg( n, x, y, xbar, ybar, xsig, ysig, yint,
     &                   slope, corr, evar, serr )

      implicit none
      
c ... arguments

      integer n
      real x( * ), y( * )
      double precision xbar, ybar, xsig, ysig, yint, slope, corr,
     &  evar, serr  
      
c ... local variables
      
      integer i
      double precision sn, cvar, xp, xvar, yp, yvar
      
      sn = dble( n )
      xbar = 0.0d0
      ybar = 0.0d0
      do i = 1, n
         xbar = xbar + dble( x( i ) )
         ybar = ybar + dble( y( i ) )
      end do
      xbar = xbar / sn
      ybar = ybar / sn

      xvar = 0.0d0
      yvar = 0.0d0
      cvar = 0.0d0
      do i = 1, n
         xp = dble( x( i ) ) - xbar
         yp = dble( y( i ) ) - ybar
         xvar = xvar + xp * xp
         yvar = yvar + yp * yp
         cvar = cvar + xp * yp
      end do

      xvar = xvar / sn
      xsig = sqrt( xvar )
      yvar = yvar / sn
      ysig = sqrt( yvar )
      cvar = cvar / sn
      slope = cvar / xvar
      yint = ybar - slope * xbar
      corr = cvar / ( xsig * ysig )
      evar = slope * cvar
      serr = sqrt( abs( yvar - evar ) )
      evar = evar / yvar
      
      end