c----------------------------------------------------------------------------
c----------------------------------------------------------------------------
         SUBROUTINE DTREND(n,D,Y,P,a,sig_a,b,sig_b,cor_ab,sig0)
c
c objet : calcul du biais et de la pente d'une serie temporelle 
c         par regression lineaire et retrait de celle-ci 
c 
c entree    n     : dimension des tableaux D,X,P
c           D     : tableau des dates
c	    P     : tableau des poids affectant les valeurs
c
c entree/sortie   : Y signal / signal sans pente 
c
c sortie    a     : pente
c           sig_a : erreur sur la pente   
c           b     : biais
c           sig_b : erreur sur le biais
c           cor_ab: correlation entre a et b
c	    sig0  : ecart type apres retrait de la tendance lineaire 	
c 
c entree/sortie   : Y signal / signal sans pente 
c N.B. Cette procedure ne permet pas d'eliminer les valeurs anormales
c  
c Code par Ch. Bizouard en FORTRAN 90 - decembre 2002
c----------------------------------------------------------------------------
      
       IMPLICIT NONE
       INTEGER n,i
       DOUBLE PRECISION D(n),Y(n),P(n)    	
       DOUBLE PRECISION a,da,b,db, sig_a,sig_b,cor_ab
       DOUBLE PRECISION s_p,s_pd,s_py,s_pdy,s_pd2,det,sig0
  
       s_p   = 0.d0
       s_pd  = 0.d0 
       s_py  = 0.d0
       s_pdy = 0.d0
       s_pd2 = 0.d0
      
                 
       DO i=1,n
	      s_p  = s_p  + P(i)
	      s_pd = s_pd + P(i) * D(i)
	      s_py = s_py + P(i) * Y(i)
 	      s_pdy= s_pdy+ P(i) * D(i) * Y(i)
	      s_pd2= s_pd2+ P(i) * D(i)**2  
       ENDDO
      
       det = s_p * s_pd2 - s_pd**2
      
       a = ( s_p  * s_pdy - s_pd * s_py) / det  
       b = (-s_pd * s_pdy + s_pd2 * s_py)/ det
      
     
C Calcul des erreurs formelles et de la correlation     
      
       ! sig0 : ecart type par rapport a : Y = a D + b            
       sig0 = 0.d0
 
       DO i = 1,n
	      sig0 = sig0 + P(i) * (Y(i) - a * D(i) - b)**2       
              Y(i) = Y(i) - ( a * D(i) + b) 
       ENDDO
      
       sig0 = dsqrt(sig0/dble(n-1))  ! ecart type 
       
       sig_a = sig0 * dsqrt(dabs(s_p  /det))
       sig_b = sig0 * dsqrt(dabs(s_pd2/det)) 
      
       cor_ab= dabs(sig0**2 * s_pd/det)/(sig_a * sig_b)
      
       RETURN  
       END
      
c----------------------------------------------------------------------------
c----------------------------------------------------------------------------
	subroutine angle(x,y,phase)  
c computation of the angle with respect to Ox of the vector (x,y)  
c       input   :  x, y .... : coordinates of a vector
c	output  : phase .... : angle of the vector en degre / 0x  
c		                     
c----------------------------------------------------------------------------

	double precision x,y,phase,pi
  
	pi = 4.d0  *  atan(1.d0)  
                
	if (x.gt.0.d0 ) phase = datan(y/x)  
          
        if ( (x.eq.0.d0 ).and.(y.gt.0.d0) ) phase =  pi/2.d0    
        if ( (x.eq.0.d0 ).and.(y.lt.0.d0) ) phase = -pi/2.d0    
  
  	if (x.lt.0.d0 ) phase = datan(y/x) + pi  
    
c conversion in degres   
  
	phase  = mod(phase  * 180.d0 /pi,360.d0 )
	if  (phase.lt.0.d0 ) phase = phase + 360.d0 
	 
	return  
	end 

c----------------------------------------------------------------------------
c----------------------------------------------------------------------------
       FUNCTION TERVOD(D,Y,n,u)
            
c   INTERPOLATION PAR DES POLYNOMES DE LAGRANGE SUR 4 POINTS
c
c   D DATES DU TABLEAU A INTERPOLER
c   Y TABLEAU A INTERPOLER
c   n DIMENSION DE CES TABLEAUX
c   u DATE A LAQUELLE ON VEUT INTERPOLER

