c ############################
c Programme integration
c Integration d'un signal y(x)
c ############################

       implicit none
       integer i,j,ii,n,nmax,col1,col2,ncom,opt
       parameter (nmax=100000)
       double precision X(nmax),TAB(20),Y(nmax),DY(nmax),dsplines,deriv      
       character*140 c_ncom,c_col1,c_col2,c_fichin,com,c_opt

               
       call getarg(1,c_fichin)  ! Entrer le nom fichier
       
       if (Len_trim(c_fichin).eq.0) then
       write(6,*) ' <<	 Derivation >>'
       write(6,*) ' Usage : deriva [file][ncom][col1][col2][opt]'
       write(6,*) ' [file]....fichier de donnees'
       write(6,*) ' [ncom]....nombre de lignes de commentaires'
       write(6,*) ' [col1]....colonne des dates'
       write(6,*) ' [col2]....colonne des valeurs'
       write(6,*) ' [opt].....0 : derivation sur trois points'
       write(6,*) '      .....1 : derivation par splines cubiques'
       stop	
       endif

       call getarg(2 ,c_ncom ) 
       call getarg(3 ,c_col1 )  
       call getarg(4 ,c_col2 )  
       call getarg(5 ,c_opt )  
   
       read(c_ncom, *) ncom
       read(c_col1, *) col1
       read(c_col2, *) col2
       read(c_opt, *) opt

       open(10,file=c_fichin,access='sequential',status='old')

       do i=1,ncom
       read(10,1020) com
       enddo
1020   format(a140)
 
       ii = 0
       do  i=1,nmax
     	 read(10,*,end=20 ) (tab(j),j=1,max(col1,col2))    
         ii = ii + 1
         x(ii) = tab(col1)
         y(ii) = tab(col2)	 
       enddo	 

20     continue     
      
       close(10)
       n= ii
      

       if (opt.eq.0) then 
     	 	write(6,1030) X(1),(Y(2) - Y(1)) / ( X(2) - X(1) )     
      		do i = 2, n-1     	
      		write(6,1030) X(i),(Y(i+1) - Y(i-1)) / (X(i+1) - X(i-1))       
      		enddo       
      		write(6,1030) X(n),(Y(n) - Y(n-1)) / (X(n) - X(n-1))               
       endif      

       if (opt.eq.1) then
      		do i = 1, n 
       		deriv= dsplines(X,Y,n,X(i))
		write(6,1030) X(i),deriv    
    		enddo
       endif
      
      
1030   format(f18.6,f18.6)     
       end

c #######################################################################
c Programme dsplines:calcul de la derivee premiere d'une serie par 
c fonctions splines cubiques completes 
c #######################################################################
	   FUNCTION dsplines(x,y,n,xs)

c entree :   x        : tableau des dates
c            xs       : date a interpolee
c            y        : tableau des valeurs a interpoler par la fonction splines
c            n        : dimensions des tableaux x et y
c            ns       : dimensions des tableaux xs et ys
c	     alpha    : valeur de la derivee du signal en t0
c            beta     : valeur de la derivee du signal en tn 
c sortie     dsplines :  valeurs interpolees a la date dates xs
c #######################################################################
	   implicit none
	   integer i,n,j,n1
	   parameter (n1=100000)		
	
	   double precision x(n),y(n)
	   double precision S(n1),T(n1),d(n1),a(n1),b(n1),c(n1)
	   double precision M(n1),h(n1)
	   double precision alpha,beta,xs,dsplines
	

c La determination des coefficients des polynomes qui constituent la fonction
c splines consiste a resoudre un systeme tridiagonal par la methode du double
c balayage de Cholesky

	   do i=1,n-1
	   h(i) = x(i+1)-x(i)
	   enddo

	   alpha =(y(2)-y(1))/h(1)
	   beta  =(y(n)-y(n-1))/h(n-1)
	
	   b(1) = h(1)/3.d0
	   c(1) = h(1)/6.d0
	   d(1)=(y(2)-y(1))/h(1)-alpha
	
	   S(1) =  d(1)/b(1)
	   T(1) = -c(1)/b(1)

	   do i=2,n-1	
	
	   a(i) = h(i-1)/6.d0
	   b(i) = (h(i-1)+h(i))/3.d0
	   c(i) = h(i)/6.d0
	   d(i) = (y(i+1)-y(i))/h(i) - (y(i)-y(i-1))/h(i-1)
	
	   S(i) = (d(i)-a(i)*S(i-1)) / (b(i)+a(i)*T(i-1))
	   T(i) =-c(i)/(b(i)+a(i)*T(i-1))	
	   enddo
	
	   a(n) = h(n-1)/6.d0
	   b(n) = h(n-1)/3.d0	
	   d(n) = beta - (y(n)-y(n-1))/h(n-1)
	
c calcul des coeeficients M(i)
		
	   M(n)=(d(n)-a(n)*s(n-1))/(a(n)*T(n-1)+b(n))

	   do i=n-1,1,-1
	   M(i)=S(i)+T(i)*M(i+1)	 
	   enddo

	   if  (xs.lt.x(1)) then
	   j = 1
		  dsplines = 1.d0 / h(j)
     & 		*( -M(j)  /2.d0 * (x(j+1)-xs)**2  
     &  	   +M(j+1)/2.d0 * (x(j  )-xs)**2 
     &  	  -( y(j  ) -M(j)   *h(j)**2 /6.d0 )  
     &   	  +  y(j+1) -M(j+1) *h(j)**2 /6.d0 )
 
	   endif
	
	   if  ( (xs.ge.x(1)).and.(xs.lt.x(n-1)) ) then
		do j=1,n-2
		if  ( (xs.ge.x(j)).and.(xs.lt.x(j+1)) ) then

		  dsplines = 1.d0 / h(j)
     & 		*( -M(j)  /2.d0 * (x(j+1)-xs)**2  
     &  	   +M(j+1)/2.d0 * (x(j  )-xs)**2 
     &  	  -( y(j  ) -M(j)   *h(j)**2 /6.d0 )  
     &   	  +  y(j+1) -M(j+1) *h(j)**2 /6.d0 )
          
		endif
		enddo
	   endif
	
	   if (xs.ge.x(n-1)) then
        j = n-1 
	 	  dsplines = 1.d0 / h(j)
     & *( -M(j)  /2.d0 * (x(j+1)-xs)**2  
     &    +M(j+1)/2.d0 * (x(j  )-xs)**2 
     &    -( y(j  ) -M(j)   *h(j)**2 /6.d0 )  
     &    +  y(j+1)-M(j+1)  *h(j)**2 /6.d0 )
          
	   endif

	
	   return

	   end    


