      PROGRAM FILTER

      PARAMETER(NTOT=14610)

      integer n
      real tim(ntot),xx(ntot),yy(ntot)
      real        xf(ntot),yf(ntot),wf(ntot)
      real        xl(ntot),yl(ntot)
      real        xs(ntot),ys(ntot)
      real g1(ntot),g2(ntot)
      real axx,bxx,ayy,byy
 

      XHIG=40. 
      XLOW=15.

      open(10,file='../tor_nhe/rnhe.corr',
     *         form='formatted')
      open(20,file='../tor_nhe/rnhe.bp',
     *         form='formatted')

      pi=acos(-1.)

      II=10
      IO=II+10
 
      xmea=0.
      do n=1,NTOT
      read(II,*)tim(n),xx(n)
       xmea=xmea+xx(n)/float(ntot)
      enddo
      do n=1,NTOT
       xx(n)=xx(n)-xmea
      enddo

      call dtrend(xx,axx,bxx,ntot)

      LIS=XHIG/(tim(2)-tim(1))
      LAS=XLOW/(tim(2)-tim(1))

C  FONCTION DE TRANSFER

      DY=tim(2)-tim(1)

      IF(II.EQ.10)CALL TRANSFERT(LIS,LAS,NTOT,DY)

      CALL FILTRE(xx,g1,LIS,N)
      CALL FILTRE(xx,g2,LAS,N)

      do i=1,ntot
      xf(I)=g2(i)-g1(i)             
      write(IO,*)tim(i),xf(i)       
      enddo

      close(II)
      close(IO)

 1    CONTINUE

      stop
      end


      SUBROUTINE FILTRE(F,F1,LAS,N)
      PARAMETER(KOR=10,JOR=10,MOR=100)
      DIMENSION F(N),W(-MOR:mor),G(-MOR:mor),F1(N)
      PI=ACOS(-1.)
      FC=1./FLOAT(LAS)
      CALL KISER(G,MOR)
      DO 1 I=-MOR,MOR
      IF (I.EQ.0) THEN
      W(I)=2.*FC
      ELSE
      W(I)=SIN(2.*PI*FC*FLOAT(I))/(PI*FLOAT(I))*G(I)
      ENDIF
   1  CONTINUE

      DO 2 I=1,N
      F1(I)=0.
      AT=0.

      IF(I.LE.KOR) THEN
      L1=-KOR
      L2=I-1
      ENDIF

      IF((I.GE.KOR+1).AND.(I.LE.MOR)) THEN
      L1=-I+1
      L2=I-1
      ENDIF

      IF((I.GE.MOR+1).AND.(I.LE.N-MOR)) THEN
      L1=-MOR
      L2=MOR
      ENDIF

      IF((I.GE.N-MOR+1).AND.(I.LE.N-JOR)) THEN
      L1=-N+I
      L2=N-I
      ENDIF

      IF(I.GE.N-JOR+1) THEN
      L1=-N+I
      L2=JOR
      ENDIF

      DO 3 K=L1,L2
      F1(I)=F1(I)+W(K)*F(I-K)
      AT=AT+W(K)
   3  CONTINUE
      F1(I)=F1(I)/AT
   2  CONTINUE

      RETURN
      END

      SUBROUTINE KISER(W,MOR)
      PARAMETER(LOR=100)
      DIMENSION W(-MOR:MOR),CO(-LOR:LOR)

      A=40.

      IF(A.LE.21.) THEN
      ALPHA=0.
      ENDIF

      IF((A.LT.50.).AND.(A.GT.21.)) THEN
      ALPHA=0.5842*(A-21.)**0.4+0.07886*(A-21.)
      ENDIF

      IF(A.GE.50.) THEN
      ALPHA=0.1102*(A-8.7)
      ENDIF

      DO 2 I=-MOR,MOR
      CO(I)=ALPHA*SQRT(1.-(FLOAT(I)/FLOAT(MOR))**2)

      N=0
      AS=1.
      AU=1.
      AS1=1.
      AU1=1.

      DO 10 K=1,300
      N=N+1

      AU=AU*((CO(I)/2.)/FLOAT(N))**2
      AS=AS+AU

      AU1=AU1*((ALPHA/2.)/FLOAT(N))**2
      AS1=AS1+AU1
  10  CONTINUE
      W(I)=AS/AS1
   2  CONTINUE
      RETURN
      END


       subroutine dtrend(x,ax,bx,n)
       real xxm,iim,xim,xm,im
       real x(n),det,ax,bx

       xxm=0.
       iim=0.
       xim=0.
       xm=0.
       im=0.
       do i=1,n
         xxm=xxm+x(i)*x(i)/float(n)
         iim=iim+float(i*i)/float(n)
         xim=xim+x(i)*float(i)/float(n)
         xm=xm+x(i)/float(n)
         im=im+float(i)/float(n)
       enddo

       det=iim-im**2
       ax=(xim-xm*im)/det
       bx=(iim*xm-xim*im)/det

       do i=1,n
          x(i)=x(i)-ax*float(i)-bx
       enddo

       return
       end

      SUBROUTINE TRANSFERT(LIS,LAS,N,DY)
      PARAMETER(KOR=10,JOR=10,MOR=100)
      DIMENSION W1(-MOR:mor),G(-MOR:MOR),W2(-MOR:MOR)
      PI=ACOS(-1.)

      FC=1./FLOAT(LIS)
      CALL KISER(G,MOR)
      AT=0.
      DO 1 I=-MOR,MOR
      IF (I.EQ.0) THEN
      W1(I)=2.*FC
      ELSE
      W1(I)=SIN(2.*PI*FC*FLOAT(I))/(PI*FLOAT(I))*G(I)
      ENDIF
      AT=AT+W1(I)
   1  CONTINUE
      DO 10 I=-MOR,MOR
      W1(I)=W1(I)/AT
   10 CONTINUE

      FC=1./FLOAT(LAS)
      CALL KISER(G,MOR)
      AT=0.
      DO 2 I=-MOR,MOR
      IF (I.EQ.0) THEN
      W2(I)=2.*FC
      ELSE
      W2(I)=SIN(2.*PI*FC*FLOAT(I))/(PI*FLOAT(I))*G(I)
      ENDIF
      AT=AT+W2(I)
   2  CONTINUE
      DO 20 I=-MOR,MOR
      W2(I)=W2(I)/AT
   20 CONTINUE


      ZLENGTH=(2.*MOR+1)*DY
      FE=2.*PI/DY
      DO J=1,MOR
      FF=J*2.*PI/ZLENGTH/4.
      ZHF=0.
      DO I=-MOR,MOR
      ZHF=ZHF+(W2(I)-W1(I))*COS(2.*PI*I*FF/FE)
C     PRINT *,I*DY,W2(I)-W1(I)
      ENDDO
      PRINT *,FF/2./PI,ZHF
      ENDDO

      RETURN
      END



