      program vondrak

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Ce programme est une refonte de "vorace"
c avec entre autre un appel a la subroutine choix_eps 
c qui permet de ne plus trop se poser de questions 
c quant au choix de epsilon (voir publis Vondrack)
c et le rajout des possibilites d'un filtre passe-haut
c passe-bas , et passe-bande .............
c j'ai aussi effectue un decoupage qui permet de gagner de la memoire
c et de reutiliser les differentes parties du programme vorace
c en vue de faire le programme convolv....
c
c               Sylvain Loyer (1995)
c Derniere version : Mars 2006 - Christian Bizouard.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c en cas de probleme de taille de fichier modifier ndimension
c dans toutes les subroutines ...
      implicit none
      integer ndimension,nbcmax
      parameter(ndimension=100000,nbcmax=20) 
      double precision x(ndimension),y(ndimension),
     &  	       xp(ndimension),yl(ndimension),
     &                 ylisse(20,ndimension)
      double precision a1(nbcmax), P0, tr             
      integer n, itp, ichoix
      integer i,j,it1, nkom, cold, col, ifiltre, itrend, ncol
      character*80 fich, c_fich , c_nkom, c_cold, c_col, c_tr, c_P0, 
     &             c_ifiltre, c_trend, komment


c ######################       
c Arguments du programme       
c ######################
       
      	call getarg(1,c_fich     )  ! Entrer le nom fichier
        
	if (Len_trim(c_fich).eq.0) then
	write(6,*) '<<            Vondrak filter                >>'
	write(6,*) 'Usage : vondrak'
        write(6,*) '[file][ncom][cold][col][TR][P0][filtre][trend]'
        write(6,*) ' [file] ...... nom du fichier                  '
        write(6,*) ' [ncom] ...... nombre de lignes de commentaire '
        write(6,*) ' [cold] ...... colonne des  dates              '
        write(6,*) ' [col] ....... colonne des valeurs             '
        write(6,*) ' [TR]  ....... Coefficient de transfer en %    '       
        write(6,*) ' [P0]  ....... Periode pour TR % du signal     '       
        write(6,*) ' [filtre] .... garde hautes periodes T>Po -> 1 '
	write(6,*) '          .... garde basses periodes T<Po -> 2 '         
	write(6,*) '          .... filtre passe bande    T=Po -> 3 '
	write(6,*) '          .... excluer la bande      T=Po -> 4 ' 
	write(6,*) ' [trend]  .... pas de retrait             -> 0 ' 
	write(6,*) '          .... retrait parabole           -> 1 ' 
	write(6,*) '          .... retrait valeur moyenne     -> 2 ' 
        write(6,*) ' Sortie ecran - rediriger vers fichier par >   '   
	stop     
	endif
	
	call getarg(2 ,c_nkom    )  ! Nombre de lignes de commentaires
        call getarg(3 ,c_cold    )  ! Numero de la colonne des dates
        call getarg(4 ,c_col     )  ! Numero de la colonne des valeurs
        call getarg(5 ,c_tr     )   ! coefficient de transfer
        call getarg(6 ,c_P0      )  ! Periode de coupure a mi-hauteur
        call getarg(7 ,c_ifiltre )  ! Nature du filtrage 
        call getarg(8 ,c_trend )    ! retrait d'une tendance 
   
	fich=c_fich
	read(c_nkom,*)      nkom
	read(c_cold,*)      cold
	read(c_col, *)      col
	read(c_tr,  *)      tr
	read(c_P0,  *)      P0
	read(c_ifiltre,*)   ifiltre     
	read(c_trend,*)     itrend     
	
c ------------------
c lecture du fichier 
c ------------------
      open(10,file=fich,status='old')
      do i=1,nkom
        read(10,'(a80)') komment
      enddo

      ncol = max(cold,col)
	      
      do i = 1, ndimension
      read(10,*,end=90) (a1(j),j=1,ncol)
      x(i)=a1(cold)   
      y(i)=a1(col)   
      enddo
 90   n=i-1
 
      if(n.lt.4.or.n.gt.ndimension) stop 'stop convvor2 : n trop grand'
      itp=1
 
      itp=0
      
      do i=1,n
      xp(i)=1.e0
      enddo
      
      close(10)