c ----------------------------------------------------------------------------
      IMPLICIT NONE	
      INTEGER m, n,i
      double precision D(n),Y(n),u, tervod,a1,a2,a3,a4

      m = n-1
      
      DO i=3,m
      IF (D(i).gt.u) GO TO 10    
      ENDDO
      
      IF (u.eq.D(n-1)) i = m
      IF (u.gt.D(n-1)) i = m      
              
 10   a1     = ( D(i-2) - D(i-1) ) * ( D(i-2) - D(i)   ) * ( D(i-2) - D(i+1) )
      a2     = ( D(i-1) - D(i-2) ) * ( D(i-1) - D(i)   ) * ( D(i-1) - D(i+1) )
      a3     = ( D(i)   - D(i-2) ) * ( D(i)   - D(i-1) ) * ( D(i)   - D(i+1) )
      a4     = ( D(i+1) - D(i-2) ) * ( D(i+1) - D(i-1) ) * ( D(i+1) - D(i)   )
 
      tervod = ( u - D(i-1) ) * ( u - D(i)   ) * ( u - D(i+1)) * Y(i-2) / a1
     .       + ( u - D(i-2) ) * ( u - D(i)   ) * ( u - D(i+1)) * Y(i-1) / a2
     .       + ( u - D(i-2) ) * ( u - D(i-1) ) * ( u - D(i+1)) * Y(i)   / a3    
     .       + ( u - D(i-2) ) * ( u - D(i-1) ) * ( u - D(i)  ) * Y(i+1) / a4
         
        
      RETURN
      END
      

c ----------------------------------------------------------------------------   
c ----------------------------------------------------------------------------   
 	 SUBROUTINE AJUS1D(T,X,SX,n,PER,nper,t0,npoly,XPARA,DXPARA,COR,XR)
c estimation de termes periodiques dans une serie temporelle a 1 dimension

c ENTREE :
c   T ........ dates de la serie
c   X ........ tableau des valeurs
c   SX ........tableau des incertitudes
c   n .........nombre de dates
c   t0 ....... date de reference 
c   PER ...... tableau des periodes
c   nper ..... nombre de periodes a estimer
c   t0 ....... date de reference
c   npoly .... degre du polynome a ajuster 

c SORTIE
c   XPARA .....tableau des parametres : les "npoly" parametres du polynome sont les "npoly" derniers ranges par degre corissant : biais - degre 1 -- 
c   DXPARA ....tableau des incertitudes sur les parametres
c   COR .......tableau des coefficients de correlation    
c   XR ........tableau des residus  
c ----------------------------------------------------------------------------   
       
	implicit none       
	integer n,i,j,nes,nper,npoly,k,l  
        DOUBLE PRECISION T(n),X(n),SX(n),PER(nper),t0	       
	DOUBLE PRECISION  A(nper),xpi,chi2
	DOUBLE PRECISION ZA(n,2*nper+npoly+1),ZP(n),PP(n)
	DOUBLE PRECISION XPARA(2*nper+npoly+1),DXPARA(2*nper+npoly+1),COR(2*nper+npoly+1,2*nper+npoly+1),XR(n) 
		                                             
	xpi=3.141592654       
      	do i = 1, n      
                                
c ---------------------       
c arguments periodiques    
c ---------------------       
		do k=1,nper 
		A(k) = 2.d0 * xpi / PER(k) * (T(i)-t0)  
		A(k)= dmod(A(k),2.D0 *xpi) 	
		enddo
		 
c --------------------------------        
c observations matrix zp  et poids       
c --------------------------------        		  
		  ZP(i) = X(i)             
		  PP(i) = 1.d0 /(SX(i)**2)     
	                
c -----------------------------------------------        
c observations relation matrix za :   
c          x = \sum [m_c cos a(j) + m_s sin a(j)]            
c -----------------------------------------------                  
     	          do j =1,nper       
	          ZA(i,2*j-1) =  dcos(A(j))  	   
	          ZA(i,2*j  ) =  dsin(A(j))  		                
		  enddo            

	          ! estimation d'un polynome de degre npoly		  
		  do j=1,npoly+1 
                  ZA(i,2*nper + j) = (T(i) -t0)**(j-1)  	  
		  enddo
	                 	       
        enddo              
	         
        l = 2*nper + npoly+1
	          
c########################        
c least square treatement        
c########################          
		          
	call mcarre(n,l,ZA,ZP, PP,XPARA,XR,DXPARA,COR,chi2,nes,0)        

	END        

c #############################################################################        
 	 SUBROUTINE AJUSPOISSON(T,X,SX,n,PER,nper,t0,npoly,XPARA,DXPARA,COR,XR)
c estimation de termes periodiques / de Poisson / polynomiaux dans une serie temporelle a 1 dimension

