#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 5 ]
then
  echo " "
  echo "            *****     FOUFIL     *****      "
  echo " "
  echo "      The program filters the data"
  echo " "
  echo "Syntax: foufil <file1> <file2> <int1> <int2> <int3> <int4>"
  echo " "
  echo "      <file1> - input file 1"
  echo "      <file2> - input file 2"
  echo "      <int1> - window's smoothed lower limit"
  echo "      <int2> - window's lower limit"
  echo "      <int3> - window's upper limit"
  echo "      <int4> - window's smoothed upper limit"
  echo "  "
  exit 0
else
 echo "file1="$1
 echo "file2="$2
 echo "Window's smoothed lower limit= "$3
 echo "Window's lower limit="$4
 echo "Window's upper limit="$5
 echo "Window's smoothed upper limit="$6
 cp $1 fort.1
fi
cat > foufil.f <<EOF
      PROGRAM PFOUFIL
      IMPLICIT REAL (A-H,O-Z)

      PRINT *,' '
      PRINT *,'    *****     FOUFIL     *****      '
      PRINT *,'   '
      NCHAN=1
      CALL PARLES (NCHAN,NDIM,NSET)
      CALL FOUFIL (NDIM,NSET) 
      STOP 
      END
C******************************************************************
      SUBROUTINE PARLES (NCHAN,NDIM,NSET)
      IMPLICIT REAL (A-H,O-Z)
C
      INTEGER I1,I2,I3

      REWIND NCHAN
      READ(NCHAN,END=10) ITIME,INAME,ILEV,ISIZE
      REWIND NCHAN

      NSET=0
      NDIM=0  
    1 READ (NCHAN,END=100) I1,I2,I3,I4
      NDIM=MAX0(NDIM,I4)
      READ (NCHAN,END=101) 
      NSET=NSET+1
      GOTO1
      
  101 PRINT *,' MORE HEADER THEN DATA '     
  100 PRINT 1000,NCHAN,NDIM,NSET
      REWIND NCHAN
      RETURN    
   10 PRINT *,'INPUT FILE ON TAPE',NCHAN,' IS EMPTY'
      STOP 'MISTAKE'
 1000 FORMAT (1X,' TAPE=',I2,',  NDIM=',I5,',  NSET=',I5)
      END

C*****************************************************************
      SUBROUTINE FOUFIL (IE,NTS)
      IMPLICIT REAL (A-H,O-Z)
C
      PARAMETER (IE77=15000,MAXF=100,UNDEF=0.9E+10,UNDEFE=0.99*UNDEF,
     &           NTS77=15000)
      DIMENSION ID(4)

C      DIMENSION X(IE), C(IE,MAXF), S(IE,MAXF) ,ID(4), FILT(NTS),
C     &          DC(IE,MAXF),DS(IE,MAXF),XM(IE), COUNT(IE)
C      POINTER (NPX,X),(NPC,C),(NPS,S),(NPFILT,FILT),(NPDC,DC),
C     &       (NPDS,DS),(NPXM,XM),(NPCOUNT,COUNT)
C      ALLOCATE (NPX,IE)
C      ALLOCATE (NPC,IE*MAXF)
C      ALLOCATE (NPS,IE*MAXF)
C      ALLOCATE (NPFILT,NTS)
C      ALLOCATE (NPDC,IE*MAXF)
C      ALLOCATE (NPDS,IE*MAXF)
C      ALLOCATE (NPXM,IE)
C      ALLOCATE (NPCOUNT,IE)

      REAL, DIMENSION(:), ALLOCATABLE :: X
      REAL, DIMENSION(:,:), ALLOCATABLE :: C
      REAL, DIMENSION(:,:), ALLOCATABLE :: S
      REAL, DIMENSION(:), ALLOCATABLE :: FILT
      REAL, DIMENSION(:,:), ALLOCATABLE :: DC
      REAL, DIMENSION(:,:), ALLOCATABLE :: DS
      REAL, DIMENSION(:), ALLOCATABLE :: XM
      REAL, DIMENSION(:), ALLOCATABLE :: COUNT

      ALLOCATE (X(IE),C(IE,MAXF),S(IE,MAXF),FILT(NTS))
      ALLOCATE (DC(IE,MAXF),DS(IE,MAXF),XM(IE),COUNT(IE))