c --------
c vorace 
c --------   
      call vorac(x,y,xp,n,tr,P0,ifiltre,itrend,yl) 

c ------------------
c ecriture resultats 
c -------------------
40    continue 
 
c      open(11,file='vondrak.dat')
      do i = 1,n 
      write(6,1020) x(i),yl(i)       
1020  format(f15.4,e18.8)
      enddo
c      close(11)

      end 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vorac(x,y,xp,n,tr,P0,ifiltre,itrend,yl)  
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c             (mange le bruit)
c
c             lissage de vondrak ameliore. le degre de lissage ne depend
c             ni de la ponderation,ni des intervalles. il est uniquement
c             fixe par epsilon. toutefois, si la mesure du temps change,
c             il faut se rappeler que epsilon a la dimension t**-6.
c
c             argument x(i), donnee y(i), poids p(i)
c             coefficient de lissage eps
c             resultat du lissage : yl(i)
c             nombre de donnees : n
c             perequation des poids sur : 2*kmoy + 1 valeurs
c
c Entrees :
c  x(i) .... tableau des dates des valeurs brutes 
c  y(i) .... liste des valeurs brutes de la serie a lisser
c  xp(i) ... poids des valeurs brutes individuelles
c  n ....... nombre de valeurs de la serie a lisser (2*kmoy+1 < n < 110000)
c  tr ...... coefficient de transfer en %
c  P0 ...... periode de demi-coupure
c  ifiltre . nature du filtrage
c
c Sorties :
c  yl ...... valeurs lissees pour les dates x(i)
c
c ==========================================================================
      implicit none
      integer ndimension
      parameter(ndimension=100000)
      double precision x(ndimension),y(ndimension),yl(ndimension)
      double precision xp(ndimension)
      double precision ypar(ndimension),yl2(ndimension)
      double precision xd,xparam(3),vmoy,val1,val2
      double precision ecmoy,xret,sr,tr,P0,eps1,eps2
      integer it 
      integer kmoy,n,itrend,ifiltre,i

c fonctions
      double precision xecmoy

c***** quelques stats
      sr=x(n)-x(1) 
      ecmoy=xecmoy(x,n)    
      call decal(x,n,xd)

c***** preparation *****
      call choix_eps(sr, ifiltre, tr, P0, eps1, eps2)
      eps1 = eps1/dble(n-3)
      eps2 = eps2/dble(n-3)
	
      it = 1  	
      if (it.eq.1)  then 
        kmoy=10
      call pereq(kmoy,x,xp,n) ! perequation des poids
      endif 
    
      call ajustpar(x,y,n,ypar,xparam) ! ajustement de la parabole
       
c***** resolution du systeme  *****
      call etasolv(x,y,xp,yl,eps1,n,sr) 
 
      if (ifiltre.ge.3) then
         call etasolv(x,y,xp,yl2,eps2,n,sr)
         do i=1,n
           yl(i)=yl(i) - yl2(i)
           if (ifiltre.eq.4) yl(i)=y(i)-yl(i)
         enddo 
      endif 

c**** calcul de la valeur moyenne ******
        call xvalmoy(y,n,vmoy)
    	
        if (itrend.eq.0) then 
	val1 = 0.d0   
        val2 = 0.d0 
	endif
	
	if (itrend.eq.1) then   ! retrait de la parabole
	val1 = 1.d0  
        val2 = 0.d0  
	endif
	
	if (itrend.eq.2) then  ! retrait moyenne
	val1 = 0.d0  
        val2 = 1.d0  
	endif
      
        do i=1,n 
     
        x(i)=x(i)+xd
        if (ifiltre.eq.1) then   
        yl(i) = yl(i) - (ypar(i) * val1 + val2 * vmoy)
        endif
	
        if (ifiltre.eq.2) yl(i)=y(i)-yl(i)  ! selectionne les hautes frequences
      enddo 
   
      return 
      end 


ccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine xvalmoy(y, n, vmoy)
ccccccccccccccccccccccccccccccccccccccccccccccccc
c calcul valeurs moyennes sur y et y+ypar 
c et affiche valeur de la parabole 
c (fonctionne avec ajustpar)
cccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none 
      integer ndimension
      parameter(ndimension=100000)
      double precision y(ndimension)
      double precision ypar(ndimension)
      double precision valmoy,xdim,xparam(3)
      double precision vmoy
      integer n,i,k
      
      xdim=dble(n)
      
      valmoy = 0.d0
      do i=1,n
        valmoy = y(i)+valmoy
      enddo
       
      vmoy = valmoy/xdim

      end 
 
ccccccccccccccccccccccccccccccccccccccccccccc
      subroutine ajustpar(x,y,n,ypar,xparam)
c Ajustement initial d'une parabole sur y(x) 
c xparam contient les parametres ajustes
c et ypar les valeurs tabulees de la parabole 
ccccccccccccccccccccccccccccccccccccccccccccc 
      implicit none 
      integer ndimension
      parameter(ndimension=100000)
      double precision x(ndimension)
      double precision y(ndimension)
      double precision ypar(ndimension)
      double precision gx(ndimension)
      double precision gk(ndimension)
      double precision gpx2,gpx3,gpk1,gpx4,gpk2
      double precision h,q,v,xparam(3),xdim
      integer n,i 
      xdim=float(n)
       do 7 i=1,4
       gx(i)=0.e0
7      gk(i)=0.e0
       do 8 i=1,n
       gx(1)=gx(1)+x(i)
       gx(2)=gx(2)+x(i)*x(i)
       gx(3)=gx(3)+x(i)**3
       gx(4)=gx(4)+x(i)**4
       gk(3)=gk(3)+y(i)
       gk(1)=gk(1)+y(i)*x(i)
8      gk(2)=gk(2)+y(i)*x(i)*x(i)
       gpx2=gx(2)-gx(1)*gx(1)/xdim
       gpx3=gx(3)-gx(1)*gx(2)/xdim
       gpk1=gk(1)-gk(3)*gx(1)/xdim
       gpx4=gx(4)-gx(2)*gx(2)/xdim
       gpk2=gk(2)-gk(3)*gx(2)/xdim
       h=(gpk2-gpk1*gpx3/gpx2)/(gpx4-gpx3*gpx3/gpx2)
       q=(gpk1-h*gpx3)/gpx2
       v=(gk(3)-gx(2)*h-gx(1)*q)/xdim
       do i=1,n
       ypar(i)=v+q*x(i)+h*x(i)*x(i)
       enddo
       xparam(1)=h
       xparam(2)=q
       xparam(3)=v
       return 
       end 

cccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine pereq(kmoy,x,xp,n)
c fait la perequation des poids (voir Guinot )
c et ecrase xp avec les nouveaux poids
cccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none 
      integer ndimension 
      parameter(ndimension=100000)
      double precision x(ndimension),xp(ndimension)
      double precision pm,sr,xpu(ndimension)
      integer n,kmoy,kmoi
      integer kp1,nmk,i,ir,is,iq

      kmoi=(kmoy*2)+1
c       if(kmoi.ge.n) write(6,604) kmoy,n
604    format(/10x,' kmoy=',i3,'  n= ',i5)
       if(kmoi.ge.n) go to 1000

       kp1=kmoy+1
       nmk=n-kmoy
       do 5 i=1,n
        pm=0.e0
        if(i.lt.kp1) go to 4
        if(i.gt.nmk) go to 6
        ir=i-kmoy
        go to 10
4       ir=1
        go to 10
6       ir=n-2*kmoy
10      is=ir+2*kmoy
        do 11 iq=ir,is
11       pm=pm+xp(iq)/(2.*kmoy+1.)
c        if (pm.eq.0.)write(6,605) kmoy,iq,x(iq),x(1)
5      xpu(i)=xp(i)/pm
       do i=1,n 
        xp(i)=xpu(i)
       enddo 

       return 
1000   write(6,*) " Echec de Vorace dans sub. pereq..."
       stop  
