      program vonda

c---------------------------------------------------------------------------
c Filtrage de Vondrak de plusieurs tableaux de valeurs ayant les memes dates
c Derniere version : Novembre 2008 - Christian Bizouard.
c --------------------------------------------------------------------------

      implicit none
      integer ndim,nbcmax
      parameter(ndim=100000,nbcmax=20) 
      double precision X(ndim) ,Y(ndim) ,Y1(ndim,nbcmax),
     &  	       XP(ndim),YL(ndim),YL1(ndim,nbcmax)
 
      double precision a1(nbcmax), p1, p2, tr1,tr2             
      integer n, itp
      integer i,j,k,it1, nkom, cold, col, coli, colf, ifiltre, itrend, ncol
      character*80 fichin, c_fich , c_fichin,c_nkom, 
     &                   c_cold, c_coli, c_colf, c_tr1, c_p1, c_tr2, c_p2,
     &                   c_ifiltre, c_trend, komment


c ######################       
c Arguments du programme       
c ######################
       
      	call getarg(1,c_fichin     )  ! Entrer le nom fichier
        
	if (Len_trim(c_fichin).eq.0) then
	write(6,*) '<<            Filtre de Vondrak multicolonne    >>'
	write(6,*) ' Usage : vonda [filein][ncom][cold][coli][colf][TR1][P1][TR2][P2][filtre][trend]'
        write(6,*) ' [filein] ..... nom du fichier a filtrer'
        write(6,*) ' [ncom] ....... nombre de lignes de commentaire '
        write(6,*) ' [cold] ....... colonne des  dates              '
        write(6,*) ' [coli] ....... 1ere colonne des valeurs        '
        write(6,*) ' [colf] ....... derniere colonne des valeurs    '
        write(6,*) ' [TR1]  ....... Coefficient de transfer en %         (0 par defaut)'	 
        write(6,*) ' [P1]  ........ Periode pour TR1 % du signal         (0 par defaut)'	 
        write(6,*) ' [TR2]  ....... Coefficient de transfer en %         (0 par defaut)'       
        write(6,*) ' [P2]  ........ Periode pour TR2 % du signal P2 > P1 (0 par defaut)'       
        write(6,*) ' [filtre] ..... selectionne periodes [P1,infini]        -> 1'
	write(6,*) '          ..... selectionne periodes [0,P2]             -> 2 '         
	write(6,*) '          ..... selectionne periodes [P1,P2]            -> 3 '
	write(6,*) '          ..... selectionne periodes [0,P1]U[P2,infini] -> 4 ' 
	write(6,*) ' [trend]  ..... pas de retrait de tendance -> 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_coli    )  ! Numero de la 1ere colonne des valeurs
        call getarg(5 ,c_colf    )  ! Numero de la derniere colonne des valeurs
        call getarg(6 ,c_tr1     )  ! coefficient de transfer
        call getarg(7 ,c_p1      )  ! Periode de coupure a mi-hauteur
        call getarg(8 ,c_tr2     )  ! coefficient de transfer
        call getarg(9,c_p2      )  ! Periode de coupure a mi-hauteur
        call getarg(10,c_ifiltre )  ! Nature du filtrage 
        call getarg(11,c_trend   )  ! retrait d'une tendance 
   
	fichin  = c_fichin
	
	read(c_nkom, *) nkom
	read(c_cold, *) cold
	read(c_coli, *) coli
	read(c_colf, *) colf
	read(c_tr1,  *) tr1
	read(c_p1,   *) p1
	read(c_tr2,  *) tr2
	read(c_p2,   *) p2
	read(c_ifiltre,*) ifiltre     
	read(c_trend,*)   itrend     
	
c ------------------
c lecture du fichier 
c ------------------
       open(10,file=fichin,status='old')
       do i=1,nkom
        read(10,'(a80)') komment
       enddo

       ncol = max(cold,colf)
       	      
       do i = 1, ndim
       read(10,*,end=90) (A1(j),j=1,ncol)
       X(i) = A1(cold)
            
       	  do j = coli, colf        
          Y1(i,j) = A1(j)   
          enddo
       enddo
       
 90    n=i-1
 
       if(n.lt.4)    stop 'stop vondrakmulti : n < 4 '
       if(n.gt.ndim) stop 'stop vondrakmulti : n trop grand'

       itp=1
 
       itp=0
      
       do i=1,n
       XP(i)=1.d0
       enddo
      
       close(10)


