      subroutine linfit ( npts, x, y, w, p, chisq )

      implicit none

      integer*4 npts
      real*4 x(*), y(*), w(*), p(2), chisq
      integer*4 i
      real*4 S0, Sx, Sy, Sxx, Syy, Sxy

      S0 = 0
      Sx = 0
      Sy = 0
      Sxx = 0
      Syy = 0
      Sxy = 0

      do i=1,npts
	S0 = S0 + w(i)
	Sx = Sx + w(i)*x(i)
	Sy = Sy + w(i)*y(i)
	Sxx = Sxx + w(i)*x(i)**2
	Syy = Syy + w(i)*y(i)**2
	Sxy = Sxy + w(i)*x(i)*y(i)
      end do

      p(1) = (Sy*Sxx - Sx*Sxy)/(S0*Sxx - Sx*Sx)
      p(2) = (S0*Sxy - Sx*Sy)/(S0*Sxx - Sx*Sx)

C      write (*,*) "S0,Sx,Sy,Sxx,Syy,Sxy = ",S0,Sx,Sy,Sxx,Syy,Sxy

C -- And a little something extra from the derivation:

      chisq = p(1)*p(1)*S0 + 2.0*p(1)*p(2)*Sx + p(2)*p(2)*Sxx 
     >      - 2.0*p(1)*Sy - 2.0*p(2)*Sxy + Syy

      return

      end