C AND ACTIVATE THE FOLLOWING FOUR LINES
C      PARAMETER (IE77=3000,MAXF=5,UNDEF=UVALUE,UNDEFE=0.99*UNDEF,
C     &           NTS=3000)
C      DIMENSION X(IE77), C(IE77,MAXF), S(IE77,MAXF) ,ID(4), FILT(NTS),
C     &          DC(IE77,MAXF),DS(IE77,MAXF),XM(IE77), COUNT(IE77)
C NOW YOU HAVE DONE ALL CHANGES FOR STANDARD FORTRAN.
      DATA NOOUT/3/, N/0/                  
      PI = 4.*ATAN(1.)
      TWOPI=2.*PI

      DO 100 I=1,IE
        COUNT(I)=0
        XM(I)=0.0
  100 CONTINUE

      PRINT *,'ENTER FILTER PERIODICITIES: PMIN,P2,P1,PMAX'
      READ(5,*,END=9925) PMIN,P2,P1,PMAX
      PRINT*,' PMIN=',PMIN,', P2=',P2,', P1=',P1,', PMAX=',PMAX
      IF(PMIN.EQ.0 .AND. PMAX.EQ.0) NOOUT = 2

C
C     * check the input file TAPE1 and copy it to TAPE3
C
      REWIND 1
      REWIND 3
      READ(1,END=9915) ITIME,INAME,ILEV,ISIZE
      PRINT *, 'first record:', ITIME, INAME, ILEV, ISIZE
      IF(ISIZE.GT.IE77) GOTO 9920
      REWIND 1
 1100   READ(1,END=1000) ID
        IF(INAME.NE.ID(2) .OR. ILEV.NE.ID(3) .OR. ISIZE.NE.ID(4))   
     &  GOTO 9935
        READ(1) (X(J),J=1,ISIZE)
        N = N + 1
        DO 1194 I=1,ISIZE
          IF(X(I).GT.UNDEFE) GOTO 1194
            COUNT(I) = COUNT(I)+1
            XM(I) = XM(I)+X(I)
 1194     CONTINUE
        WRITE(NOOUT) ID
        WRITE(NOOUT) (X(J),J=1,ISIZE)
        GOTO 1100
C
C   e-o-f FOUND
c
 1000 IF(NOOUT.EQ.2) STOP
      IF(N.GT.NTS) GOTO 9916
      DO 1195 I=1,ISIZE
        IF(COUNT(I).NE.0) XM(I)=XM(I)/COUNT(I)
 1195   IF(COUNT(I).EQ.0) XM(I)=UNDEF

      NF=N/2
      IF(PMAX.LT.1.) PMAX=FLOAT(N)
      IF(P1.LT.1.) P1=FLOAT(N)
      IF(P2.LT.1.) P2=2.
      IF(PMIN.LT.1.) PMIN=2.
      IF(.NOT.((PMAX.GE.P1).AND.(P1.GE.P2).AND.
     *         (P2.GE.PMIN))) THEN
               PRINT *,' ** MUST HAVE PMAX>=P1>=P2>=PMIN ** '
               STOP
               ENDIF 
      LMIN=N/PMIN
      L2=N/P2
      L1=N/P1
      LMAX=N/PMAX
      IF(L1.LT.1) L1=1
      IF(LMAX.LT.1) LMAX=1
      IF(L2.GT.NF) L2=NF
      IF(LMIN.GT.NF) LMIN=NF
      PRINT *, ' FILTER: ',PMIN,P2,P1,PMAX, ' delta T'
C
C     calculating the No. of frequencies to be treated
C     NFW1+NFW2: No. of frequencies to be substracted
C     NFB: No. of frequencies to bo constracted
C
      
      FMIN=1./PMIN
      F2=1./P2
      F1=1./P1
      FMAX=1./PMAX
      NFW1=NF-L2+1
      NFW2=L1-1+1
      IF(L1.EQ.1) NFW2=0
      IF(L2.EQ.NF) NFW1=0
      NFW=NFW1 + NFW2     
      NFB=LMIN-LMAX+1                     
       

      IF(NFW.LE.NFB) THEN
C
C*    1. MAIN LOOP (SUBSTRACTION)