c ENTREE :
c   T ........ dates de la serie
c   X ........ tableau des valeurs
c   SX ........tableau des incertitudes
c   n .........nombre de dates
c   t0 ....... date de reference 
c   PER ...... tableau des periodes
c   nper ..... nombre de periodes a estimer
c   t0 ....... date de reference
c   npoly .... degre du polynome a ajuster 

c SORTIE
c   XPARA .....tableau des parametres
c   DXPARA ....tableau des incertitudes sur les parametres
c   COR .......tableau des coefficients de correlation    
c   XR ........tableau des residus  
c ##################################################################         
       
	implicit none       
	integer n,i,j,nes,nper,npoly,k,l  
        DOUBLE PRECISION T(n),X(n),SX(n),PER(nper)	       
	DOUBLE PRECISION  A(nper)
	DOUBLE PRECISION ZA(n,4*nper+npoly+1),ZP(n),PP(n)
	DOUBLE PRECISION XPARA(4*nper+npoly+1),DXPARA(4*nper+npoly+1),COR(4*nper+npoly+1,4*nper+npoly+1),XR(n) 
	DOUBLE PRECISION t0,dpi,chi2
	                                             
	dpi=8.d0*datan(1.d0)             
           
      	do i = 1, n      
                                
c ---------------------       
c arguments periodiques    
c ---------------------       
		do k=1,nper 
		A(k) = dpi / PER(k) * (T(i)-t0) 
		A(k)= dmod(A(k),dpi) 
		enddo
		 
c --------------------------------        
c observations matrix zp  et poids       
c --------------------------------        		  
		  ZP(i) = X(i)                 
		  PP(i) = 1.d0 /(SX(i)**2)     
	                
c -----------------------------------------------        
c observations relation matrix za :   
c          x = \sum [m_c cos a(j) + m_s sin a(j)]            
c -----------------------------------------------                  
     	          do j =1,nper       
	          ZA(i,4*j-3) =  dcos(A(j))  
		  ZA(i,4*j-2) =  (T(i) - t0) * dcos(A(j))  	   
	          ZA(i,4*j-1) =  dsin(A(j))
		  ZA(i,4*j  ) =  (T(i) - t0) * dsin(A(j)) 		                
	          enddo            

	          ! estimation d'un polynome de degre npoly		  
		  do j=1,npoly+1 
                  ZA(i,4*nper + j) = (T(i) -t0)**(j-1)  	  
	          enddo
	                 	       
       enddo              
	         
        l = 4*nper + npoly+1
	          
c########################        
c least square treatement        
c########################          
		          
	CALL MCARRE(n,l,ZA,ZP, PP,XPARA,XR,DXPARA,COR,chi2,nes,0)        

	END 

c ####################################################################        
 	 SUBROUTINE AJUSPOLY(T,X,SX,n,t0,npoly,XPARA,DXPARA,COR,XR)
c estimation d'un polynome de dimension npoly

c ENTREE :
c   T ........ dates de la serie
c   X ........ tableau des valeurs
c   SX ........tableau des incertitudes
c   n .........nombre de dates
c   t0 ....... date de reference 
c   npoly .....degre du polynonme a ajuster  

c SORTIE
c   XPARA .....tableau des parametres
c   DXPARA ....tableau des incertitudes sur les parametres
c   COR .......tableau des coefficients de correlation    
c   XR ........tableau des residus  
c ##################################################################         
       
	implicit none       
	integer n,i,j,nes,npoly,k,l
	
        DOUBLE PRECISION T(n),X(n),SX(n)	       
	DOUBLE PRECISION ZA(n,npoly+1),ZP(n),PP(n)
	DOUBLE PRECISION XPARA(npoly+1),DXPARA(npoly+1),COR(npoly+1,npoly+1),XR(n) 
	DOUBLE PRECISION t0,dpi,chi2
	                                             
	dpi=8.d0*datan(1.d0)               
      	do i = 1, n      
                                		 
c --------------------------------        
c observations matrix zp  et poids       
c --------------------------------        		  
		  ZP(i) = X(i)                 
		  PP(i) = 1.d0 /(SX(i)**2)     
	                		
c -------------------------------------------        
c observations relation matrix za :   
c          x = \sum [m_c cos a(j) + m_s sin a(j)]            
c -----------------------------------------------                  
     	          ZA(i,1) = 1.d0
		  do j =2,npoly+1     
                  ZA(i,j) = (T(i) -t0)**(j-1) 
	          enddo
	                 	       
       enddo       
         
	l = npoly+1    
        
       
c########################        
c least square treatement        
c########################          
		          
	CALL MCARRE(n,l,ZA,ZP, PP,XPARA,XR,DXPARA,COR,chi2,nes,0)        

	END                 
       