c------------------------
c boucle sur les colonnes
c -----------------------
	do k = coli,colf
	
		do i=1,n
		Y(i) = Y1(i,k) 
		enddo	   
c --------
c vorace 
c --------   
     		call vorac(X,Y,XP,n,tr1,p1,tr2,p2,ifiltre,itrend,YL) 
                 
         	do i=1,n
		YL1(i,k) = YL(i) 
		enddo
	 
	 enddo
         
c ------------------
c ecriture resultats 
c -------------------
   		 
     		  do i = 1,n 
     		  write(6,1020) X(i),(YL1(i,k), k=coli,colf)       
1020  		  format(f15.4,20e18.8)
      		  enddo     	
       
      end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vorac(X,Y,XP,n,tr1,p1,tr2,p2,ifiltre,itrend,YL)  
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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... .... tableau des dates des valeurs brutes 
c  Y ....... liste des valeurs brutes de la serie a lisser
c  XP ...... poids des valeurs brutes individuelles
c  n ....... nombre de valeurs de la serie a lisser (2*kmoy+1 < n < 110000)
c  tr1 ..... coefficient de transfer en % pur la periode p1
c  p1 ...... periode coupure a TR1 %
c  tr2 ..... coefficient de transfer en % pour la periode p2
c  p2 ...... periode coupure a TR2 %
c  ifiltre . nature du filtrage
c
c Sorties :
c  YL ...... valeurs filtrees / lissees pour les dates x(i)
c
c ==========================================================================
       implicit none
       integer ndim
       parameter(ndim=100000)
       double precision X(ndim),Y(ndim),YL(ndim),YL1(ndim),YL2(ndim)
       double precision xp(ndim),YPAR(ndim),XPARAM(3)
       double precision xd,vmoy,val1,val2
       double precision ecmoy,xret,tr1,p1,tr2,tr22,p2,eps1,eps2,pi
       integer it 
       integer kmoy,n,itrend,ifiltre,i
	       
       pi = 4.d0 * datan(1.d0)

c decalage des dates pour qu'elles commencent a zero avant filtrage
       call decal(X,n,xd)

c perequation des poids 
	it = 1  	
       	if (it.eq.1)  then 
        kmoy=10
       	call pereq(kmoy,X,XP,n)  
       	endif 
	
        ! Selectionne periodes dans l'intervalle [P1,infini]  avec coefficient de transfert TR1
        if (ifiltre.eq.1) then                    	     
        	eps1 = (2.d0*pi/p1)**6 * tr1 / (100.d0 - tr1)                      
       		eps1 = eps1/dble(n-3)    	   	
       
c lissage / selection des periodes > P1 
       		call etasolv(X,Y,XP,YL1,eps1,n)  ! filtrage laissant passer periodes > p1 avec un taux > tr1 % --> YL              	    	
c ajustement de la parabole    
                call ajustpar(X,Y,n,YPAR,XPARAM) 

c calcul de la valeur moyenne 
        	call valmoy(Y,n,vmoy)
    	
	        if (itrend.eq.0) then  ! pas de retrait 
		val1 = 0.d0   
	        val2 = 0.d0 
		endif
	
		if (itrend.eq.1) then  ! retrait de parabole
		val1 = 1.d0  
	        val2 = 0.d0  
		endif
	
		if (itrend.eq.2) then  ! retrait de valeur moyenne 
		val1 = 0.d0  
	        val2 = vmoy   
		endif
                              
	        do i=1,n 
		X(i)  = X(i) + xd  ! on reforme le tableau initiales 
		YL(i) = YL1(i) - (YPAR(i) * val1 + val2)  
	        enddo
	endif

	! Selectionne periodes  dans l'intervalle [0,P2] avec coefficient de transfert TR2 pour P2
        if (ifiltre.eq.2) then 
       		tr22 = 100.d0 - tr2
       		eps2 = (2.d0*pi/p2)**6 * tr22 / (100.d0 - tr22)       
       		eps2 = eps2 / dble(n-3)
		
		call etasolv(X,Y,XP,YL2,eps2,n)  ! filtrage laissant passer periodes > p2 avec un taux > tr2 % --> YL2
       			do i=1,n   
       			YL(i) = Y(i) - YL2(i) 
       			enddo 
		do i=1,n 
		X(i)  = X(i) + xd  ! on reforme le tableau initiales 	
		enddo	
        endif
	
	! Selectionne periodes P  dans l'intervalle [P1,P2] P2 > P1 
	! avec coefficients de transfert TR1 pour P1 et TR2 pour P2 
        if ( (ifiltre.eq.3).or.(ifiltre.eq.4) ) then 
		! periodes superieures a > P1
        	eps1 = (2.d0*pi/p1)**6 * tr1 / (100.d0 - tr1)                      
       		eps1 = eps1/dble(n-3)    	      		        
                ! lissage / selection des periodes > P1 
       		call etasolv(X,Y,XP,YL1, eps1,n)  ! filtrage laissant passer periodes > p1 avec un taux > tr1 % --> YL       	          
		! filtrage des periodes < P2
		tr22 = 100.d0 - tr2 
                eps2 = (2.d0*pi/p2)**6 * (tr22) / (100.d0 - tr22)       
                eps2 = eps2/dble(n-3)
		call etasolv(X,YL1,XP,YL2,eps2,n)  ! filtrage laissant passer periodes > p2 avec un taux > tr2 % --> YL2
       		do i=1,n   
       		X(i)  = X(i) + xd  ! on reforme le tableau initiales 
		if (ifiltre.eq.3) YL(i) = YL1(i) - YL2(i) 
       		if (ifiltre.eq.4) YL(i) = Y(i) - ( YL1(i) - YL2(i)) 
		enddo 		
	endif
   
      return 
      end 


