        program PBFT
C  PARABOLIC FITTING PROGRAM
C  Y=A+BX+CX**2
C  WEIGHTING:W(I)=1/SIGY(I)**2   ENTER SIGY(I)
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  w is the weight array
C  N IS THE NUMBER OF POINTS: INPUT DATA
C  R IS A three LINES ARRAY WHERE:R(1)=A,R(2)=B, R(3)=C:OUTPUT DATA
C  SIG2 IS THE VARIANCE OF THE FIT:SIG2=SUM((Y OBS-Y CALC)**2)*P/(N-2):
C    OUTPUT DATA
C  CM IS THE 3x3 COVARIANCE MATRIX:OUTPUT DATA
C************************************************************************
C  Written in 6 jun 86 S.O. Kepler
C  Last modified ** 6 jun 86 **
C************************************************************************
        implicit double precision(a-h,p-z)
        dimension x(1000),y(1000),w(1000),sigy(1000),ycalc(1000)
        dimension R(3),A(3,3),B(3),CM(3,3),sigma(3)
        DIMENSION P(20),Q(20),RR(20)
        CHARACTER*20 infile
  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)

10       DO 28 I=1,3
       DO 78 J=1,3
       A(I,J)=0.0d0
       B(I)=0.0d0
       R(I)=0.0d0
       CM(I,J)=0.0d0
78       CONTINUE
28       CONTINUE

c  read header***********************************

c  *****************************************
c  read in data - time (x), intensity (y) and sigma (sigy)
        DO 325 I=1,N
  325   READ(2,*) X(I),Y(I),sigy(i)
c  *****************************************
        iii=2
        write(*,521) iii

       DO 77 I=1,N
       XI=X(I)
       YI=Y(I)
       W(I)=1/SIGY(I)**2
       WI=W(I)
       A(1,2)=A(1,2)+WI*XI
       A(1,3)=A(1,3)+WI*XI**2
       A(2,3)=A(2,3)+WI*XI**3
       A(3,3)=A(3,3)+WI*XI**4
       B(1)=B(1)+WI*YI
       B(2)=B(2)+WI*YI*XI
       B(3)=B(3)+WI*YI*XI**2
       A(1,1)=A(1,1)+WI
77       CONTINUE
       A(2,1)=A(1,2)
       A(3,1)=A(1,3)
       A(2,2)=A(1,3)
       A(3,2)=A(2,3)
       NN=3
C       CALL SYMIN(A,NN)
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.0d0
        ONE = 1.0d0
        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 Y=R(1)+R(2)X+R(3)X**2
       DO 199 I=1,3
       DO 299 J=1,3
       R(I)=R(I)+A(I,J)*B(J)
299       CONTINUE
199       CONTINUE
       DO 1 I=1,N
       YCALC(I)=R(1)+R(2)*X(I)+R(3)*X(I)**2
1       CONTINUE
C  FINDING SIG2
38       CONTINUE
       RSIG=0.0d0
       DO 1333 I=1,N
       RSIG=RSIG+((Y(I)-R(1)-R(2)*X(I)-R(3)*X(I)**2)**2)*W(I)
1333       CONTINUE
       SIG2=RSIG/(N-3)
       SIG2NEW=RSIG/PTOTAL*N/(N-3)
       DO 1001 I=1,3
       DO 1002 J=1,3
       CM(I,J)=A(I,J)*SIG2NEW
1002       CONTINUE
1001       CONTINUE

        do 3915 i=1,3
        sigma(i)=dsqrt(cm(i,i))
3915    continue
100     format('Fitting a line y=a+bx+cx**2'/,'a = ',1pe12.5,/'b = ',1pe
     &12.5,/'c = ',1pe12.5,/'sig2= ',1pe12.3)

c... output
c
c  write header*******************************
c
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,3)
110     format('siga =',1pe12.3,/'sigb =',1pe12.3,/'sigc =',1pe12.3)
        write(3,107)((cm(i,j),j=1,3),i=1,3)
107     format(31x,'Correlation Matrix'/,
     &29x,'a, b, c'/,
     &1x,3(3(0pg18.6,2x)/,1x))
101     continue

c **********************
        stop
        end
