      program linfit
      real*8 sum,sumx,sumy,sumx2,sumxy,sumy2
      real*8 x1,y1,weight,delta,varnce
      real*8 x(1000),y(1000),sigmay(1000)
      real*8 a,sigmaa,b,sigmab,r
      integer*4 ndata,mode
      character data*80
c
   10 write(6,*)' Enter file name with x,y,dy !'
      read(5,1)data
    1 format(a80)
      open(1,file=data,status='old',err=10)
      do i=1,1000
       read(1,*,end=111)x(i),y(i),sigmay(i)
      enddo
  111 ndata=i-1
      write(6,*)'Enter mode!'
      write(6,*)'(+1: instrumental; 0: no weighting; -1: statistical)'
      read(5,*)mode
      call linefit(x,y,sigmay,ndata,mode,a,sigmaa,b,sigmab,r)
      write(6,*)'Fit is to: y = a + bx'
      write(6,*)' a=',a,'sigmaa=',sigmaa,' b=',b,' sigmab=',sigmab
      end
c---------------------------------------------------------------------
      subroutine linefit(x,y,sigmay,npts,mode,a,sigmaa,b,sigmab,r)
      real*8 sum,sumx,sumy,sumx2,sumxy,sumy2
      real*8 x1,y1,weight,delta,varnce
      real*8 x(*),y(*),sigmay(*)
      real*8 a,sigmaa,b,sigmab,r
      integer*4 c,npts,mode
   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
   53 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
     1 -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)/
     1 dsqrt(delta*(sum*sumy2-sumy*sumy))
      return
      end