ccccccccccccccccccccccccccccccccccc
      subroutine valmoy(Y, n, vmoy)

c valeur moyenne du tableau Y
ccccccccccccccccccccccccccccccccccc
 
      implicit none 
      integer ndim
      parameter(ndim=100000)
      double precision Y(ndim)
      double precision vmoy
      integer n,i
      
      vmoy =0.d0
          
      do i=1,n
        vmoy = Y(i) + vmoy   
      enddo     
  
      vmoy= vmoy / dble(n) 
      end 
 
ccccccccccccccccccccccccccccccccccccccccccccc
      subroutine ajustpar(X,Y,n,YPAR,XPARAM)
c Ajustement initial d'une parabole sur y(x) 
c entrees 
c X ....... tableau des dates 
c Y ....... tableau des valeurs  
c n ....... dimension des tableaux X et Y
c
c sorties :
c XPARAM(3)... contient les parametres ajustes de la parabole : 
c              ypar = XPARAM(3) * x(i)^2 + XPARAM(2) * x(i) + XPARAM(1)  
c YPAR ....... valeurs tabulees de la parabole 

ccccccccccccccccccccccccccccccccccccccccccccc 
      implicit none 
      integer ndim
      parameter(ndim=100000)
      double precision x(ndim)
      double precision y(ndim)
      double precision ypar(ndim)
      double precision gx(ndim)
      double precision gk(ndim)
      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 ndim 
      parameter(ndim=100000)
      double precision X(ndim),XP(ndim)
      double precision pm,XPU(ndim)
      integer n,kmoy,kmoi
      integer kp1,nmk,i,ir,is,iq

      kmoi = (kmoy*2)+1

604    format(/10x,' kmoy=',i3,'  n= ',i5)
       if(kmoi.ge.n) go to 1000

       kp1 = kmoy + 1
       nmk = n - kmoy
       do i=1,n
          pm=0.d0
          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 iq = ir,is
          pm = pm + XP(iq) / (2.d0 * dble(kmoy) + 1.d0)
	  enddo

          XPU(i) = XP(i)/pm
       enddo
       
       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  date initiale mise a 0 et tests sur le tableau x
cccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none 
      integer ndim 
      parameter(ndim=100000)
      double precision X(ndim)
      double precision xd,dxi
      integer n,i
      xd   = X(1)
      X(1) = 0.d0
      
      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 

c ===========================================
       subroutine etasolv(X,Y,P,W,eps,n)
c       
c      etablissement et resolution du systeme
c ===========================================
       implicit none
       integer ndim 
       parameter (ndim=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(ndim),Y(ndim)        
       double precision P(ndim),AA(7,ndim),W(ndim)
       double precision XPP(4,4)
      
      sr = X(n) - X(1) 
       do j=1,3
        do i=1,4
         XPP(i,j) = 0.d0
        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.d0 
         enddo  ! l
        else
         delx=x(i+2)-x(i+1)
         aux=(6.d0*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*P(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 * P(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 ndim 
                  parameter(ndim=100000) 
                  double precision xmin,xmax,ymin,ymax
                  integer i,imax
 	double precision x(ndim),y(ndim)
	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 
