      PROGRAM FIT                                                                                                                   
C     PROGRAM FIT (SERIES,       INPUT,       OUTPUT,                   )       H2                                                  
C    1       TAPE1=SERIES, TAPE5=INPUT, TAPE6=OUTPUT)                                                                               
C     -----------------------------------------------                           H2                                                  
C                                                                               H2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       H2                                                  
C     JUL 22/92 - E. CHAN  (REPLACE UNFORMATTED I/O WITH I/O ROUTINES)                                                              
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII)                                                                    
C     SEP 13/85 - B.DUGAS.                                                                                                          
C                                                                               H2                                                  
CFIT     - TIME SERIES POLYNOMIAL REGRESSION FIT,DISPLAY & STAT. TEST   1    C  H1                                                  
C                                                                               H3                                                  
CAUTHOR  - B. DUGAS                                                             H3                                                  
C                                                                               H3                                                  
CPURPOSE - FIT DOES A SERIES OF POLYNOMIAL REGRESION AND STATISTICAL            H3                                                  
C          TESTING THEREOF. SEE (BENNET & FRANKLIN, 1962, JOHN WILEY            H3                                                  
C          PUB., P. 244, 255-257) FOR EXACT REFERENCE ON THE METHOD.            H3                                                  
C          IT ALSO PLOTS THE RESULTS.                                           H3                                                  
C          NOTE - THE PLOTTED VALUES HAVE THE MEAN OF THE INPUT RECORD          H3                                                  
C                 REMOVED FROM THEM. BUT THIS REMOVED VALUE IS ALSO             H3                                                  
C                 WRITTEN TO OUTPUT...                                          H3                                                  
C                                                                               H3                                                  
CINPUT FILE...                                                                  H3                                                  
C                                                                               H3                                                  
C      SERIES = A SET OF TIME SERIES SHORTER THAN "TSL" POINTS EACH.            H3                                                  
C                                                                                                                                   
CINPUT PARAMETERS...                                                                                                                
C                                                                               H5                                                  
C      K = THE MAXIMUM DEGREE OF THE POLYNOMIAL FIT TO THE DATA. THIS           H5                                                  
C          IS DEFAULTED TO 2. THIS NUMBER IS SMALLER OR EQUAL TO 6.             H5                                                  
C                                                                               H5                                                  
C      N = THE NUMBER OF BUNCHED POINTS TO BE USED IN THE ANALYSIS AND          H5                                                  
C          DISPLAY. THE DEFAULT VALUE OF N IS THE ACTUAL LENGTH OF ONE          H5                                                  
C          DATA RECORD. THIS IS SMALLER OR EQUAL TO "IJ".                       H5                                                  
C                                                                               H5                                                  
CEXAMPLE OF INPUT CARD...                                                       H5                                                  
C                                                                               H5                                                  
C*FIT         3  365                                                            H5                                                  
C-------------------------------------------------------------------------                                                          
C                                                                                                                                   
                                                                                                                                    
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      DIMENSION F(18528),X(18528),DX(18528),R(4)                                                                                    
      DIMENSION DR(6),A(6),P(18528,6),T(6)                                                                                          
      COMMON /ICOM/ IBUF(8),IDAT(18528)                                                                                             
      LOGICAL OK                                                                                                                    
                                                                                                                                    
      DATA MAXX/18528/                                                                                                              
C--------------------------------------------------------------                                                                     
                                                                                                                                    
      NF = 3                                                                                                                        
      CALL JCLPNT(NF,1,5,6)                                                                                                         
      REWIND 1                                                                                                                      
      UN = 1.0E0                                                                                                                    
                                                                                                                                    
C     * READ INPUT PARAMETERS K, N. DEFAULT VALUE OF N IS THE                                                                       
C     * ACTUAL LENGTH OF ONE DATA RECORD (HAS TO BE < 1001).                                                                        
                                                                                                                                    
      READ (5,'(10X,2I5)',END=900) K, N                                         H4                                                  
      IF (K.EQ.0) K = 3                                                                                                             
      IF (N.EQ.0)   THEN                                                                                                            
          CALL FBUFFIN(1,IBUF,-8,KK,LEN)                                                                                            
          IF (KK.GE.0) GOTO 999                                                                                                     
          N = IBUF(5)*IBUF(6)                                                                                                       
          REWIND 1                                                                                                                  
      ENDIF                                                                                                                         
      WRITE(6,'("0 VALUES OF K, N ARE ",I5,1X,I5)') K,N                                                                             
      IF (N  .LT.4 .OR. N.GT.1000 .OR. K.LT.1 .OR.                                                                                  
     1    N-K.LT.2 .OR. K.GT.5    )  CALL          XIT('FIT',-1)                                                                    
                                                                                                                                    
