      SUBROUTINE PREP (A,T,NI,NF,DELT,NPOLY,NORM) 
C     
C     THIS SUBROUTINE IS TAKEN FROM OCTPUS AND THE ADDING ZEROS   
C     AND COSINE WINDOW HAVE BEEN REMOVED WHILE SPLINE FITTING HAS BEEN 
C     ADDED AND THE NORMALIZATION SCHEME HAS CHANGED  
C     
C     THIS SUBROUTINE PREPARES THE DATA FOR ANALYSIS BY REMOVING TRENDS,
C     AND NORMALIZING.  
C     ***** INPUT DATA *****  
C     NPOLY AND NORM IN A FREE FORMAT, WHERE    
C     NPOLY = DEGREE OF POLYNOMIAL USED TO REMOVE TRENDS IF NPOLY >0;   
C           = NUMBER OF CHANNELS TO AVERAGE TOGETHER TO PRODUCE POINTS  
C             FOR THE SPLINE FITTING IF NPOLY<0.
C     NORM = 0, ONLY TRENDS ARE REMOVED. NO NORMALIZATION.  
C     NORM = 888 ONLY TRENDS ARE REMOVED. NO NORMALIZATION. COSINE TAPER
C          = 1, DATA NORMALIZED BY THE AVERAGE OF THE DATA. 
C          = 2, DATA NORMALIZED BY THE FITTED CURVE (SPLINE OR POLYNOMIAL)    
C          < 0, AS FOR DATA > 0 EXCEPT COSINE BELL TAPER APPLIED. 
C
      DIMENSION A(1),ALPHA(6),XPOW(6),YY(2000),XX(2000)
      DOUBLE PRECISION T(*)
      ITAPER=0    
      NORMIN=NORM 
      IF(NORM.LT.0.OR.NORM.EQ.888)THEN    
        IF(NORM.LT.0)NORM=-NORM     
        ITAPER=1  
      END IF
      NO=NF-NI+1  
      PRINT 2,NO,NI,NF  
    2 FORMAT(12X,9HPREPQ OF ,I6,18H POINTS IN Q FROM ,I6,4H TO ,I6)     
      IF(NORM.NE.999)PRINT 3,NPOLY,NORM   
    3 FORMAT(12X,8HNPOLY = ,I5,/,12X,7HNORM = ,I8)    
      AVERAG =0.0 
      DO 4 I=NI,NF
      AVERAG=AVERAG +A(I)     
    4 CONTINUE    
      AVERAG=AVERAG /NO 
      IF(NORM.NE.999)PRINT 1000,AVERAG    
 1000 FORMAT(12X,10HAVERAGE = ,E15.8)     
      IF(NPOLY.GT.0)GO TO 15  
C     
C     IF NPOLY=0 THEN JUST SUBTRACT THE MEAN    
C     
      IF(NPOLY.EQ.0)THEN
        DO 5 J=NI,NF    
          A(J)=A(J)-AVERAG    
          IF(NORM.GT.0)A(J)=A(J)/AVERAG   
    5   CONTINUE  
        GO TO 50  
      END IF
C     
C     SPLINE FITTING SECTION  
C     
      READ *,NPTS,NSLIDE     
