program llsl
C  THIS program FITS Y=A+BX
C  X IS THE INDEPENDENT VARIABLE ARRAY WITH N POINTS:INPUT DATA
C  Y IS THE DEPENDENT VARIABLE ARRAY WITH N POINTS:INPUT DATA
C  N IS THE NUMBER OF POINTS: INPUT DATA
C  R IS A two LINES ARRAY WHERE:R(1)=A,R(2)=B:OUTPUT DATA
C  SIG2 IS THE VARIANCE OF THE FIT:SIG2=SUM((Y OBS-Y CALC)**2)/(N-2)
C  CM IS THE 2X2 COVARIANCE MATRIX:OUTPUT DATA
C  USES SUBROUTINE SYMIN
C************************************************************************
C  Adapted by S. O. Kepler, this version written 16 DEC 85
C  Last modified ** 16 DEC 1985 **
C************************************************************************
        implicit double precision(a-h,p-z)
        common /x/ X(8000)
        common /y/ Y(8000)
        dimension R(2),A(2,2),B(2),CM(2,2),sigma(2)
        DIMENSION P(20),Q(20),RR(20)
        CHARACTER*20 infile
        pi=3.1415926536
  321   write(*,322)
  322   format(1x,'input file is: ', $)
        read(*,323) infile
  323   format(a20)
        open(2,file=infile,status='old')
c  ***********************************
c  read in the number of data points
        read(2,*)N
c  ***********************************
        iii=1
        write(*,521) iii
  521   format(1x,'program completed phase: ', I3)
        write(*,522) N
  522   format(1x,'N= ',i5)
        DO 1 I=1,2
        DO 1 J=1,2
        A(I,J)=0
        B(I)=0
        R(I)=0
        CM(I,J)=0
1       CONTINUE
c  read header***********************************

c  *****************************************
c  read in data - time (x) and intensity (y)
        DO 325 I=1,N
  325   READ(2,*) X(I),Y(I)
c  *****************************************
        iii=2
        write(*,521) iii
        DO 2 I=1,N
        XI=X(I)
        YI=Y(I)
        PI=1
        A(1,1)=A(1,1)+PI
        A(1,2)=A(1,2)+PI*XI
        A(2,2)=A(2,2)+PI*XI**2
        B(1)=B(1)+PI*YI
        B(2)=B(2)+PI*YI*XI
2       CONTINUE
        a(2,1)=a(1,2)
        NN=2

C        SUBROUTINE SYMIN (A,NN)
C
C        THIS SUBROUTINE INVERTS A SYMMETRIC MATRIX.
C        ***** INPUT DATA *****
C        A  =  A SYMMETRIC MATRIX
C        NN  =  THE RANK OF THE MATRIX. N MUST BE LESS THAN 20.
C        ***** OUTPUT DATA *****
C        A  =  THE INVERSE OF THE INPUT MATRIX
C
C        DIMENSION A(NN,NN)
C        DIMENSION P(20),Q(20),R(20)
1300    FORMAT(13H0SYMIN FAILED)
        ZERO = 0.0
        ONE = 1.0
        DO 1000 M=1,NN
1000    RR(M) = ONE
        DO 200 M=1,NN
        BIG = ZERO
        DO 3000 L=1,NN
        AB = DABS(A(L,L))
        IF(AB-BIG)3000,3000,4000
4000    IF(RR(L))1400,3000,1400
1400    BIG = AB
        K = L
3000    CONTINUE
        IF(BIG)6000,5000,6000
5000    WRITE (*,1300)
        GO TO 101
C       RETURN
6000    RR(K) = ZERO
        Q(K) = ONE/A(K,K)
        P(K) = ONE
        A(K,K) = ZERO
        KM1 = K-1
        IF(KM1.EQ.0) GO TO 1600
        DO 7000 L=1,KM1
        P(L) = A(L,K)
        IF(RR(L))9000,8000,9000
8000    Q(L) = A(L,K)*Q(K)
        GO TO 7000
9000    Q(L) = -A(L,K)*Q(K)
7000    A(L,K) = ZERO
1600    CONTINUE
        KP1 = K+1
        IF(KP1.GT.NN) GO TO 1700
        DO 1500 L=KP1,NN
        IF(RR(L))1200,1100,1200
1200    P(L) = A(K,L)
        GO TO 10000
1100    P(L) = -A(K,L)
10000   Q(L) = (-A(K,L))*Q(K)
1500    A(K,L) = ZERO
1700    CONTINUE
        DO 200 L=1,NN
        DO 200 K=L,NN
200     A(L,K) = A(L,K) + P(L)*Q(K)
        M = NN+1
        L = NN
        DO 27 K=2,NN
        M = M-1
        L = L-1
        DO 27 J=1,L
27      A(M,J) = A(J,M)
C       RETURN
C  NOW THE A MATRIX IS THE INVERSE OF THE INITIAL ONE
C  FINDING R(I),WHERE R(1)=A,R(2)=B
        DO 4 I=1,2
        DO 4 J=1,2
        R(I)=R(I)+A(I,J)*B(J)
4       CONTINUE
C  FINDING SIG2
        RSIG=0
        DO 5 I=1,N
        PI=1
        RSIG=RSIG+((Y(I)-R(1)-R(2)*X(I))**2)*PI
5       CONTINUE
        SIG2=RSIG/(N-2)
C  FINDIND THE CORRELATION MATRIX CM
        DO 6 I=1,2
        DO 6 J=1,2
        CM(I,J)=SIG2*A(I,J)
6        CONTINUE
        do 3915 i=1,2
        sigma(i)=dsqrt(cm(i,i))
3915    continue
100     format('Fitting a line y=a+bx'/,'a = ',1pe12.5,/'b = ',1pe12.5,
     &/'sig2= ',1pe12.3)
c  write header*******************************
c  *******************************************
  432   write(*,433) 
  433   format(1x,'output file name is: ', $)
        read(*,323) infile
        open(3,file=infile,status='new')
        write(3,99)N
99      format('npoints= ',i5)

c  output results
        write(3,100)r,sig2
        write(3,110)(sigma(i),i=1,2)
110     format('siga ='1pe12.3,/'sigb ='1pe12.3)
        write(3,107)((cm(i,j),j=1,2),i=1,2)
107     format(31x,'Correlation Matrix'/,
     &29x,'a, b'/,
     &1x,2(2(0pg18.6,2x)/,1x))
101     continue
        end