C     * CALCULATE THE P(I,J), I=1,N AND J=0,K WHEN                                                                                  
C     *                                                                                                                             
C     *                              (N**2-K**2)                                                                                    
C     *   P(Z,K+1) = P(Z,K) - K**2*---------------*P(Z,K-1) .                                                                       
C     *                             4*(4*K**2 -1)                                                                                   
                                                                                                                                    
      C1 = FLOAT(N*N-1)/12.0E0                                                                                                      
      C2 = FLOAT(3*N*N-7)/20.0E0                                                                                                    
      C3 = FLOAT(9*N*N-81)/140.0E0                                                                                                  
      C4 = C2+C3                                                                                                                    
      C5 = C3*C1                                                                                                                    
      C6 = FLOAT(4*N*N-64)/63.0E0                                                                                                   
      C7 = C4+C6                                                                                                                    
      C8 = C5+C6*C2                                                                                                                 
      DO 050 I=1,N                                                                                                                  
          XI     = FLOAT(2*I-N-1)/2.E0                                                                                              
          Z0     = UN                                                                                                               
          Z1     = Z0 * XI                                                                                                          
          Z2     = Z1 * XI                                                                                                          
          Z3     = Z2 * XI                                                                                                          
          Z4     = Z3 * XI                                                                                                          
          Z5     = Z4 * XI                                                                                                          
          P(I,1) = UN                                                                                                               
          P(I,2) = Z1                                                                                                               
          P(I,3) = Z2 - C1                                                                                                          
          P(I,4) = Z3 - Z1 * C2                                                                                                     
          P(I,5) = Z4 - Z2 * C4 + C5                                                                                                
          P(I,6) = Z5 - Z3 * C7 + Z1 * C8                                                                                           
  050 CONTINUE                                                                                                                      
                                                                                                                                    
C     * GET VARIANCE OF THE PI'S.                                                                                                   
                                                                                                                                    
      DR(1) = FLOAT(N)                                                                                                              
      DO 075 L=2,K+1                                                                                                                
          DR(L) = 0.0E0                                                                                                             
          DO 075 I=1,N                                                                                                              
              DR(L) = DR(L)+P(I,L)*P(I,L)                                                                                           
  075 CONTINUE                                                                                                                      
                                                                                                                                    