c ################################################################       
        SUBROUTINE MCARRE(n,l,ZA,ZP,PP,XS,XR,XDX,COR,chi2,NES,nopt)       
c ################################################################       
c input : n  : nombre d'observations
c         l  : nombre de parametre  
c        ZA  : matrice des relations d'observations
c        ZP  : observations
c        PP  : matrice des poids
c
c sortie XS  : parametres estimes  
c        XR  : residus  
c        XDX : erreurs sur les parametres estimes 
c        COR : matrice de correlation entre les parametres
c        chi2: estimation du chi2
c        Nes : nombre d'equation supprimees
c        nopt: supression d'equation --> 1 / pas de supression --> 0   
     
        implicit DOUBLE PRECISION (a-h,o-z)             
        integer n, l 
        DOUBLE PRECISION ZA(n,l),ZP(n),PP(n),XR(n),       
     $         XS(l),XDX(l),COR(l,l)   

        DOUBLE PRECISION ZA1(n,l),ZP1(n),ZS(n),ZB(l,l),ZC(l,l)       


        kt=0.d0        
        nes=0                                                
      
c #############################################        
c matrixs ZA1=ZA t(sqrt(ZP)) and ZP1=sqrt(ZP) L        
c #############################################        
        
 100    DO i=1,n        
        DO j=1,l        
        ZA1(I,J)=ZA(I,J)*DSQRT(PP(I))        
	ENDDO        
        ZP1(I)=ZP(I)*DSQRT(PP(I))        
	ENDDO        
       
        do i=1,l
        ZS(I)=0.D0  
        enddo
       
        DO 24 I=1,l        
        DO 25 J=1,l        
        ZB(I,J)=0.D0        
 25     ZC(I,J)=0.D0        
 24     ZC(I,I)=1.D0        
        
        DO I=1,l       
        XS(I)=0.D0        
        ENDDO        
        
c ############################        
c normal equations : matrix ZB        
c ############################        
        
        DO 20 I=1,L        
        
        DO 21 J=1,L        
        
        DO 22 K=1,N        
   22   ZB(I,J)=ZB(I,J)+ZA1(K,I)*ZA1(K,J)             
   21   CONTINUE        
        
        DO 23 K=1,N        
        
   23   ZS(I)=ZS(I)+ZP1(K)*ZA1(K,I)        
        
   20   CONTINUE        
        
       
c ################        
c ZC: ZB INVERSION        
c ################        
        DO 34 K=1,L        
        
        X=ZB(K,K)        

        DO 31 J=1,L        
        ZB(K,J)=ZB(K,J)/X        
        
   31   ZC(K,J)=ZC(K,J)/X        
        
        DO 32 I=1,L        
        
        IF(I.EQ.K) GO TO 32        
        
        Y=ZB(I,K)        
        
        DO 33 J=1,L        
        ZB(I,J)=ZB(I,J)-Y*ZB(K,J)        
        
   33   ZC(I,J)=ZC(I,J)-Y*ZC(K,J)        
        
   32   CONTINUE        
        
   34   CONTINUE        
        
c #############################        
c solutions XS and residuals XR        
c #############################        
        
        DO 40 I=1,L        
        DO 41 J=1,L        
  41    XS(I)=XS(I)+ZS(J)*ZC(I,J)        
  40    CONTINUE        
        XDS=0.D0        

        DO 42 I=1,N        
       
        XR(I)=ZP(I)        
        
        DO 43 J=1,L   
      
  43    XR(I)= XR(I)-XS(J)*ZA(I,J) 
	 
  42    XDS=XDS+(XR(I)*XR(I))*pp(i)        
        
c   CALCUL DU CHI2        
        
        chi2 = xds/dble(n-l-1)        
         
        write(6,*) 'chi2',chi2
c ##########################################################        
c dispersion XDP, precision XDX and correlation coefficients        
c ##########################################################        
        
        xnu=dble(n-l-kt)        
        Xdp2=xds/xnu        
        XDP=DSQRT(XDP2)        
        DO 50 I=1,L        
        XDX(I)=XDP*DSQRT(dabs(ZC(I,I)))        
        DO 51 J=1,L        
        cor(i,j)=dabs(zc(i,j) / dsqrt(dabs(ZC(i,i)*ZC(j,j))))        
  51    continue        
  50    continue        
        
       
c ##################        
c equations removing        
c ##################        
       
	  if (nopt.ne.0) then        
	               
      	  DO 52 I=1,N        
              XT=2.57d0*XDP        
       
	      XRM=DABS(XR(I))       
        
      	  IF(XRM.lT.XT) go to 52        
	        kt  = kt+1        
	        NES = NES+1        
       
	        DO 53 J=1,L        
       	  ZA(I,J)=0.D0        
   53 	  continue        
	        ZP(I)=0.D0        
   52 	  continue        
	        if (kt.ne.0) go to 100        
	         
	  endif       
	       
	  return              
    	  end
	  