C
C     * calculating the filter weights *FILT* for given filter parameter
C

      DO 1301 I=1,NF
      F=FLOAT(I)/FLOAT(N)
      IF(I.EQ.1) THEN
          P=FLOAT(N)
          ELSE
          P=1./F
          ENDIF
      FILT(I)=0.
      IF(P.LT.PMIN .OR. P.GT.PMAX) FILT(I)=1.
      IF(P.GE.PMIN .AND. P.LT.P2)
     1    FILT(I) = (1.-COS(PI*(F-F2)/(FMIN-F2))) / 2.
      IF(P.LE.PMAX .AND. P.GT.P1)
     1    FILT(I) = (1.+COS(PI*(F-FMAX)/(F1-FMAX))) / 2.
 1301 CONTINUE

      IF(NFW1.EQ.0) GOTO 9983

      PRINT *,'subtract periodicities at frequencies > ',
     1         L2, ' cycles per ',N,' DELTA T'
  
C     1.1 substracting periodicities at frequencies > L2    
  
      N1=1
      N2=MIN0(NFW1,N1+MAXF-1)
C
 8000 REWIND 3
      DO 2000 LL=N1,N2
        L=LL-N1+1
        DO 2000 J=1,ISIZE
          DC(J,L) = 0
          DS(J,L) = 0
          C(J,L)=0.0
 2000     S(J,L)=0.0
C
      DO 2100 K=1,N 
         READ(3) ID
         READ(3) (X(J),J=1,ISIZE)
         DO 2100 LL=N1,N2
           L=LL-N1+1
           ARG=TWOPI*(L2+LL-1)*K/N
           COSA=COS(ARG)
           SINA=SIN(ARG)
           CC = COSA**2
           SS = SINA**2
           DO 2100 J=1,ISIZE
             IF(X(J).GT.UNDEFE) GOTO 2100
             DS(J,L) = DS(J,L) + SS
             DC(J,L) = DC(J,L) + CC
             C(J,L)=C(J,L) + COSA*(X(J)-XM(J))
             S(J,L)=S(J,L) + SINA*(X(J)-XM(J))
 2100      CONTINUE

C
C     * compute the trigonometric contributions and subtract it
C     * from TAPE3. Write result on TAPE2.
C
      REWIND 3
      REWIND 2
      DO 3000 LL=N1,N2
        L=LL-N1+1
        DO 3000 J=1,ISIZE
          IF(DS(J,L).EQ.0) S(J,L)=UNDEF
          IF(DC(J,L).EQ.0) C(J,L)=UNDEF
          IF(DS(J,L).NE.0) S(J,L) = S(J,L)/DS(J,L)
 3000     IF(DC(J,L).NE.0) C(J,L) = C(J,L)/DC(J,L)
C
      DO 3100 K = 1,N 
         READ(3) ID
         READ(3) (X(J),J=1,ISIZE)
         DO 3110 LL=N1,N2
           L=LL-N1+1
           ARG=TWOPI*(L2+LL-1)*K/N
           COSA=COS(ARG)
           SINA=SIN(ARG)
           DO 3110 J=1,ISIZE
               IF(C(J,L).LT.UNDEFE .AND. X(J).LT.UNDEFE)      
     *         X(J) = X(J)-C(J,L)*FILT(L2+LL-1)*COSA 
               IF(S(J,L).LT.UNDEFE .AND. X(J).LT.UNDEFE)     
     *         X(J) = X(J)-S(J,L)*FILT(L2+LL-1)*SINA
 3110      CONTINUE
         WRITE(2) ID
 3100    WRITE(2) (X(J),J=1,ISIZE)      
C
C  IF MORE FREQUENCIES ARE CONSIDERED, NEW LOOP OF CALCULATIONS
         N1=N1+MAXF 
         IF(N1.GT.NFW1) GOTO 9983  
         N2=MIN0(NFW1,N1+MAXF-1)
         REWIND 2 
         REWIND 3
         DO 4000 K = 1,N 
           READ(2) ID
           READ(2) (X(J),J=1,ISIZE)
           WRITE(3) ID
 4000      WRITE(3) (X(J),J=1,ISIZE)      
         GOTO 8000
C
C*    1.2 substracting periodicities at frequences < L1
C

 9983 CONTINUE

      IF(NFW2.EQ.0) GOTO 9900
      PRINT *,'subtract periodicities at frequencies < ',
     1         L1, ' cycles per ',N,' DELTA T'

      N1=1
      N2=MIN0(NFW2,N1+MAXF-1)