605    format(//'  La moyenne des poids de',i5,' valeurs autour de la',
     . i5,'eme,date',f6.2,' est nulle - date origine',f8.2/)
      end 

cccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine decal(x,n,xd)
c
c  decalage des dates 
c et tests sur le tableau x
cccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none 
      integer ndimension 
      parameter(ndimension=100000)
      double precision x(ndimension)
      double precision xd,dxi
      integer n,i
      xd=x(1)
      x(1)=0.e0
      do i=2,n 
        x(i)=x(i)-xd
        dxi=x(i)-x(i-1)
        if (dxi.le.0.D0) then 
          write(6,900) i,x(i),xd
          stop
        endif 
      enddo 
      return
900   format(' Dates non croissantes',i8,f8.2,
     . ' date origine:',f15.5/)
      end 


cccccccccccccccccccccccccccccccccccccccccccccc
      double precision function xecmoy(x,n) 
c calcul l'ecart moyen entre les points 
c successifs de x(1:n) 
ccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none                         
      integer ndimension
      parameter(ndimension=100000)
      double precision x(ndimension),xdim
      integer i,n 
      xdim=float(n-1)
      xecmoy=0.e0
      do i=2,n 
       xecmoy=xecmoy+x(i)-x(i-1)
      enddo 
      xecmoy=xecmoy/xdim
      return
      end 

c ====================================================
      subroutine choix_eps(xr, ifiltre, tr, P0, eps1, eps2)
c
c permet de choisir le bon epsilon pour appeler etasolv 
c voir  publis Vondrack
c 
c entrees : xr .......... duree totale    
c           ifiltre...... nature du filtrage
c           tr .......... coefficient de transfer en %
c           P0 .......... periode de demi-coupure
c
c sorties : eps1......... parametre de filtrage
c           eps2......... parametre de filtrage
c
c ================================================================
      implicit none  
      double precision eps1,eps2,xr
      integer ifiltre
      double precision xdim, P0, pi, Tmax, Tmin, tr
   
      pi = 4.d0 * datan(1.d0)
      	           
      Tmax = tr/100.d0
      Tmin = 1.d0 - Tmax 
     
      eps1 = (2.d0*pi/P0)**6 * Tmax / (1.d0 - Tmax) 
  	  
      if (ifiltre.le.2) then
      if (P0/xr.ge.1.d0) 
     &write(6,'(" Periode de coupure > duree du signal  0.1 !")')
      eps2 = 0.d0   
      endif
   
      if (ifiltre.ge.3) then
      eps2 = (2.d0*pi/P0)**6 * Tmin / (1.d0 - Tmin) 
      endif	
     
c      write(6,'("    eps1=",d10.2,"   eps2=",d10.2)') eps1,eps2   

      return
      end

c ===========================================
       subroutine etasolv(x,y,pu,w,eps,n,sr)
