PROGRAM POLY
C
C     GENERAL POLYNOMIAL CURVE FITTING PROGRAM
C
C        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(10,10), COEF(10), XX(50), YY(50), ICAP(12), C(10), X(5
     10),DIFFC(50), TIT(6)
      CHARACTER TIT*13, TAB*1
      DIMENSION Y(50), B(10), YC(50)
C  Cricket Graph accepts columns separated by TAB characters.
      TAB= CHAR(9)
      TIT(1) = '      X(Obs)'//TAB
      TIT(2) = '      Y(Obs)'//TAB
      TIT(3) = '        f(X)'//TAB
      TIT(4) = '        f(Y)'//TAB
      TIT(5) = '     Y(Calc)'//TAB
      TIT(6) = 'Y(Cal)-Y(ob)'
C	Open files for UNIMAP, CricketGraph
      OPEN(3,FILE='cricket.dat',STATUS='UNKNOWN')
C
C     FIRST DATA CARD. MN= MAXIMUM DEGREE OF POLYNOMIAL, NN= NUMBER OF
C**** DATA POINTS, NWT= WEIGHTING(1 FOR RELATIVE ERRORS, 0 FOR ABSOLUTE
C**** ERRORS). ICAP CONTAINS TITLE OF PRINTOUT.
C
  10  WRITE(*,*)'Read Degree, No. Points, Weight(1=Rel,0=Abs)'
      READ (*,*,END=222) MN,NN,NWT
      IF (MN.EQ.0) GOTO 222
      IF (NN.EQ. 0) NN = 99
      IF (NN.LE.MN) GO TO 130
      WRITE(*,*)'READ VALUES FOR X AND Y'
      IC = 0
      DO 502 JJ=1,NN
       READ(*,*) XX(JJ), YY(JJ)
       IF(YY(JJ).EQ.0.0) GO TO 501
      IC = IC+1
  502 CONTINUE
  501 CONTINUE
      IF (NN.EQ.99) NN= IC
      WRITE(*,9002) NN
 9002 FORMAT('NUMBER OF POINTS IN SET IS ', I3)
      DO 20 JJ=1,NN
C
C     FUNCTION CARDS-APPROPRIATE FUNCTION OF THE VARIABLES IS FORMED HER
C
         X(JJ)=XX(JJ)
      Y(JJ) = YY(JJ)
   20 CONTINUE
         N = MN
C     DO 140 N=2,MN
         SUMYZ=0.0
         DO 30 J=1,N
            COEF(J)=0.0
            C(J)=0.0
         DO 30 K=1,N
   30    A(J,K)=0.0
         DO 70 JJ=1,NN
            WT=1.
            IF (NWT.NE.0) WT=1./(Y(JJ)*Y(JJ))
         DO 70 J=1,N
            IF (J.NE.1) GO TO 40
            COEF(J)=COEF(J)+WT*Y(JJ)
            GO TO 50
   40       COEF(J)=COEF(J)+(X(JJ)**(J-1))*WT*Y(JJ)
   50    DO 70 K=1,N
            IF ((K+J).NE.2) GO TO 60
            A(J,K)=A(J,K)+WT
            GO TO 70
   60       A(J,K)=A(J,K)+(X(JJ)**(K+J-2))*WT
   70    CONTINUE
         CALL MTRXIN1(A,N,WARN)
         DO 80 J=1,N
         DO 80 K=1,N
   80    C(J)=C(J)+A(J,K)*COEF(K)
         DO 100 JJ=1,NN
            WT=1.
            IF (NWT.NE.0) WT=1./(Y(JJ)*Y(JJ))
            YC(JJ)=C(1)
            DO 90 J=2,N
   90       YC(JJ)=C(J)*(X(JJ)**(J-1))+YC(JJ)
        DIFFC(JJ) = YC(JJ) - Y(JJ)
  100   SUMYZ = SUMYZ + (DIFFC(JJ)**2)*WT
         AN=NN-N
         S=SQRT(ABS(SUMYZ/AN))
 180     FORMAT(F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4)
 181     FORMAT(6F12.4)
         WRITE (*,170) N,S
         WRITE(3,171) ((TIT(J)),J=1,6)
 171     FORMAT(6A)
         WRITE(*,'(6(A12))') ((TIT(J)),J=1,6)
         WRITE (*,181) (XX(JJ),YY(JJ),XX(JJ),
     1   Y(JJ),YC(JJ),DIFFC(JJ),JJ=1,NN)
         WRITE (3,180) (XX(JJ),TAB,X(JJ),TAB,YY(JJ),TAB,
     1   Y(JJ),TAB,YC(JJ),TAB,DIFFC(JJ),JJ=1,NN)
         DO 110 J=1,N
            B(J)=(SQRT(ABS(A(J,J))))*S
            IJ=J-1
            WRITE (*,190) IJ,C(J),B(J)
  110    CONTINUE