C
 8300 REWIND 3
      DO 2300 LL=N1,N2
        L=LL-N1+1
        DO 2300 J=1,ISIZE
          DC(J,L) = 0
          DS(J,L) = 0
          C(J,L)=0.0
 2300     S(J,L)=0.0
C
      DO 2310 K=1,N 
         READ(3) ID
         READ(3) (X(J),J=1,ISIZE)
         DO 2310 LL=N1,N2
           L=LL-N1+1
           ARG=TWOPI*(1+LL-1)*K/N
           COSA=COS(ARG)
           SINA=SIN(ARG)
           CC = COSA**2
           SS = SINA**2
           DO 2310 J=1,ISIZE
             IF(X(J).GT.UNDEFE) GOTO 2310
             DS(J,L) = DS(J,L) + SS
             DC(J,L) = DC(J,L) + CC
             C(J,L)=C(J,L) + COSA*(X(J)-XM(J))
             S(J,L)=S(J,L) + SINA*(X(J)-XM(J))
 2310      CONTINUE
C
C     * compute the trigonometric contributions and subtract it
C     * from TAPE3. Write result on TAPE2.
C
      REWIND 3
      REWIND 2
      DO 3300 LL=N1,N2
        L=LL-N1+1
        DO 3300 J=1,ISIZE
          IF(DS(J,L).EQ.0) S(J,L)=UNDEF
          IF(DC(J,L).EQ.0) C(J,L)=UNDEF
          IF(DS(J,L).NE.0) S(J,L) = S(J,L)/DS(J,L)
 3300     IF(DC(J,L).NE.0) C(J,L) = C(J,L)/DC(J,L)
C
      DO 3331 K = 1,N 
         READ(3) ID
         READ(3) (X(J),J=1,ISIZE)
         DO 3310 LL=N1,N2
           L=LL-N1+1
           ARG=TWOPI*(1+LL-1)*K/N
           COSA=COS(ARG)
           SINA=SIN(ARG)
           DO 3310 J=1,ISIZE
             IF(C(J,L).LT.UNDEFE .AND. X(J).LT.UNDEFE)                  
     *       X(J) = X(J)-C(J,L)*FILT(1+LL-1)*COSA 
             IF(S(J,L).LT.UNDEFE .AND. X(J).LT.UNDEFE)                 
     *       X(J) = X(J)-S(J,L)*FILT(1+LL-1)*SINA
 3310        CONTINUE                                                  
         WRITE(2) ID
 3331    WRITE(2) (X(J),J=1,ISIZE)      
C
C  IF MORE FREQUENCIES ARE CONSIDERED, NEW LOOP OF CALCULATIONS
         N1=N1+MAXF 
         IF(N1.GT.NFW2) GOTO 9900
         N2=MIN0(NFW2,N1+MAXF-1)
         REWIND 2 
         REWIND 3
         DO 4300 K = 1,N 
           READ(2) ID
           READ(2) (X(J),J=1,ISIZE)
           WRITE(3) ID
 4300      WRITE(3) (X(J),J=1,ISIZE)      
         GOTO 8300


      ELSE      

C
C*    2. MAIN LOOP
C     2.1 constraucting  

C
C     * calculating the filter weights *FILT* for given filter parameter
C

      PRINT *,'reconstruct periodicities at frequencies ',LMAX,        
     *        ' to ', LMIN, ' cycles per ',N, ' delta T'
      DO 1331 I=1,NF
      F=FLOAT(I)/FLOAT(N)
      IF(I.EQ.1) THEN
          P=FLOAT(N)
          ELSE
          P=1./F
          ENDIF
      FILT(I)=1.
      IF(P.LT.PMIN .OR. P.GT.PMAX) FILT(I)=0.
      IF(P.GE.PMIN .AND. P.LT.P2)
     1    FILT(I) = (1.+COS(PI*(F-F2)/(FMIN-F2))) / 2.
      IF(P.LE.PMAX .AND. P.GT.P1)
     1    FILT(I) = (1.-COS(PI*(F-FMAX)/(F1-FMAX))) / 2.
 1331 CONTINUE
     
      REWIND 3
      REWIND 4
      DO 3611 K=1,N
        READ(3) ID
        READ(3)
        WRITE(4) ID
 3611   WRITE(4) (XM(J),J=1,ISIZE)

      N1=1
      N2=MIN0(NFB,N1+MAXF-1)