C
C     THIS CALCULATES THE KNOTS OF THE SPLINE BY AVERAGING
C     TOGETHER NSLIDE POINTS FOR THE FIRST KNOT. THE SECOND
C     KNOT COMES FROM THE FIRST NPTS POINTS OF THE CURVE,
C     AVERAGED TOGETHER. THE THIRD ONE COMES FROM POINTS 
C     NSLIDE+1 TO NSLIDE+NPTS ETC. AND THE LAST POINT IS
C     MADE FROM THE LAST NSLIDE POINTS.
C
      NEWNO=0
      LAST=0
    6 NEWNO=NEWNO+1
      IF(NEWNO.GT.2000)THEN    
        PRINT *,'REDIMENSION XX,YY ARRAYS IN PREP LARGER'   
        STOP 'STORAGE OVERFLOW IN PREP'   
      END IF
      XX(NEWNO)=0.0
      YY(NEWNO)=0.0
      IF(NEWNO.EQ.1)THEN
        IS=NI
        IE=MIN(NF,NI+NSLIDE-1)
      ELSE IF(NEWNO.EQ.2)THEN
        IS=NI
        IE=NI+NPTS-1
      ELSE
        IS=IS+NSLIDE
        IE=IE+NSLIDE
      END IF
      IF(IE.GT.NF)THEN
        IS=NF-NSLIDE+1
        IE=NF
        LAST=1
      END IF
      DO 7 J=IS,IE
        XX(NEWNO)=XX(NEWNO)+T(J)
        YY(NEWNO)=YY(NEWNO)+A(J)
    7 CONTINUE
      XX(NEWNO)=XX(NEWNO)/(IE-IS+1)
      YY(NEWNO)=YY(NEWNO)/(IE-IS+1)
      IF(LAST.EQ.0)GO TO 6
      PRINT 8,(XX(JJJ),JJJ=1,NEWNO)
    8 FORMAT(12X,27HTHE SPLINE FITTING ARRAYS :/
     *       (14X,10F11.1))   
      PRINT 9,(YY(JJJ),JJJ=1,NEWNO)
    9 FORMAT(14X,10F11.1)
      DO 10 J=NI,NF     
        TT=T(J+1-NI)
        call intep(tt,yt,xx,yy,newno,ier)
        if (ier.eq.2) then
           write(6,12) tt
   12      format(1x,'### WARNING: Extrapolated value at ',f8.1,
     *            ' set to spline end values')
        endif
        A(J)=A(J)-yt
        IF(NORM.EQ.1)A(J)=A(J)/AVERAG     
        IF(NORM.EQ.2)A(J)=A(J)/yt
   10 CONTINUE    
      GO TO 50    
C     
C     POLYNOMIAL FITTING SECTION    
C     
   15 CALL FITTER(A,T,NI,NF,NO,NPOLY,ALPHA,NORM)  
      NN=NPOLY+1  
      DO 40 I=1,NO
        XPOW(1)=1.0 
        DO 20 J=2,6 
          XPOW(J)=XPOW(J-1)*(T(I+NI-1)-T(NI))          
   20   CONTINUE    
        DA=0.0
        DO 30 J=1,NN
          DA=DA+ALPHA(J)*XPOW(J)  
   30   CONTINUE    
        A(NI+I-1)=A(NI+I-1)-DA  
   35   IF(NORM.EQ.1)A(NI+I-1)=A(NI+I-1)/AVERAG   
        IF(NORM.EQ.2)A(NI+I-1)=A(NI+I-1)/DA 
   40 CONTINUE    
   50 VAR=0.0     
      DO 80 I=NI,NF     
        VAR=VAR+A(I)*A(I) 
   80 CONTINUE    
      VAR=VAR/NO  
      IF(NORM.NE.999)PRINT 90,VAR   
   90 FORMAT(12X,11HVARIANCE = ,E15.8)    
  120 CONTINUE    
      NORM=NORMIN 
      IF(ITAPER.EQ.0)RETURN   
      CALL COSWIN (A,NI,NF,1) 
C     THIS CORRECTS FOR THE EFFECTS OF THE 10 PERCENT COSINE BELL 
C         WINDOW ON THE VARIANCE.   
      FAC=SQRT(1.0/0.875)     
      DO 130 I=NI,NF    
      A(I)=A(I)*FAC     
  130 CONTINUE    
      RETURN
      END   
      SUBROUTINE FITTER (A,T,NI,NF,NO,N,ALPHA,NORM)     