c ################################################################       
        SUBROUTINE MCARRE_FCLORA(n1,l1,n,l,ZA,ZP,PP,XS,XR,XDX,COR,chi2,NES,nopt)       
c ################################################################       
c input : n1 : dimension des tableau d'entree figee pour la compilation avec FCLORA
c         l1 : dimension des tableau d'entree figee 
c         n  : nombre d'observations
c         l  : nombre de parametre  
c        ZA  : matrice des relations d'observations
c        ZP  : observations
c        PP  : matrice des poids
c
c sortie XS  : parametres estimes  
c        XR  : residus  
c        XDX : erreurs sur les parametres estimes 
c        COR : matrice de correlation entre les parametres
c        chi2: estimation du chi2
c        Nes : nombre d'equation supprimees
c        nopt: supression d'equation --> 1 / pas de supression --> 0   
     
        implicit DOUBLE PRECISION (a-h,o-z)             
        integer n, l
	 
        DOUBLE PRECISION ZA(n1,l1),ZP(n1),PP(n1),XR(n1),       
     $         XS(l1),XDX(l1),COR(l1,l1)   

        DOUBLE PRECISION ZA1(n,l),ZP1(n),ZS(n),ZB(l,l),ZC(l,l)       


        kt=0.d0        
        nes=0                                                
      
c #############################################        
c matrixs ZA1=ZA t(sqrt(ZP)) and ZP1=sqrt(ZP) L        
c #############################################        
        
 100    DO i=1,n        
        DO j=1,l        

        ZA1(I,J)=ZA(I,J)*DSQRT(PP(I))        
	ENDDO        
        ZP1(I)=ZP(I)*DSQRT(PP(I))        
	ENDDO        
       
        do i=1,l
        ZS(I)=0.D0  
        enddo
       
        DO 24 I=1,l        
        DO 25 J=1,l        
        ZB(I,J)=0.D0        
 25     ZC(I,J)=0.D0        
 24     ZC(I,I)=1.D0        
        
        DO I=1,l       
        XS(I)=0.D0        
        ENDDO        
        
c ############################        
c normal equations : matrix ZB        
c ############################        
        
        DO 20 I=1,L        
        
        DO 21 J=1,L        
        
        DO 22 K=1,N        
   22   ZB(I,J)=ZB(I,J)+ZA1(K,I)*ZA1(K,J)             
   21   CONTINUE        
        
        DO 23 K=1,N        
        
   23   ZS(I)=ZS(I)+ZP1(K)*ZA1(K,I)        
        
   20   CONTINUE        
        
       
c ################        
c ZC: ZB INVERSION        
c ################        
        DO 34 K=1,L        
        
        X=ZB(K,K)        

        DO 31 J=1,L        
        ZB(K,J)=ZB(K,J)/X        
        
   31   ZC(K,J)=ZC(K,J)/X        
        
        DO 32 I=1,L        
        
        IF(I.EQ.K) GO TO 32        
        
        Y=ZB(I,K)        
        
        DO 33 J=1,L        
        ZB(I,J)=ZB(I,J)-Y*ZB(K,J)        
        
   33   ZC(I,J)=ZC(I,J)-Y*ZC(K,J)        
        
   32   CONTINUE        
        
   34   CONTINUE        
        
c #############################        
c solutions XS and residuals XR        
c #############################        
        
        DO 40 I=1,L        
        DO 41 J=1,L        
  41    XS(I)=XS(I)+ZS(J)*ZC(I,J)        
  40    CONTINUE        
        XDS=0.D0        

        DO 42 I=1,N        
       
        XR(I)=ZP(I)        
        
        DO 43 J=1,L   
      
  43    XR(I)= XR(I)-XS(J)*ZA(I,J) 
	 
  42    XDS=XDS+(XR(I)*XR(I))*pp(i)        
        
c   CALCUL DU CHI2        
        
        chi2=xds/dble(N-L-1)        
         
        write(6,*) 'chi2',chi2
c ##########################################################        
c dispersion XDP, precision XDX and correlation coefficients        
c ##########################################################        
        
        xnu=dble(n-l-kt)        
        Xdp2=xds/xnu        
        XDP=DSQRT(XDP2)        
        DO 50 I=1,L        
        XDX(I)=XDP*DSQRT(dabs(ZC(I,I)))        
        DO 51 J=1,L        
        cor(i,j)=dabs(zc(i,j) / dsqrt(dabs(ZC(i,i)*ZC(j,j))))        
  51    continue        
  50    continue        
        
       