C  140 CONTINUE
      GO TO 10
  130 WRITE (*,200)
      GO TO 10
C
 222  STOP
  160 FORMAT (2F12.5)
  170 FORMAT (///1X, 49HLINEAR LEAST SQUARES FIT TO POLYNOMIAL OF DEGRE
     1E ,I1,/1X,'S.E. OF Y IS ',E10.5/1X,30(2H -)//)
  190 FORMAT (/1X, 16HCOEFFICIENT OF X,I1,  3H = ,E12.5,  9H   +/-   ,E
     112.5)
  200 FORMAT (1X,36H INSUFFICIENT DATA-EXECUTION DELETED)
C
      END
      SUBROUTINE MTRXIN1(A,N,WARN)
      DIMENSION A(10,10), IPV(10,10)
      DO 10 J=1,N
   10 IPV(J,3)=0
      DO 110 I=1,N
         AMAX=0.
         DO 40 J=1,N
            IF (IPV(J,3).EQ.1) GO TO 40
            DO 30 K=1,N
               IF (IPV(K,3)-1) 20,30,150
   20          IF (AMAX.GE.ABS(A(J,K))) GO TO 30
               IROW=J
               ICOLUM=K
               AMAX=ABS(A(J,K))
   30       CONTINUE
   40    CONTINUE
C
C     ILL CONDITIONED MATRIX TEST NORMALISATION
C
         IF (I.GT.1) GO TO 60
         AMAX12=AMAX
         DO 50 J=1,N
         DO 50 K=1,N
   50    A(J,K)=A(J,K)/AMAX
   60    IPV(ICOLUM,3)=IPV(ICOLUM,3)+1
         IPV(I,1)=IROW
         IPV(I,2)=ICOLUM
         IF (IROW.EQ.ICOLUM) GO TO 80
         DO 70 LL=1,N
            SWAP=A(IROW,LL)
            A(IROW,LL)=A(ICOLUM,LL)
   70    A(ICOLUM,LL)=SWAP
   80    PIVOT=A(ICOLUM,ICOLUM)
         A(ICOLUM,ICOLUM)=1.0
         DO 90 LL=1,N
   90    A(ICOLUM,LL)=A(ICOLUM,LL)/PIVOT
      DO 110 L1=1,N
         IF (L1.EQ.ICOLUM) GO TO 110
         T=A(L1,ICOLUM)
         A(L1,ICOLUM)=0.0
         DO 100 LL=1,N
  100    A(L1,LL)=A(L1,LL)-A(ICOLUM,LL)*T
  110 CONTINUE
      DO 130 I=1,N
         LL=N-I+1
         IF (IPV(LL,1).EQ.IPV(LL,2)) GO TO 130
         JROW=(IPV(LL,1))
         JCOLUM=(IPV(LL,2))
         DO 120 K=1,N
            SWAP=A(K,JROW)
            A(K,JROW)=A(K,JCOLUM)
  120    A(K,JCOLUM)=SWAP
  130 CONTINUE
      DO 140 J=1,N
      DO 140 K=1,N
  140 A(J,K)=A(J,K)/AMAX12
      RETURN
  150 WRITE (6,160)
      RETURN
C
  160 FORMAT (1X,30H MATRIX IS SINGULAR-CHECK DATA/1X,60(2H *))
C
      END