c       
c      etablissement et resolution du systeme
c ===========================================
       implicit none
       integer ndimension 
       parameter (ndimension=100000)
       integer n,i,j,ndep,nlim,l,ls
       integer i1,j1,i2,j2,il1,il,jl,in
       double precision eps,coef
       double precision sr,delx,aux
       double precision x(ndimension),y(ndimension)        
       double precision pu(ndimension),
     $ aa(7,ndimension),w(ndimension)
      double precision xpp(4,4)
      do j=1,3
        do i=1,4
         xpp(i,j)=0.
        enddo
      enddo 
      nlim=n-2
      ls=n-3
      do i=1,n
        if (i.gt.ls) then
         do l=1,4
           xpp(l,4)=0. 
         enddo  ! l
        else
         delx=x(i+2)-x(i+1)
         aux=(6.e0*sqrt(delx))/sqrt(sr)
         xpp(1,4)=aux/((x(i)-x(i+1))*(x(i)-x(i+2))*(x(i)-x(i+3)))
         xpp(2,4)=aux/((x(i+1)-x(i))*(x(i+1)-x(i+2))*(x(i+1)-x(i+3)))
         xpp(3,4)=aux/((x(i+2)-x(i))*(x(i+2)-x(i+1))*(x(i+2)-x(i+3)))
         xpp(4,4)=aux/((x(i+3)-x(i))*(x(i+3)-x(i+1))*(x(i+3)-x(i+2)))
        endif
        aa(1,i)=xpp(1,1)*xpp(4,1)
        aa(2,i)=xpp(1,2)*xpp(3,2)+xpp(2,1)*xpp(4,1)
        aa(3,i)=xpp(1,3)*xpp(2,3)+xpp(2,2)*xpp(3,2)+xpp(3,1)*xpp(4,1)
        aa(4,i)=eps*pu(i)  +xpp(1,4)*xpp(1,4)+xpp(2,3)*xpp(2,3)
     .  +xpp(3,2)*xpp(3,2)+xpp(4,1)*xpp(4,1)
        aa(5,i)=xpp(1,4)*xpp(2,4)+xpp(2,3)*xpp(3,3)+xpp(3,2)*xpp(4,2)
        aa(6,i)=xpp(1,4)*xpp(3,4)+xpp(2,3)*xpp(4,3)
        aa(7,i)=xpp(1,4)*xpp(4,4)
        do j=1,3
          do l=1,4 
            xpp(l,j)=xpp(l,j+1) 
          enddo 
        enddo 
        w(i)=eps*pu(i)*y(i)
      enddo   
      ndep=0
  707 i1=4
      if(ndep) 702,701,702
  701 j1=1
      ndep=1
      go to 703
  702 j1=j1+1
      if(j1-nlim) 703,708,703
  703 continue
      i2=3
      j2=j1+1
  704 coef=aa(i2,j2)/aa(i1,j1)
      ls=i2+3
      il1=3
      do il=i2,ls
        il1=il1+1
        aa(il,j2)=aa(il,j2)-aa(il1,j1)*coef 
      enddo 
      w(j2)=w(j2)-w(j1)*coef
      i2=i2-1
      j2=j2+1
      if(i2) 704,706,704
  706 continue
      go to 707

  708 continue
      i1=4
      j1=n-2
      i2=3
      j2=n-1
  709 ls=i2+2
      coef=aa(i2,j2)/aa(i1,j1)
      il1=3
      do 710 il=i2,ls
      il1=il1+1
  710 aa(il,j2)=aa(il,j2)-aa(il1,j1)*coef
      w(j2)=w(j2)-w(j1)*coef
      i2=i2-1
      j2=j2+1
      if(i2-1) 709,711,709
711   continue
      j1=n-1
      coef=aa(3,n)/aa(4,j1)
      aa(4,n)=aa(4,n)-aa(5,j1)*coef
      w(n)=w(n)-w(j1)*coef
      w(n)=w(n)/aa(4,n)
      w(j1)=(w(j1)-aa(5,j1)*w(n))/aa(4,j1) 
      j2=n-2
      w(j2)=(w(j2)-aa(5,j2)*w(j1)-aa(6,j2)*w(n))/aa(4,j2)
      jl=n-3
      do j=1,jl
      in=j2-j
      w(in)=(w(in)-aa(5,in)*w(in+1)-aa(6,in)*w(in+2)-aa(7,in)*
     .w(in+3))/aa(4,in)
      enddo 
      end
c ===================================================
	subroutine maxmin(x,y,imax,xmin,xmax,ymin,ymax)
c calcul du minimum et maximun de x et y
c ===================================================
                  implicit none
                  integer ndimension 
                  parameter(ndimension=100000) 
                  double precision xmin,xmax,ymin,ymax
                  integer i,imax
 	double precision x(ndimension),y(ndimension)
	xmin=x(1)
	xmax=x(1)
	ymin=y(1)
	ymax=y(1)
	do 16 i=2,imax
	if(x(i).lt.xmin) xmin=x(i)
	if(x(i).gt.xmax) xmax=x(i)
	if(y(i).lt.ymin) ymin=y(i)
	if(y(i).gt.ymax) ymax=y(i)
16	continue
	 write(6,30) xmin,xmax,ymin,ymax
30	 format(/1x,'xmin=',f18.6/1x,'xmax=',f18.6
     .   /1x,'ymin=',f18.6/1x,'ymax=',f18.6)
	return
	end 