c ##################        
c equations removing        
c ##################        
       
	  if (nopt.ne.0) then        
	               
      	  DO 52 I=1,N        
              XT=2.57d0*XDP        
       
	      XRM=DABS(XR(I))       
        
      	  IF(XRM.lT.XT) go to 52        
	        kt  = kt+1        
	        NES = NES+1        
       
	        DO 53 J=1,L        
       	  ZA(I,J)=0.D0        
   53 	  continue        
	        ZP(I)=0.D0        
   52 	  continue        
	        if (kt.ne.0) go to 100        
	         
	  endif       
	       
	  return              
    	  end            
          
c #############################################################
         SUBROUTINE REGRESSION(n,D,Y,P,a,sig_a,b,sig_b,cor_ab)
c
c objet : calcul du biais et de la pente d'une serie temporelle 
c         par regression lineaire  
c 
c entree    n  : dimension des tableaux D,X,P
c           D  : tableau des dates
c	    Y  : tableau des valeurs
c	    P  : tableau des poids affectant les valeurs
c
c sortie    a     : pente
c           sig_a : erreur sur la pente   
c           b     : biais
c           sig_b : erreur sur le biais
c           cor_ab: correlation entre a et b
c 
c N.B. Cette procedure ne permet pas d'eliminer les valeurs anormales
c  
c Code par Ch. Bizouard en FORTRAN 90 - decembre 2002
c ##############################################################
      
      IMPLICIT NONE
      INTEGER n,i
      DOUBLE PRECISION D(n),Y(n),P(n)    	
      DOUBLE PRECISION a,da,b,db,cor_ab,sig_a,sig_b
      DOUBLE PRECISION s_p,s_pd,s_py,s_pdy,s_pd2,det,sig0
  
      s_p   = 0.d0
      s_pd  = 0.d0 
      s_py  = 0.d0
      s_pdy = 0.d0
      s_pd2 = 0.d0
      
                 
      DO i=1,n
	      s_p  = s_p  + P(i)
	      s_pd = s_pd + P(i) * D(i)
	      s_py = s_py + P(i) * Y(i)
 	      s_pdy= s_pdy+ P(i) * D(i) * Y(i)
	      s_pd2= s_pd2+ P(i) * D(i)**2  
      ENDDO
      
      det = s_p * s_pd2 - s_pd**2
      
      a = ( s_p * s_pdy - s_pd * s_py)   / det  
      b = (-s_pd * s_pdy + s_pd2 * s_py) / det
      
C Calcul des erreurs formelles et de la correlation     
      
      ! sig0 : ecart type par rapport a : Y = a D + b            
      sig0 = 0.d0
 
      DO i = 1,n
	      sig0 = sig0 + P(i) * (Y(i) - a * D(i) - b)**2       
      ENDDO
      
      sig0 = dsqrt(sig0/dble(n-2))
       
      sig_a = sig0 * dsqrt(dabs(s_p  /det))
      sig_b = sig0 * dsqrt(dabs(s_pd2/det)) 
      
      cor_ab= dabs(sig0**2 * s_pd/det)/(sig_a * sig_b)
      
      RETURN  
      END

c ###################################################################################                                                                
         SUBROUTINE MOCACH(NEQ,NINC)
c neq .. nombre d'observation
c ninc ....nombre de parametres
c        mocach pour 10000 mesures 

         DOUBLE PRECISION XX,Y
	 
	 PARAMETER (ne=50) 
         DOUBLE PRECISION W(50,ne),U(ne,ne),B(ne) 
	 
	 INTEGER NEQ,NINC
	 
         DOUBLE PRECISION T, X, SIGMA, R, SIG, A, CC     
         COMMON/CACH/T(10000,50),X(10000),SIGMA,R(10000),SIG(50),A(50),CC(50,50)
c entrees :
C     T(NEQ,NINC)....MATRICE DES EQUATIONS D'OBSERVATION
C     X .............VECTEUR DES OBSERVATIONS
C     B(NINC)........SECONDS MEMBRES DE CES EQUATIONS
C     W(NEQ,NEQ).... MATRICE NORMALE
C     U(NEQ,NEQ).... MATRICE INVERSE
C     CC ........... MATRICE DE CORRELATION

C
C sorties
C     A .............VECTEUR SOLUTION
C     R ............ RESIDUS
C     SIG..........  INCERTITUDES SUR LA SOLUTION
C     SIGMA .......  RESIDU QUADRATIQUE MOYEN
C
              
       DO i=1,ninc
             B(I)=0.d0
       ENDDO       