C
C     * compute the trigonometric contributions
C
 8500 REWIND 3
      DO 2500 LL=N1,N2
        L=LL-N1+1
        DO 2500 J=1,ISIZE
          DC(J,L) = 0.
          DS(J,L) = 0.
          C(J,L) = 0.0
 2500     S(J,L) = 0.0
C
      DO 2510 K=1,N 
         READ(3) ID
         READ(3) (X(J),J=1,ISIZE)
         DO 2510 LL=N1,N2
           L=LL-N1+1
           ARG=TWOPI*(LMAX+LL-1)*K/N
           COSA=COS(ARG)
           SINA=SIN(ARG)
           CC = COSA**2
           SS = SINA**2
           DO 2510 J=1,ISIZE
             IF(X(J).GT.UNDEFE) GOTO 2510
             DS(J,L) = DS(J,L) + SS
             DC(J,L) = DC(J,L) + CC
             C(J,L)=C(J,L) + COSA*(X(J)-XM(J))
             S(J,L)=S(J,L) + SINA*(X(J)-XM(J))
 2510      CONTINUE
C
      DO 3510 LL=N1,N2
        L=LL-N1+1
        DO 3510 J=1,ISIZE
          IF(XM(J).GT.UNDEFE) THEN
            S(J,L) = UNDEF
            C(J,L) = UNDEF
          ELSE 
            S(J,L) = S(J,L)/DS(J,L)
            C(J,L) = C(J,L)/DC(J,L)
          ENDIF
 3510     CONTINUE
C
C     * construct the filtered data 
C     * Write result on TAPE2.
C
      REWIND 4
      REWIND 2
      DO 3520 K = 1,N 
         READ(4) ID
         READ(4) (X(J),J=1,ISIZE)
         DO 3530 LL=N1,N2
           L=LL-N1+1
           ARG=TWOPI*(LMAX+LL-1)*K/N
           COSA=COS(ARG)
           SINA=SIN(ARG)
           DO 3530 J=1,ISIZE
             IF(XM(J).LT.UNDEFE) THEN
               IF((LMAX+LL-1).EQ.NF .AND. 2*NF.EQ.N) THEN
               X(J) = X(J) +
     &         C(J,L)*FILT(LMAX+LL-1)*COSA                                
               ELSE
               X(J) = X(J) +
     &         C(J,L)*FILT(LMAX+LL-1)*COSA + S(J,L)*FILT(LMAX+LL-1)*SINA
               ENDIF
             ENDIF
 3530    CONTINUE
         WRITE(2) ID
 3520    WRITE(2) (X(J),J=1,ISIZE)      
C
C  IF MORE FREQUENCIES ARE CONSIDERED, NEW LOOP OF CALCULATIONS
         N1=N1+MAXF 
         IF(N1.GT.NFB) GOTO 9900
         N2=MIN0(NFB,N1+MAXF-1)
         REWIND 2 
         REWIND 4
         DO 4771 K=1,N
           READ(2) ID
           READ(2) (X(J),J=1,ISIZE)
           WRITE(4) ID
 4771      WRITE(4) (X(J),J=1,ISIZE)

        GOTO 8500
C
      ENDIF

 9900 PRINT *,'REGULAR END!'
      STOP
C
 9915 PRINT *, 'TAPE1 FILE IS EMPTY'
         GOTO 9910
 9916 PRINT *, 'TIME SERIES TOO LONG: ', N, ' > ', NTS
         GOTO 9910 
 9920 PRINT *, 'ARRAY TOO LONG (IE = ',IE77,')'
         GOTO 9910
 9925 PRINT *, 'INPUT FILE IS EMPTY'
         GOTO 9910
 9935 PRINT *, 'DATA ARE MIXED'
         GOTO 9910
C
C        * ABNORMAL TERMINATION
 9910    PRINT *, "foufil: FATAL ERROR"
         STOP
      END
EOF
lf95 foufil.f -o foufil.x
foufil.x $1 $2 << M
$3
$4
$5
$6
M
cp fort.2 $2
rm fort*
rm foufil.x foufil.f
exit