C     * GET DATA.                                                                                                                   
                                                                                                                                    
      NR = 0                                                                                                                        
  100 CALL GETFLD2 (1,F,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                                   
      IF (.NOT.OK)   THEN                                                                                                           
          WRITE(6,'("0 FIT PROCESSED",I5,"RECORDS.")') NR                                                                           
          CALL                                     XIT('FIT',0)                                                                     
      ENDIF                                                                                                                         
                                                                                                                                    
      NWDS = IBUF(5)*IBUF(6)                                                                                                        
                                                                                                                                    
C     * REMOVE MEAN.                                                                                                                
                                                                                                                                    
      F0 = 0.0E0                                                                                                                    
      DO 102 I=1,NWDS                                                                                                               
          F0 = F0 +F(I)                                                                                                             
  102 CONTINUE                                                                                                                      
      F0 = F0/FLOAT(NWDS)                                                                                                           
                                                                                                                                    
      DO 103 I=1,NWDS                                                                                                               
          F(I) = F(I)-F0                                                                                                            
  103 CONTINUE                                                                                                                      
                                                                                                                                    
C     * BUNCH DATA INTO N AVERAGES.                                                                                                 
                                                                                                                                    
      IF (NWDS.LT.N) CALL                          XIT('FIT',-2)                                                                    
      NW  = NWDS/N                                                                                                                  
      XNW = FLOAT(NW)                                                                                                               
      DO 110 I=1,N                                                                                                                  
          X(I) = 0.0E0                                                                                                              
          N0   = (I-1)*NW                                                                                                           
          DO 105 J=1,NW                                                                                                             
              X(I) = X(I)+F(J+N0)                                                                                                   
  105     CONTINUE                                                                                                                  
          X(I) = X(I)/XNW                                                                                                           
  110 CONTINUE                                                                                                                      
                                                                                                                                    
C     * GET THE EXPANSION COEFFICIENTS.                                                                                             
                                                                                                                                    
      DO 120 L=1,K+1                                                                                                                
          A(L) = 0.0E0                                                                                                              
          DO 115 I=1,N                                                                                                              
              A(L) = A(L) + X(I)*P(I,L)                                                                                             
  115     CONTINUE                                                                                                                  
          A(L) = A(L)/DR(L)                                                                                                         
  120 CONTINUE                                                                                                                      
                                                                                                                                    
C     * EVALUATE S.                                                                                                                 
                                                                                                                                    
      Y2 = 0.0E0                                                                                                                    
      DO 125 I=1,N                                                                                                                  
          Y2 = Y2 + X(I)*X(I)                                                                                                       
  125 CONTINUE                                                                                                                      
                                                                                                                                    
      IF (Y2.EQ.0.0E0)  THEN                                                                                                        
          NR = NR + 1                                                                                                               
          WRITE(6,'("1 RECORD NO ",I5," IS NULL.")') NR                                                                             
          GOTO 100                                                                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
      S = Y2                                                                                                                        
      DO 130 L=1,K+1                                                                                                                
          S = S - A(L)*A(L)*DR(L)                                                                                                   
  130 CONTINUE                                                                                                                      
      S = SQRT( S/FLOAT(N-K-1))                                                                                                     
                                                                                                                                    
C     * EVALUATE TEST STATISTIC.                                                                                                    
                                                                                                                                    
      T(1) = 0.0E0                                                                                                                  
      DO 135 L=2,K+1                                                                                                                
          T(L) = A(L)*SQRT(DR(L))/S                                                                                                 
  135 CONTINUE                                                                                                                      
                                                                                                                                    
C     * WRITE OUT RESULTS.                                                                                                          
                                                                                                                                    
      NR = NR + 1                                                                                                                   
      WRITE(6,'("1 FOR RECORD ",I5)') NR                                                                                            
      WRITE(6,'("0 BUNCHED X "/(6(I5,E15.6)))') (I,X(I),I=1,N)                                                                      
      WRITE(6,'("0      COEFFICIENTS AND TEST STATISTIC VALUES"/                                                                    
     1         ("0",I5.1,2E15.5))') (I-1,A(I),T(I),I=1,K+1)                                                                         
      WRITE(6,'("0 STANDART DEVIATION IS   ",E15.5)') S                                                                             
      WRITE(6,'("0 MEAN OF INPUT RECORD IS ",E15.5)') F0                                                                            
                                                                                                                                    
      IF (N.LE.500)   THEN                                                                                                          
          D  = 1.E0                                                                                                                 
          NT = N                                                                                                                    
      ELSE                                                                                                                          
          D  = (FLOAT(N)-1.E0)/499.E0                                                                                               
          NT = 500                                                                                                                  
      ENDIF                                                                                                                         
                                                                                                                                    
      DX(1) = 1.E0                                                                                                                  
      DO 140 I=1,NT-1                                                                                                               
          DX(I+1)   = DX(I)+D                                                                                                       
          IX        = INT(DX(I+1))                                                                                                  
          X(I+1)    = X(IX)                                                                                                         
          X(I+NT+1) = A(1)                                                                                                          
          DO 140 L=2,K+1                                                                                                            
              X(I+NT+1) = X(I+NT+1) + A(L)*P(IX,L)                                                                                  
  140 CONTINUE                                                                                                                      
                                                                                                                                    
      XNT1 = A(1)                                                                                                                   
      DO 145 L=2,K+1                                                                                                                
          XNT1 = XNT1 + A(L)*P(1,L)                                                                                                 
  145 CONTINUE                                                                                                                      
      X(NT+1) = XNT1                                                                                                                
                                                                                                                                    
      R(1) = 1.E0                                                                                                                   
      R(2) = DX(NT)                                                                                                                 
      R(3) = 0.E0                                                                                                                   
      R(4) = 0.E0                                                                                                                   
      CALL USPLO( DX, X, NT, NT, 2, 1,                                                                                              
     1     'INPUT FIELD (MINUS MEAN) AND APPROXIMATION TO IT.',49,                                                                  
     2     'POSITION IN SERIES',18,                                                                                                 
     3     'INPUT AND APPROXIMATION',23,                                                                                            
     4     R, 'IA',1,IER)                                                                                                           
                                                                                                                                    
      GOTO 100                                                                                                                      
                                                                                                                                    
C     * PREMATURE EOF ON INPUT.                                                                                                     
                                                                                                                                    
  900 CALL                                         XIT('FIT',-3)                                                                    
                                                                                                                                    
C     * PREMATURE EOF ON UNIT 1.                                                                                                    
                                                                                                                                    
  999 CALL                                         XIT('FIT',-4)                                                                    
                                                                                                                                    
C---------------------------------------------------------------------                                                              
      END                                                                                                                           