C             EQUATIONS NORMALES
C
       DO  i=1,ninc
             B(i) = 0.d0
       ENDDO
       
       DO i = 1,neq
             DO ij=1,ninc
                  B(ij)=B(ij)+T(i,ij)*X(i)
             ENDDO
       ENDDO

      
       !MOINDRES CARRES
     
       DO i=1,ninc
             DO j=1,ninc
                  U(i,j) = 0.D0
                  W(i,j) = 0.D0
                  U(i,i) = 1.D0
             ENDDO
       ENDDO 
       	     
       DO 12 I=1,NINC
       DO 12 J=1,NINC
       DO 12 K=1,NEQ
   12  W(I,J)=W(I,J)+T(K,I)*T(K,J)
       
       DO k=1,ninc
             XX = W(k,k)
             IF(XX.EQ.0.d0)  GO TO 117
       
             DO j=1,ninc
                  W(k,j) = W(k,j)/XX
                  U(k,j) = U(k,j)/XX
             ENDDO
	     
             DO 15 I=1,NINC
                  IF(I.EQ.K) GO TO 15
                  Y = W(I,K)
       
                  DO 14 J=1,NINC
                       W(i,j) = W(i,j)-Y*W(k,j)
   14                  U(i,j) = U(i,j)-Y*U(k,j)
   15        CONTINUE
       ENDDO
       
       DO 104 ij=1,ninc
             SIG(ij) = 0.d0
               A(ij) = 0.d0
       
       DO 104 LN=1,NINC
  104  A(IJ)=A(IJ)+U(IJ,LN) * B(LN)

       
       ! RESIDUS...DISPERSION...

       SQ=0
       
       DO 105 I=1,NEQ
       CALC=0.d0
       
       DO 106 J=1,NINC
  106  CALC=CALC+A(J)*T(I,J)
       R(I)=X(I)-CALC
  105  SQ=SQ+R(I)*R(I)
       SIGMA=DSQRT(SQ/dble(NEQ-NINC))
       
       DO 107 J=1,NINC
  107  SIG(J)=SIGMA*DSQRT(U(J,J))
       
       DO 110 I=1,NINC
       DO 110 J=1,NINC
 110   CC(I,J)=U(I,J)/DSQRT(U(I,I)*U(J,J))
       GOTO 119
       
  117  WRITE(6,118)  K,((W(I,J),I=1,NINC),J=1,NINC)
  118  FORMAT(//' DIVISEUR NUL DANS MOCACH',I6,/12(12F10.5/))
  119  CONTINUE
       RETURN
       END

c #################################################################################
       SUBROUTINE PENTAN(DAT,P,n,dat0,ninc)
c
c          pentan pour dim mesures
c
c    AJUSTEMENT DE CONSTANTE, DERIVE LINEAIRE, TERMES ANNUEL ET SEMI ANNUEL.
c    DAT0 EST LA DATE origine pour l'ajustement du terme lineaire.
c
c      
c        entrees   DAT .....DATES EN MJD
c                  P ...... POIDS
c                  n .....  NOMBRE DE MESURES
c                  DAT0 ....DATE ORIGINE
C                  ninc ....si ninc = 1  CALCUL BIAIS ET PENTE
c
C        sorties dans le common de MOCACH 
C                 X ....... BIAIS, PENTE, SIN(AN.), COS(AN.), SIN(SEMI-AN.), COS(SEMI-AN.)
C                 EPS...... RESIDU QUADRATIQUE MOYEN
C                 R ....... TABLEAU DES RESIDUS


      DOUBLE PRECISION DAT(n),P(n),DAT0     	
      DOUBLE PRECISION pi
      INTEGER n

      ! common de mocach
      DOUBLE PRECISION A, UT, eps, R, SIG, X, CC      
      COMMON/CACH/A(10000,50),UT(10000),eps,R(10000),SIG(50),X(50),CC(50,50)

      INTEGER ninc
      
      pi = 4.d0 * datan(1.d0)
      
    
      DO j=1,n
      P(j)   = DSQRT(P(j))
   
      A(j,1) = 1.d0 * P(J)
     
      A(j,2) = DSIN( 2.d0 * pi * (DAT(j)-dat0)/365.242 ) * P(j)
      A(j,3) = DCOS( 2.d0 * pi * (DAT(j)-dat0)/365.242 ) * P(j)
      A(j,4) = DSIN( 4.d0 * pi * (DAT(j)-dat0)/365.242 ) * P(j)
      A(j,5) = DCOS( 4.d0 * pi * (DAT(j)-dat0)/365.242 ) * P(j)
      
      A(j,6) = (DAT(j) - dat0) /365.242 * P(J)

C CAS OU NINC EST 1 
      IF(NINC.EQ.1) A(j,2) = A(j,6)

      UT(j) = UT(j) * P(j)
      P(j)  = P(j)**2 
      enddo
    
      MINC=NINC
      
      IF(NINC.EQ.1) MINC=2

      CALL MOCACH(n,MINC)

      SQ = 0.d0
      SP = 0.d0
      
      DO K=1,n
      SQ = SQ + R(K)**2  ! ecart quadratique des residus
      SP = SP + P(K)  
      R(K) = R(K) / A(K,1)  
      enddo
      
      IF(NINC.EQ.1) THEN
        X(6)   = X(2)
        SIG(6) = SIG(2)
        X(2)   = 0.d0
        SIG(2) = 0.d0
      ENDIF
      
      EPS = DSQRT(DBLE(n) * SQ / (SP*dble(n-MINC) ) )
      RETURN
      END

c --------------------------------------------------------------------------
c --------------------------------------------------------------------------	
       SUBROUTINE MEMCOF(X,n,m,PM,COF)
c input   X(1:n) .... tableau de  valeurs
c         m ......... entier > 0 
c
c output  COF(1:m)... m coefficients  de la prediction lineaire
c         PM ........ ecart moyen des donnees a la prediction        
c --------------------------------------------------------------------------	
       
        implicit none
        integer n,m,j,k,i
        double precision X(n),COF(m),WK1(n),WK2(n),WKM(m)
        double precision p,pm,pneum,denom
       
       
        P=0.d0
        DO j=1,n
        	P = P + X(j)**2
	enddo

        PM       = P / dble(n)
        WK1(1)   = X(1)
        WK2(N-1) = X(N)
       
        do j = 2, n-1
        WK1(j)   = X(j)
        WK2(j-1) = X(j)
	enddo
       
        DO K = 1, m
		PNEUM = 0.d0
        	DENOM = 0.d0
	
        	do J=1,N-K
         		PNEUM = PNEUM + WK1(J) * WK2(J)
          		DENOM = DENOM + WK1(J)**2 + WK2(J)**2
		enddo

                COF(K)=2.d0*PNEUM/DENOM
                PM=PM*(1.d0-COF(K)**2)
	
                IF (K.NE.1) THEN        
	        	do I=1,K-1
                		COF(I) = WKM(I) - COF(K) * WKM(K-I)
			enddo        
	        ENDIF
        
	        IF(K.EQ.M) RETURN
        
		do I=1,K
          		WKM(I)=COF(I)
		enddo
        
		do J=1,N-K-1
         	 	WK1(J) = WK1(J)   - WKM(K) * WK2(J)
         		WK2(J) = WK2(J+1) - WKM(K) * WK1(J+1)
		enddo

	enddo
   
        END

c --------------------------------------------------------------------------
c --------------------------------------------------------------------------	
	 subroutine PACF(X,AC,PAC,n,m)
c
c 	input  X(1:n) .... signal
c              n ......... nombre de valeurs
c 
c	output AC(1:n) ... coefficient autocorélation 
c	       PAC(1:n) .. coefficients autocorélation partiels
c
c --------------------------------------------------------------------------
	
	 implicit none
	 integer i,j,n,m,h,k
	 double precision X(n),AC(0:m),PAC(m),r(m,m)	
	 double precision s1,s2,AC0,moy
	 
	 ! calcul de la moyenne
	 moy =0.d0
	 
	 do i=1,n
	 	moy =moy + X(i)
	 enddo
	 moy =moy/dble(n)
	 	
	 do h=0,m
	 	AC(h)=0.d0
	 	do i=h+1,n
	 		AC(h) = AC(h) + (X(i)-moy) * (X(i-h)-moy) 
	 	enddo
	 enddo
	
	 AC0=AC(0)
	 
	 do h=1,m
	 	AC(h) = AC(h) / AC0 
	 enddo
	 
	 r(1,1)=AC(1)
	 
	 do k=2,m
	 	s1 = AC(k)
	 	s2 = 1.d0
	 	do j=1,k-1
	 		s1 = s1 - R(k-1,j) * AC(k-j)
	 		s2 = s2 - R(k-1,j) * AC(j)
		enddo
	  
	 	r(k,k)= s1/s2
	 
	 	do i=1,k-1
	 		R(k,i) = R(k-1,i) - R(k,k) * R(k-1,k-i)
	 	enddo
	 
	 enddo
	 
	 do k=1,m
	 	PAC(k) = R(k,k)
	 enddo
	 
	 return
	 end	
