subroutine linfit
      implicit double precision (a-h,o-z)
      common x,y,sigmay,npts,mode,a,sigmaa,b,sigmab,r
      dimension x(100),y(100),sigmay(100)
11    sum=0.
      sumx=0.
      sumy=0.
      sumx2=0.
      sumxy=0.
      sumy2=0.
21    do 50 i=1,npts
      xi=x(i)
      yi=y(i)
      if (mode) 31,36,38
31    if (yi) 34,36,32
32    weight=1./yi
      goto 41
34    weight=1./(-yi)
      goto 41
36    weight=1.
      goto 41
38    weight=1./sigmay(i)**2
41    sum=sum+weight
      sumx=sumx+weight*xi
      sumy=sumy+weight*yi
      sumx2=sumx2+weight*xi*xi
      sumxy=sumxy+weight*xi*yi
      sumy2=sumy2+weight*yi*yi
50    continue
51    delta=sum*sumx2-sumx*sumx
      a=(sumx2*sumy-sumx*sumxy)/delta
      b=(sumxy*sum-sumx*sumy)/delta
61    if (mode) 62,64,62
62    varnce=1.
      goto 67
64    c=npts-2
      varnce=(sumy2+a*a*sum+b*b*sumx2-2.*(a*sumy+b*sumxy-a*b*sumx))/c
67    sigmaa=dsqrt(varnce*sumx2/delta)
68    sigmab=dsqrt(varnce*sum/delta)
71    r=(sum*sumxy-sumx*sumy)/dsqrt(delta*(sum*sumy2-sumy*sumy))
      return
      end