C     
C     THIS SUBROUTINE FITS A POLYNOMIAL UP TO FIFTH DEGREE BY LEAST     
C     SQUARES AND RETURNS ITS COEFFICIENTS IN ALPHA.  
C     A = REAL DATA ARRAY WITH BAD POINTS SET TO LE -1.0E9. 
C     NO = NUMBER OF DATA POINTS IN A     
C     N = DEGREE OF POLYNOMIAL. IE. FOF LINEAR FIT, SET N=1.
C     
      DIMENSION A(1),XX(6,6),ALPHA(6),XY(6),XPOW(11)  
      DOUBLE PRECISION T(1)
      NN=N+1
      DO 20 I=1,6 
        DO 10 J=1,6 
          XX(I,J)=0.0 
   10   CONTINUE    
        ALPHA(I)=0.0
        XY(I)=0.0   
   20 CONTINUE    

      IF(NN.EQ.6)GO TO 45     
      NNN=NN+1    
      DO 40 I=NNN,6     
        XX(I,I)=1.0 
   40 CONTINUE    

   45 DO 80 K=1,NO
        XPOW(1)=1.0 
        DO 50 L=2,2*NN-1
          XPOW(L)=XPOW(L-1)*(T(K+NI-1)-T(NI))     
   50   CONTINUE    
        DO 70 I=1,NN
          DO 60 J=1,NN
            XX(I,J)=XX(I,J)+XPOW(I+J-1)   
   60     CONTINUE    
          XY(I)=XY(I)+XPOW(I)*A(K+NI-1) 
   70   CONTINUE    
   80 CONTINUE    
      CALL SYMINV(XX,6) 
      DO 100 I=1,NN     
        DO 90 J=1,NN
          ALPHA(I)=ALPHA(I)+XX(I,J)*XY(J)     
   90   CONTINUE    
  100 CONTINUE    

      IF(NORM.NE.999)PRINT 110
  110 FORMAT(14X,23HPOLYNOMIAL COEFFICIENTS)    
      DO 130 I=1,NN     
        IF(NORM.NE.999)PRINT 120,I,ALPHA(I) 
  120   FORMAT(16X,6HALPHA(,I1,4H) = ,E15.8)
  130 CONTINUE    
      RETURN
      END   
      SUBROUTINE SYMINV (A,N) 
      DIMENSION A(N,N)  
      DIMENSION P(20),Q(20),R(20)   
   13 FORMAT(15HOSYMINV FAILED )    
      ZERO=0.0    
      ONE=1.0     
      DO 1 M=1,N  
    1 R(M)=ONE    
      DO 2 M=1,N  
      BIG=ZERO    
      DO 3 L=1,N  
      AB=ABS (A(L,L))   
      IF(AB-BIG)3,3,4   
    4 IF(R(L))14,3,14   
   14 BIG=AB
      K=L   
    3 CONTINUE    
      IF(BIG)6,5,6
    5 PRINT 13    
      RETURN
    6 R(K)=ZERO   
      Q(K)=ONE/A(K,K)   
      P(K)=ONE    
      A(K,K)=ZERO 
      KM1=K-1     
      IF(KM1.EQ.0)GO TO 16    
   18 DO 7 L=1,KM1
      P(L)=A(L,K) 
      IF(R(L))9,8,9     
    8 Q(L)=A(L,K)*Q(K)  
      GO TO 7     
    9 Q(L)=-A(L,K)*Q(K) 
    7 A(L,K)=ZERO 
   16 CONTINUE    
      KP1=K+1     
      IF(KP1.GT.N)GO TO 17    
   19 DO 15 L=KP1,N     
      IF(R(L))12,11,12  
   12 P(L)=A(K,L) 
      GO TO 10    
   11 P(L)=-A(K,L)
   10 Q(L)=(-A(K,L))*Q(K)     
   15 A(K,L)=ZERO 
   17 CONTINUE    
      DO 2 L=1,N  
      DO 2 K=L,N  
    2 A(L,K)=A(L,K)+P(L)*Q(K) 
      M=N+1 
      L=N   
      DO 27 K=2,N 
      M=M-1 
      L=L-1 
      DO 27 J=1,L 
   27 A(M,J)=A(J,M)     
      RETURN
      END   
