!-------------------------------------------------------------------------------------
!					MODULES
!--------------------------------------------------------------------------------------
!    Constantes du programme							1
	module constante
	   integer,parameter			:: nd1=20
	   integer,parameter			:: nd2=2
	   integer,parameter			:: ncm=3	! Nb options calcul
	   integer,parameter			:: nca=22	! Longueur chaine car.
	end module constante
!    ----------------------------------------------------------------------------------
!    Parametres de calcul							2
	module parametre
	   use constante,only: ncm,nca
	   integer				:: n,m,nopt,kx,ky,c,b,nv
	   real(8)				:: beta,wmin,wmax,u,mmu,remp
	   character(len=nca)			:: ficpara='lsvd.dat'
	   character(len=nca)			:: ficin,ficout,ficgr
	   character(len=nca)			:: ficcov
	   character(len=1),dimension(ncm)	:: ynmes
	   integer,dimension(2,20)		:: vect
	   real(8),dimension(20)		:: alph
	end module parametre
!--------------------------------------------------------------------------------------
!				   PROGRAMME PRINCIPAL
!--------------------------------------------------------------------------------------
	Program bryanmem
!--------------------------------------------------------------------------------------
! Nom	     :	meshmem	[MAIN]							3
!	
! Version    : 	(10 octobre 1998)
!
! Fonction   :	Methode de l'entropie maximum selon l'algorithme
!		de Bryan
!
! Auteur(s)  : 	David Poulin
!
! Revision   :  
!
! Adresse    : 	Centre de recherche en physique du solide (CRPS)
!		Departement de physique, Universite de Sherbrooke
!		Sherbrooke, Quebec, Canada, J1K 2R1
!
! References :  [1] Bryan, in Maximum entropy and Bayesian methods, Fougere, 1990.
!		[2] D. Poulin, www.physique.usherb.ca/~tremblay/david
!		[3] Jarrell et Gubernatis, Physics Reports 269 (1996)
!
! Notes: 	- Version avec le moyennage sur le parametre de regularisation (alpha)
!		- Traitement des donnees et du kernel avec la matrice de covariance
!		  (optionnel)
!		- Les donnees qui s'etalent de 0 a TRAN-DTAU (TRAN valeurs) sont stockees
!		  dans g(k), k=1..TRAN  (TRAN = n)
!		- Exteriorisation des parametres de calcul (ficpara)
!		- UNIT 2 : Fichier des parametres de calcul
!		  UNIT 3 : Fichier des donnees
!		  UNIT 5 : Fichier des resultats intermediares
!		  UNIT 7 : Fichier des resultats
!		- Options de calcul
!		  YNMES(1) : Matrice de covariance
!                 YNMES(2) : Imposition de la symétrie particule-trou
!		  YNMES(3) : Doubler le nombre de mesures
!
!--------------------------------------------------------------------------------------
	use constante,only: nd1,nd2
	use parametre
	implicit none
!     Variables scalaires 
	real(8)				:: ss0,ss1,ss2             
	real(8)				:: ti,tf,td
	real(8)				:: dw,dtau,mu,chi2,delta
	real(8)				:: mu1,mu2,pi,ep
	integer				:: i,j,k,l,iterv
	integer				:: s,crit,qx,qy,xx,yy
!     Fonction
	character(8)			:: date
	integer				:: mclock
!     Tableaux
	integer,allocatable		:: it(:)
	real(8),allocatable		:: a(:),w(:)
	real(8),allocatable		:: t(:),g(:),err(:),f(:)
	real(8),allocatable		:: ek(:,:),ekt(:,:),v(:,:),sigma(:)
	real(8),allocatable		:: covm(:,:)
!--------------------------------------------------------------------------------------
!
	ti=mclock()
!
!Lecture parametres de calcul
	call lecture_para
	
	do 200 iterv=1,nv
	kx=vect(1,iterv)
	ky=vect(2,iterv)

	
!
!Ouverture des fichiers d'ecriture
	open(unit=7,file=ficout,position='append')
	open(unit=5,file=ficgr,position='append')
	if(iterv==1) then
	   rewind(unit=7)
	   rewind(unit=5)
	end if
	write(5,*)'#-------------------------------------------------------------'
	write(5,*)'#kx=',kx,'	ky=',ky
	write(7,*)'#-------------------------------------------------------------'
	write(7,*)'#kx=',kx,'	ky=',ky
!
!Declaration dynamique des tableaux
	allocate (t(n),it(n),g(n),err(n))
	allocate (a(m),w(m))
	allocate (ek(n,m),ekt(m,n),sigma(n),v(n,n))
	allocate (covm(n,n),f(n))
!
!Lecture donnees
	call lecture_donnees(it,g,err,iterv)
	
	delta=1
	do i=1,n
	   if(err(i)/g(i)<delta) then
	      delta=err(i)/g(i)
	   end if
	end do
!
!Paramètres de calcul
	dw=wmax-wmin
	dw=dw/(m-1)
	dtau=beta/(n)
	write(*,*) 'dtau=',dtau
	
	do i=1,m
	   w(i)=wmin+(i-1)*dw
	end do 
	
	do k=1,n
           t(k)=it(k)*dtau
	end do
   
!Traitement des donnees par la matrice de covariance
	if(ynmes(1)=='y'.or.ynmes(1)=='Y') then
	   call lecture_cova(n,covm)
	   call invmg(n,n,covm)
	end if 

!Calcul des moments
	pi=acos(-1.d0)
	if(c==1) then
	   ep=-2*cos(pi*2*ky/b)
	end if
	if(b==1) then
	   ep=-2*cos(pi*2*kx/c)
	end if
	if(b.ne.1.and.c.ne.1) then
	   ep=-2*cos(pi*2*kx/c)-2*cos(pi*2*ky/b)
	end if
	
	ep=0
	
	
	mu1=ep-mmu-u/2+0.5*remp*u
	mu2=(ep-mmu-u/2)**2+2*u*(ep-mmu-u/2)*remp/2 &
	   &+u*u*remp/2
		
!noyau d'interraction
	!if(ep-mmu>-1d-6.and.ep-mmu<1d-6) then
	!   do i=1,m
	!      do k=1,n
	!         ek(k,i)=cosh((beta/2-t(k))*w(i))/(2*cosh(beta*w(i)/2))&
	!            &/(1+dexp((-u+dabs(w(i)))*beta/12))
	!         ekt(i,k)=ek(k,i)
	!      end do
	!   end do
	!else
	   do i=1,m
	      do k=1,n
	         ek(k,i)=dexp(-w(i)*t(k))/(1.d0+dexp(-beta*w(i)))&
	            &/(1+dexp((-u+dabs(w(i)))*beta/12)) 
	         ekt(i,k)=ek(k,i)
	      end do
	   end do
	!end if

!Diagonalisation SVD de la transpose du noyau
	call svdcmp(ekt,m,n,m,n,sigma,v)
	s=n
	do i=n,1,-1
	   if(sigma(i)/sigma(1)<delta) then
	      s=i-1 
	      sigma(i)=0.d0
	   end if
	end do
	
	do i=1,m
	  a(i)=0.d0
	  do l=1,s
	    do k=1,n
	      a(i)=a(i)+ekt(i,l)*v(k,l)*g(k)/sigma(l)
	    end do
	  end do
	end do
	
	if(ynmes(1)=='y'.or.ynmes(1)=='Y') then
	   call chicov(n,chi2,g,ek,a,m,covm,dw,f)
	else
	   call chi(n,chi2,g,ek,a,m,err,dw,f)	
	end if
	
	write(5,35) mu1,mu2
	
	a=a/dw

	do i=1,m
	   ss0=ss0+a(i)
	   ss1=ss1+a(i)*w(i)
	   ss2=ss2+a(i)*w(i)*w(i)
	end do
	ss0=ss0*dw
	ss1=ss1*dw
	ss2=ss2*dw
	

	write(5,*) 'ordre 0:', ss0,' 1'
	write(5,*) 'ordre 1:', ss1,' ',mu1
	write(5,*) 'ordre 2:', ss2,' ',mu2
	write(5,*) 's=',s
	
	call sortie(w,a,f,chi2,n)
	if(ynmes(2)=='m'.or.ynmes(2)=='Y') then
	   QX=MOD(KX+C/2,C)
	   QY=MOD(KY+B/2,B)
	   IF(QX==0) QX=C
	   IF(QY==0) QY=B
	   IF(C==1) QX=1
	   IF(B==1) QY=1
	   write(7,*)'#-------------------------------------------------------------'
	   write(7,*)'#kx=',qx,'	ky=',qy
	   xx=c-qx
	   yy=b-qy
	   do i=1,m
	     write(7,44) w(i),a(m+1-i)-xx*5.d-1
	   end do
	end if

	tf=mclock()
	td=(tf-ti)/6000
	write(*,*) td,'Temps de calcul (min CPU)'
!
	close(unit=5)
	close(unit=7)
	close(unit=9)
!
!	deallocation dynamique des tableaux
	
	deallocate (t,it,g,err)
	deallocate (a,w)
	deallocate (ek,ekt,sigma,v)
	deallocate (covm,f)
	
200	end do
!
33	format(3(1x,f10.5),3x,i3.3)
34	format(3(1x,f10.5),3x,i3.3)
44	format(2(1x,f10.5))
35	format(34x,2(f10.5,1x))
	end program bryanmem
!---------------------------------------------------------------------------
!				Routines
!---------------------------------------------------------------------------
!	Nom:		CHI						4
!	Fonction: 	Calcul la fonction de Green d'essai et le chi2 
!	Auteur:		David Poulin
!---------------------------------------------------------------------------
	SUBROUTINE CHI(N,CHI2,G,EK,A,M,ERR,DW,F)
	REAL(8)				:: CHI2,DW
	INTEGER				:: N,M
	REAL(8),DIMENSION(N)		:: G,ERR,F
	REAL(8),DIMENSION(M)		:: A
	REAL(8),DIMENSION(N,M)		:: EK
	INTEGER				:: I,L
	
	CHI2=0.D0
	DO I=1,N
	   F(I)=0.D0
	   DO L=1,M
	      F(I)=F(I)+EK(I,L)*A(L)
	   END DO
	   F(I)=F(I)*DW
	END DO
	DO I=1,N
	   CHI2=CHI2+((F(I)-G(I))/ERR(I))**2
	END DO
	END SUBROUTINE CHI
!---------------------------------------------------------------------------
!	Nom:		CHICOV						5
!	Fonction: 	Calcul la fonction de Green d'essai et le chi2 
!	Auteur:		David Poulin
!---------------------------------------------------------------------------
	SUBROUTINE CHICOV(N,CHI2,G,EK,A,M,COV,DW,F)
	REAL(8)				:: CHI2,DW
	INTEGER				:: N,M
	REAL(8),DIMENSION(N)		:: G,F
	REAL(8),DIMENSION(M)		:: A
	REAL(8),DIMENSION(N,M)		:: EK
	REAL(8),DIMENSION(N,N)		:: COV
	INTEGER				:: I,L
	
	CHI2=0.D0
	DO I=1,N
	   F(I)=0.D0
	   DO L=1,M
	      F(I)=F(I)+EK(I,L)*A(L)
	   END DO
	   F(I)=F(I)*DW
	END DO
	DO I=1,N
	   DO L=1,N
	      CHI2=CHI2+((F(I)-G(I))*(F(L)-G(L)))*COV(I,L)
	   END DO
	END DO
	END SUBROUTINE CHICOV
!---------------------------------------------------------------------------
!	Nom:		SVDCMP						6
!	Fonction:	Single value decomposition
!	Auteur:		Numerical Recipes
!---------------------------------------------------------------------------
     SUBROUTINE SVDCMP(A,M,N,MP,NP,W,V)
      implicit real*8 (a-h,o-z)
      PARAMETER (NMAX=200)
      DIMENSION A(MP,NP),W(NP),V(NP,NP),RV1(NMAX)
      G=0.0
      SCALE=0.0
      ANORM=0.0
      DO 25 I=1,N
        L=I+1
        RV1(I)=SCALE*G
        G=0.0
        S=0.0
        SCALE=0.0
        IF (I.LE.M) THEN
          DO 11 K=I,M
            SCALE=SCALE+ABS(A(K,I))
11        CONTINUE
          IF (SCALE.NE.0.0) THEN
            DO 12 K=I,M
              A(K,I)=A(K,I)/SCALE
              S=S+A(K,I)*A(K,I)
12          CONTINUE
            F=A(I,I)
            G=-SIGN(SQRT(S),F)
            H=F*G-S
            A(I,I)=F-G
            IF (I.NE.N) THEN
              DO 15 J=L,N
                S=0.0
                DO 13 K=I,M
                  S=S+A(K,I)*A(K,J)
13              CONTINUE
                F=S/H
                DO 14 K=I,M
                  A(K,J)=A(K,J)+F*A(K,I)
14              CONTINUE
15            CONTINUE
            ENDIF
            DO 16 K= I,M
              A(K,I)=SCALE*A(K,I)
16          CONTINUE
          ENDIF
        ENDIF
        W(I)=SCALE *G
        G=0.0
        S=0.0
        SCALE=0.0
        IF ((I.LE.M).AND.(I.NE.N)) THEN
          DO 17 K=L,N
            SCALE=SCALE+ABS(A(I,K))
17        CONTINUE
          IF (SCALE.NE.0.0) THEN
            DO 18 K=L,N
              A(I,K)=A(I,K)/SCALE
              S=S+A(I,K)*A(I,K)
18          CONTINUE
            F=A(I,L)
            G=-SIGN(SQRT(S),F)
            H=F*G-S
            A(I,L)=F-G
            DO 19 K=L,N
              RV1(K)=A(I,K)/H
19          CONTINUE
            IF (I.NE.M) THEN
              DO 23 J=L,M
                S=0.0
                DO 21 K=L,N
                  S=S+A(J,K)*A(I,K)
21              CONTINUE
                DO 22 K=L,N
                  A(J,K)=A(J,K)+S*RV1(K)
22              CONTINUE
23            CONTINUE
            ENDIF
            DO 24 K=L,N
              A(I,K)=SCALE*A(I,K)
24          CONTINUE
          ENDIF
        ENDIF
        ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I))))
25    CONTINUE
      DO 32 I=N,1,-1
        IF (I.LT.N) THEN
          IF (G.NE.0.0) THEN
            DO 26 J=L,N
              V(J,I)=(A(I,J)/A(I,L))/G
26          CONTINUE
            DO 29 J=L,N
              S=0.0
              DO 27 K=L,N
                S=S+A(I,K)*V(K,J)
27            CONTINUE
              DO 28 K=L,N
                V(K,J)=V(K,J)+S*V(K,I)
28            CONTINUE
29          CONTINUE
          ENDIF
          DO 31 J=L,N
            V(I,J)=0.0
            V(J,I)=0.0
31        CONTINUE
        ENDIF
        V(I,I)=1.0
        G=RV1(I)
        L=I
32    CONTINUE
      DO 39 I=N,1,-1
        L=I+1
        G=W(I)
        IF (I.LT.N) THEN
          DO 33 J=L,N
            A(I,J)=0.0
33        CONTINUE
        ENDIF
        IF (G.NE.0.0) THEN
          G=1.0/G
          IF (I.NE.N) THEN
            DO 36 J=L,N
              S=0.0
              DO 34 K=L,M
                S=S+A(K,I)*A(K,J)
34            CONTINUE
              F=(S/A(I,I))*G
              DO 35 K=I,M
                A(K,J)=A(K,J)+F*A(K,I)
35            CONTINUE
36          CONTINUE
          ENDIF
          DO 37 J=I,M
            A(J,I)=A(J,I)*G
37        CONTINUE
        ELSE
          DO 38 J= I,M
            A(J,I)=0.0
38        CONTINUE
        ENDIF
        A(I,I)=A(I,I)+1.0
39    CONTINUE
      DO 49 K=N,1,-1
        DO 48 ITS=1,30
          DO 41 L=K,1,-1
            NM=L-1
            IF ((ABS(RV1(L))+ANORM).EQ.ANORM)  GO TO 2
            IF ((ABS(W(NM))+ANORM).EQ.ANORM)  GO TO 1
41        CONTINUE
1         C=0.0
          S=1.0
          DO 43 I=L,K
            F=S*RV1(I)
            IF ((ABS(F)+ANORM).NE.ANORM) THEN
              G=W(I)
              H=SQRT(F*F+G*G)
              W(I)=H
              H=1.0/H
              C= (G*H)
              S=-(F*H)
              DO 42 J=1,M
                Y=A(J,NM)
                Z=A(J,I)
                A(J,NM)=(Y*C)+(Z*S)
                A(J,I)=-(Y*S)+(Z*C)
42            CONTINUE
            ENDIF
43        CONTINUE
2         Z=W(K)
          IF (L.EQ.K) THEN
            IF (Z.LT.0.0) THEN
              W(K)=-Z
              DO 44 J=1,N
                V(J,K)=-V(J,K)
44            CONTINUE
            ENDIF
            GO TO 3
          ENDIF
          IF (ITS.EQ.30) PAUSE 'No convergence in 30 iterations SVD'
          X=W(L)
          NM=K-1
          Y=W(NM)
          G=RV1(NM)
          H=RV1(K)
          F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
          G=SQRT(F*F+1.0)
          F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X
          C=1.0
          S=1.0
          DO 47 J=L,NM
            I=J+1
            G=RV1(I)
            Y=W(I)
            H=S*G
            G=C*G
            Z=SQRT(F*F+H*H)
            RV1(J)=Z
            C=F/Z
            S=H/Z
            F= (X*C)+(G*S)
            G=-(X*S)+(G*C)
            H=Y*S
            Y=Y*C
            DO 45 NM=1,N
              X=V(NM,J)
              Z=V(NM,I)
              V(NM,J)= (X*C)+(Z*S)
              V(NM,I)=-(X*S)+(Z*C)
45          CONTINUE
            Z=SQRT(F*F+H*H)
            W(J)=Z
            IF (Z.NE.0.0) THEN
              Z=1.0/Z
              C=F*Z
              S=H*Z
            ENDIF
            F= (C*G)+(S*Y)
            X=-(S*G)+(C*Y)
            DO 46 NM=1,M
              Y=A(NM,J)
              Z=A(NM,I)
              A(NM,J)= (Y*C)+(Z*S)
              A(NM,I)=-(Y*S)+(Z*C)
46          CONTINUE
47        CONTINUE
          RV1(L)=0.0
          RV1(K)=F
          W(K)=X
48      CONTINUE
3       CONTINUE
49    CONTINUE
      RETURN
      END
!-------------------------------------------------------------------------------------
! Nom		: lecture_para							7
! Fonction	: Lecture des parametres de calcul
! Auteur(s)	: S. Allen
! Revision	: 10 Octobre 1998 D. Poulin
! Notes		:
!--------------------------------------------------------------------------------------

	SUBROUTINE LECTURE_PARA
        USE CONSTANTE
        USE PARAMETRE
	IMPLICIT NONE
	INTEGER				:: I,I1,I2,I3,I4,I5,I6,I7,I8,I9
	CHARACTER(LEN=NCA)		:: TP1,TP2,TP3,TP4,TP5,TP6,TP7
	CHARACTER(LEN=NCA)		:: PREFIX='svd',NANUM,NANUM2,BUFFER
	CHARACTER(LEN=10)		:: NUMTOASCII
        OPEN(UNIT=4,FILE=FICPARA,STATUS='OLD',ACCESS='SEQUENTIAL',POSITION=&
            &'REWIND',ERR=110)
        READ(4,*) TP2
        READ(4,*) BETA
        READ(4,*) REMP
        READ(4,*) MMU
        READ(4,*) N
        READ(4,*) U
	READ(4,*) C
	READ(4,*) B
	READ(4,*) M
	READ(4,*) WMIN
	READ(4,*) WMAX
	READ(4,*) NV
	READ(4,*) NOPT
	DO I=1,NOPT
	    READ(4,*) YNMES(I)
	END DO
	READ(4,*) BUFFER
	READ(4,*) BUFFER
	READ(4,*) BUFFER
	DO I=1,NV
	   READ(4,*) VECT(1,I),VECT(2,I)
	END DO
        CLOSE(UNIT=4)
!
	TP1=NUMTOASCII(INT(ABS(U)))
	TP3=NUMTOASCII(C)
	TP4=NUMTOASCII(B)
	TP5=NUMTOASCII(INT(BETA))
	TP6=NUMTOASCII(KX)
	TP7=NUMTOASCII(KY)
	I1=INDEX(TP1,' ')
	I2=INDEX(TP2,' ')
	I3=INDEX(TP3,' ')
	I4=INDEX(TP4,' ')
	I5=INDEX(TP5,' ')
	I6=INDEX(PREFIX,' ')
	I8=INDEX(TP6,' ')-1
	I9=INDEX(TP7,' ')-1
	NANUM=TP2(1:I2-1)//'.'//PREFIX(1:I6-1)//TP1(1:I1-1)//&
	    &TP3(1:I3-1)//TP4(1:I4-1)//TP5(1:I5-1)//' '
	NANUM2=TP2(1:I2-1)//'.'//'sae'//TP1(1:I1-1)//&
	    &TP3(1:I3-1)//TP4(1:I4-1)//TP5(1:I5-1)//' '
	I7=INDEX(NANUM,' ')-1
	FICIN=NANUM2(1:I7)//'gi.dat'
	FICOUT=NANUM(1:I7)//'ps'//'.dat'
	FICCOV=NANUM2(1:I7)//'cv.dat'
	FICGR=NANUM(1:I7)//'inv.dat'
!
        RETURN
110     WRITE(*,*) 'Ne peut ouvrir le fichier FICIN'; STOP
        END SUBROUTINE LECTURE_PARA
!
!--------------------------------------------------------------------------------------
! Nom		: lecture_donnees	[sub]					8
! Fonction	: Lecture des donnees
! Auteur(s)	: D. Poulin
! Notes		:
!--------------------------------------------------------------------------------------

	SUBROUTINE LECTURE_DONNEES(IT,G,ERR,ITERV)
	USE CONSTANTE
	USE PARAMETRE
	IMPLICIT NONE
	INTEGER					:: QX,QY,KP,PX,PY,RX,RY,iterv
	INTEGER,INTENT(OUT),DIMENSION(N+1)	:: IT
	REAL(8),DIMENSION(N+1,C,B)		:: GR
	REAL(8),DIMENSION(N)			:: G,ERR
	REAL(8),DIMENSION(N+1,C,B)		:: ERRT
	REAL(8),DIMENSION(N+1)			:: GSYM
	CHARACTER(LEN=NCA)			:: LIGNE
	CHARACTER(LEN=1),DIMENSION(NCM)		:: YNMS
	CHARACTER				:: BUFFER
	INTEGER					:: NUMY,L,I,J,K,NO,MESC,BL,P
	REAL(8)					:: TEMP,SUM,SUM2,GO

	OPEN(UNIT=11,FILE=FICIN,STATUS='OLD',ACCESS='SEQUENTIAL',POSITION=&
	    &'REWIND',ERR=100)


	READ(11,*) LIGNE
	READ(11,*) NO
	READ(11,*) MESC
	DO I=1,6
	    READ(11,*) BUFFER
	END DO
	
        DO J=1,C
            DO I=1,B
                DO K=1,N+1
                     READ(11,*) NUMY,NUMY,IT(K),GR(K,J,I),ERRT(K,J,I)
                END DO
            END DO
        END DO
        
	IF(YNMES(3)=='Y'.OR.YNMES(3)=='y') THEN
        IF(YNMES(2)=='Y'.OR.YNMES(2)=='y') THEN
	   QX=MOD(KX+C/2,C)
	   QY=MOD(KY+B/2,B)
	   IF(QX==0) QX=C
	   IF(QY==0) QY=B
	   IF(C==1) QX=1
	   IF(B==1) QY=1
	   PX=C-KX
	   PY=B-KY
	   IF(PX==0) PX=C
	   IF(PY==0) PY=B
	   RX=C-QX
	   RY=B-QY
	   IF(RX==0) RX=C
	   IF(RY==0) RY=B
	   DO K=1,N
	      KP=N+2-K
	      G(K)=(GR(K,KX,KY)+GR(K,PX,PY)&
	         +GR(KP,QX,QY)+GR(KP,RX,RY))/4
	      SUM=(GR(K,KX,KY)**2+GR(K,PX,PY)**2 &
	         & +GR(KP,QX,QY)**2+GR(KP,RX,RY)**2)/4
	      SUM2=(ERRT(K,KX,KY)**2+ERRT(K,PX,PY)**2 &
	         & +ERRT(KP,QX,QY)**2+ERRT(KP,RX,RY)**2)/4
	      ERR(K)=DSQRT(SUM+SUM2-G(K)**2)
	   END DO
	ELSE
	   QX=C-KX
           QY=B-KY
           IF(QX==0) QX=C
           IF(QY==0) QY=B
           DO K=1,N
              G(K)=(GR(K,KX,KY)+GR(K,QX,QY))/2
              ERR(K)=(ERRT(K,KX,KY)**2+ERRT(K,QX,QY)**2 &
                 &+GR(K,KX,KY)**2+GR(K,QX,QY)**2)/2-G(K)**2
              ERR(K)=SQRT(ERR(K))
           END DO
        END IF
        ELSE
        IF(YNMES(2)=='Y'.OR.YNMES(2)=='y') THEN
	   QX=MOD(KX+C/2,C)
	   QY=MOD(KY+B/2,B)
	   IF(QX==0) QX=C
	   IF(QY==0) QY=B
	   IF(C==1) QX=1
	   IF(B==1) QY=1
	   DO K=1,N
	      KP=N+2-K
	      G(K)=(GR(K,KX,KY)+GR(KP,QX,QY))/2
	      SUM=(GR(K,KX,KY)**2+GR(KP,QX,QY)**2)/2
	      SUM2=(ERRT(K,KX,KY)**2+ERRT(KP,QX,QY)**2)/2
	      ERR(K)=DSQRT(SUM+SUM2-G(K)**2)
	   END DO
	ELSE
           DO K=1,N
              G(K)=GR(K,KX,KY)
              ERR(K)=ERRT(K,KX,KY)
           END DO
        END IF
        END IF
              
!
	IF(ITERV==1) THEN
	   WRITE(7,*) '#',LIGNE
	   WRITE(7,*) '#',NO,'        No. simulation'
	   WRITE(7,*) '#',MESC,'      Nb mesures effectuees'
	   WRITE(7,*) '# ',YNMES(1),'	MATRICE DE COVARIANCE'
	   WRITE(7,*) '# ',YNMES(2),'	SYMÉTRIE DE G'
	   WRITE(7,*) '# ',YNMES(3),'	DOUBLER LE NOMBRE DE MESURES'
	   WRITE(7,*) '# '
	END IF
	CLOSE(UNIT=11)
	RETURN
!    Messages d'erreurs
100	WRITE(*,*) 'Erreur : Ne peut ouvrir le fichier ',FICIN; STOP
	END SUBROUTINE LECTURE_DONNEES
!--------------------------------------------------------------------------------------
! Nom           : lecture_cova          [sub]                                   9
! Fonction      : Lecture de la matrice de covariance
! Auteur(s)     : H.Touchette
! Revision      : 1998 (Poulin)
! Notes         : - UNIT 3
!--------------------------------------------------------------------------------------
        subroutine lecture_cova(n,covm)
        use parametre,only: ficcov,kx,ky,c,b
        implicit none
        integer                                 :: i,j,k1,k2,l1,l2,p,lim
        integer,intent(in)                      :: n
        real(8),dimension (n,n),intent(out)     :: covm
        character                               :: buffer
!
        open(unit=3,file=ficcov)
        rewind(unit=3)
        do i=1,5
           read(3,*) buffer
        end do
        lim=(kx-1)*b
        do p=1,lim
          do i=1,n
            do j=1,n
              read(3,*) buffer
            end do
          end do
        end do
        do p=1,ky-1
          do i=1,n
            do j=1,n
              read(3,*) buffer
            end do
          end do
        end do
        do i=1,n
           do j=1,n
              read(3,*)  k1,k2,l1,l2,covm(i,j)
           end do
        end do
        end subroutine lecture_cova
!--------------------------------------------------------------------------------------
! Nom		: sortie		[sub]					10
! Fonction	: Impression des resultats
! Auteur(s)	: D. Poulin
! Revision	: 
! Notes		:
!--------------------------------------------------------------------------------------
	subroutine sortie(w,a,f,chi2,n)
	use parametre,only: m,kx,ky,b,c
	implicit none
	integer					:: i,j,xx,yy,n
	real(8),dimension (m),intent(in)	:: w,a
	real(8),dimension (n),intent(in)	:: f
	real(8)					:: chi2
!
	xx=mod(kx,c)
	yy=mod(ky,b)
	do i=1,m
	  write(7,44) w(i),a(i)-xx*5.d-1
	end do
	write(5,*) 'chi2=',chi2
	do i=1,n
	   write(5,45) i-1,f(i)
	end do
44	format(3(1x,f10.5))
45	format(i3,f10.5)
	end subroutine sortie
!--------------------------------------------------------------------------------------
! Nom		: funcd 		[sub]					11
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
	subroutine funcd(m,d,alpha,entr,falpha,df)
	implicit real*8(a-h,o-z)
	dimension d(m)
	falpha=2.d0*alpha*entr
	df=2.d0*entr
	do i=1,m
	   falpha=falpha-d(i)/(alpha+d(i))
	   df=df+d(i)/(alpha+d(i))**2
	end do
	end subroutine funcd
!--------------------------------------------------------------------------------------
! Nom		: TQLI			[sub]					12
! Fonction	:
! Auteur(s)	: Numerical recipes
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      SUBROUTINE TQLI(D,E,N,NP,Z)
      implicit real*8(a-h,o-z)
      DIMENSION D(NP),E(NP),Z(NP,NP)
      DIMENSION WKSP(Np),IWKSP(Np)
      IF (N.GT.1) THEN
      DO 11 I=2,N
      E(I-1)=E(I)
11    CONTINUE
      E(N)=0.d0
      DO 15 L=1,N
      ITER=0
1     DO 12 M=L,N-1
      DD=dABS(D(M))+dABS(D(M+1))
      IF (dABS(E(M))+DD.EQ.DD) GO TO 2
12    CONTINUE
      M=N
2     IF(M.NE.L)THEN
      IF(ITER.EQ.1500)PAUSE 'too many iterations TQLI'
      ITER=ITER+1
      G=(D(L+1)-D(L))/(2.d0*E(L))
      R=dSQRT(G**2+1.d0)
      G=D(M)-D(L)+E(L)/(G+dSIGN(R,G))
      S=1.d0
      C=1.d0
      P=0.d0
      DO 14 I=M-1,L,-1
      F=S*E(I)
      B=C*E(I)
      IF(dABS(F).GE.dABS(G))THEN
      C=G/F
      R=dSQRT(C**2+1.d0)
      E(I+1)=F*R
      S=1.d0/R
      C=C*S
      ELSE
      S=F/G
      R=dSQRT(S**2+1.d0)
      E(I+1)=G*R
      C=1.d0/R
      S=S*C
      ENDIF
      G=D(I+1)-P
      R=(D(I)-G)*S+2.d0*C*B
      P=S*R
      D(I+1)=G+P
      G=C*R-B
      DO 13 K=1,N
      F=Z(K,I+1)
      Z(K,I+1)=S*Z(K,I)+C*F
      Z(K,I)=C*Z(K,I)-S*F
13    CONTINUE
14    CONTINUE
      D(L)=D(L)-P
      E(L)=G
      E(M)=0.d0
      GO TO 1
      ENDIF
15    CONTINUE
      ENDIF
      call SORT3(N,np,D,Z,WKSP,IWKSP)
      RETURN
      END
!--------------------------------------------------------------------------------------
! Nom		: TRED			[sub]					13
! Fonction	:
! Auteur(s)	:
! Revision	:Numerical recipes
! Notes		:
!--------------------------------------------------------------------------------------
      SUBROUTINE TRED(A,N,NP,D,E)
      implicit real*8(a-h,o-z)
      DIMENSION A(NP,NP),D(NP),E(NP)
      IF(N.GT.1)THEN
      DO 18 I=N,2,-1
      L=I-1
      H=0.d0
      SCALE=0.d0
      IF(L.GT.1)THEN
      DO 11 K=1,L
      SCALE=SCALE+dABS(A(I,K))
11    CONTINUE
      IF(SCALE.EQ.0.d0)THEN
      E(I)=A(I,L)
      ELSE
      DO 12 K=1,L
      A(I,K)=A(I,K)/SCALE
      H=H+A(I,K)**2
12    CONTINUE
      F=A(I,L)
      G=-dSIGN(dSQRT(H),F)
      E(I)=SCALE*G
      H=H-F*G
      A(I,L)=F-G
      F=0.d0
      DO 15 J=1,L
      A(J,I)=A(I,J)/H
      G=0.d0
      DO 13 K=1,J
      G=G+A(J,K)*A(I,K)
13    CONTINUE
      IF(L.GT.J)THEN
      DO 14 K=J+1,L
      G=G+A(K,J)*A(I,K)
14    CONTINUE
      ENDIF
      E(J)=G/H
      F=F+E(J)*A(I,J)
15    CONTINUE
      HH=F/(H+H)
      DO 17 J=1,L
      F=A(I,J)
      G=E(J)-HH*F
      E(J)=G
      DO 16 K=1,J
      A(J,K)=A(J,K)-F*E(K)-G*A(I,K)
16    CONTINUE
17    CONTINUE
      ENDIF
      ELSE
      E(I)=A(I,L)
      ENDIF
      D(I)=H
18    CONTINUE
      ENDIF
      D(1)=0.d0
      E(1)=0.d0
      DO 23 I=1,N
      L=I-1
      IF(D(I).NE.0.d0)THEN
      DO 21 J=1,L
      G=0.d0
      DO 19 K=1,L
      G=G+A(I,K)*A(K,J)
19    CONTINUE
      DO 20 K=1,L
      A(K,J)=A(K,J)-G*A(K,I)
20    CONTINUE
21    CONTINUE
      ENDIF
      D(I)=A(I,I)
      A(I,I)=1.d0
      IF(L.GE.1)THEN
      DO 22 J=1,L
      A(I,J)=0.d0
      A(J,I)=0.d0
22    CONTINUE
      ENDIF
23    CONTINUE
      RETURN
      END
!--------------------------------------------------------------------------------------
! Nom		: SORT3			[sub]					14
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      SUBROUTINE SORT3(N,np,RA,RB,WKSP,IWKSP)
      implicit real*8(a-h,o-z)
      DIMENSION RA(Np),RB(Np,Np),WKSP(Np),IWKSP(Np)
      CALL INDEXX(N,np,RA,IWKSP)
      DO 11 J=1,N
      WKSP(J)=RA(J)
11    CONTINUE
      DO 12 J=1,N
      RA(J)=WKSP(IWKSP(J))
12    CONTINUE
      do i=1,n
      DO 13 J=1,N
      WKSP(J)=RB(i,J)
13    CONTINUE
      DO 14 J=1,N
      RB(i,J)=WKSP(IWKSP(J))
14    CONTINUE
      end do
      RETURN
      END
!--------------------------------------------------------------------------------------
! Nom		: INDEXX		[sub]					15
! Fonction	:
! Auteur(s)	:
! Revision	:
! Notes		:
!--------------------------------------------------------------------------------------
      SUBROUTINE INDEXX(N,np,ARRIN,INDX)
      implicit real*8(a-h,o-z)
      DIMENSION ARRIN(Np),INDX(Np)
      DO 11 J=1,N
      INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
      IF(L.GT.1)THEN
      L=L-1
      INDXT=INDX(L)
      Q=ARRIN(INDXT)
      ELSE
      INDXT=INDX(IR)
      Q=ARRIN(INDXT)
      INDX(IR)=INDX(1)
      IR=IR-1
      IF(IR.EQ.1)THEN
      INDX(1)=INDXT
      RETURN
      ENDIF
      ENDIF
      I=L
      J=L+L
20    IF(J.LE.IR)THEN
      IF(J.LT.IR)THEN
      IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
      ENDIF
      IF(Q.LT.ARRIN(INDX(J)))THEN
      INDX(I)=INDX(J)
      I=J
      J=J+J
      ELSE
      J=IR+1
      ENDIF
      GO TO 20
      ENDIF
      INDX(I)=INDXT
      GO TO 10
      return
      END
!--------------------------------------------------------------------------------------
! Nom		: INVMG 	[SUB]						16
! Fonction	: Inversion d'une matrice carree
! Auteurs	: H.Touchette
! Revision	: 31 juillet 97 (Touchette)
! Notes		: A -> A^(-1)=AI
!		  Provient de Numerical Recipes, p.40.
! Dependances	: LUDCMP,LUBKSB (externes)
!--------------------------------------------------------------------------------------
	SUBROUTINE INVMG(NP,NM,A)
	IMPLICIT NONE
	REAL(8)					:: D		! Var. locale
	INTEGER					:: I,J		! var. locale
	INTEGER,INTENT(IN)			:: NP,NM
	INTEGER,DIMENSION (NP)			:: INDX   	! Var. locale
	REAL(8),DIMENSION (NP,NP)		:: Y		! Var. locale
	REAL(8),DIMENSION (NP,NP),INTENT(INOUT)	:: A
!
	DO I=1,NM
	   DO J=1,NM
	      Y(I,J)=0.D0
	   END DO
	   Y(I,I)=1.D0
	END DO
	CALL LUDCMP(A,NM,NP,INDX,D)
	DO J=1,NM
	   CALL LUBKSB(A,NM,NP,INDX,Y(1,J))
	END DO
	A=Y
	END SUBROUTINE INVMG
!--------------------------------------------------------------------------------------
! Nom		: LUDCMP 	[SUB]						17
! Fonction	: Decomposition gaussienne A = LU
! Auteurs	:
! Revision	:
! Notes		: Provient entierement de Numerical Recipes
!--------------------------------------------------------------------------------------
      	SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      	IMPLICIT REAL(8) (A-H,O-Z)
      	PARAMETER (NMAX=1000,TINY=1.0E-20)
      	DIMENSION A(NP,NP),INDX(N),VV(NMAX)
      	D=1.D0
      	DO 12 I=1,N
           AAMAX=0.D0
           DO 11 J=1,N
              IF(ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11         CONTINUE
           IF(AAMAX.EQ.0.) PAUSE '--- Erreur : Singular matrix LUDCMP'
           VV(I)=1.D0/AAMAX
12     	CONTINUE
      	DO 19 J=1,N
           IF(J.GT.1) THEN
              DO 14 I=1,J-1
                 SUM=A(I,J)
                 IF(I.GT.1) THEN
                    DO 13 K=1,I-1
                       SUM=SUM-A(I,K)*A(K,J)
13                  CONTINUE
                    A(I,J)=SUM
                 ENDIF
14            CONTINUE
           ENDIF
           AAMAX=0.D0
           DO 16 I=J,N
              SUM=A(I,J)
              IF(J.GT.1) THEN
                 DO 15 K=1,J-1
                    SUM=SUM-A(I,K)*A(K,J)
15               CONTINUE
                 A(I,J)=SUM
              ENDIF
              DUM=VV(I)*ABS(SUM)
              IF(DUM.GE.AAMAX) THEN
                 IMAX=I
                 AAMAX=DUM
              ENDIF
16         CONTINUE
           IF(J.NE.IMAX) THEN
              DO 17 K=1,N
                 DUM=A(IMAX,K)
                 A(IMAX,K)=A(J,K)
                 A(J,K)=DUM
17            CONTINUE
              D=-D
              VV(IMAX)=VV(J)
           ENDIF
           INDX(J)=IMAX
           IF(J.NE.N) THEN
             IF(A(J,J).EQ.0.) A(J,J)=TINY
             DUM=1./A(J,J)
             DO 18 I=J+1,N
                A(I,J)=A(I,J)*DUM
18           CONTINUE
           ENDIF
19    	CONTINUE
      	IF(A(N,N).EQ.0.) A(N,N)=TINY
      	RETURN
      	END
!--------------------------------------------------------------------------------------
! Nom		: LUBKSB 	[SUB]						18
! Fonction	: Backsubstitution (pour l'inversion)
! Auteurs	:
! Revision	:
! Notes		: Provient de Numerical Recipes
!--------------------------------------------------------------------------------------
	SUBROUTINE LUBKSB(A,N,NP,INDX,B)
	INTEGER N,NP,INDX(N)
	REAL(8) A(NP,NP),B(N)
	INTEGER I,II,J,LL
	REAL(8) SUM
	II=0
	DO I=1,N
	   LL=INDX(I)
	   SUM=B(LL)
	   B(LL)=B(I)
	   IF(II.NE.0) THEN
	      DO J=II,I-1
	         SUM=SUM-A(I,J)*B(J)
	      END DO
	   ELSE IF(SUM.NE.0) THEN
	      II=I
	   END IF
	   B(I)=SUM
	END DO
	DO I=N,1,-1
	   SUM=B(I)
	   DO J=I+1,N
	      SUM=SUM-A(I,J)*B(J)
	   END DO
	   B(I)=SUM/A(I,I)
	END DO
	RETURN
	END
!--------------------------------------------------------------------------------------
! Nom		:NUMTOASCII							19
! Fonction	:Transforme en CHARACTER
! Auteur	:C. Foisi
!--------------------------------------------------------------------------------------
	CHARACTER(LEN=10) FUNCTION NUMTOASCII(I)
	IMPLICIT NONE
	CHARACTER(LEN=10)		:: TMP
	INTEGER,INTENT(IN)		:: I
	INTEGER				:: J,K
	J=I
	K=1
	NUMTOASCII(1:10)='          '
	TMP=NUMTOASCII
	DO WHILE (J.GT.0.OR.K.EQ.1)
	    TMP(K:K)=CHAR(ICHAR('0')+MOD(J,10))
	    K=K+1
	    J=J/10
	END DO
	J=INDEX(TMP,' ')
	DO K=1,J-1
	    NUMTOASCII(K:K)=TMP(J-K:J-K)
	END DO
	END FUNCTION NUMTOASCII
!--------------------------------------------------------------------------------------
! Nom           : GGUBFS        [FUNCTION]                                      20
! Fonction      : Generateur de nombres pseudo-aleatoires
! Auteurs       :
! Revision      : 15 decembre 1996
! References    : [1] Numerical Recipes, Fortran 90, p.1142.
! Notes         : Remplace le generateur de P. L'Ecuyer.
!
!--------------------------------------------------------------------------------------
        FUNCTION GGUBFS(DSEED)
        IMPLICIT NONE
        INTEGER,PARAMETER               :: K4B=SELECTED_INT_KIND(9)
        INTEGER(K4B),INTENT(INOUT)      :: DSEED
        REAL(8)                         :: GGUBFS
        INTEGER(K4B),PARAMETER          :: IA=16807,IM=2147483647,IQ=127773,IR=2836
        REAL(8),SAVE                    :: AM
        INTEGER(K4B),SAVE               :: IX=-1,IY=-1,K
!
        IF(DSEED<=0.OR.IY<0) THEN
           AM=NEAREST(1.0,-1.0)/IM
           IY=IOR(IEOR(888889999,ABS(DSEED)),1)
           IX=IEOR(777755555,ABS(DSEED))
           DSEED=ABS(DSEED)+1
        ENDIF
        IX=IEOR(IX,ISHFT(IX,13))
        IX=IEOR(IX,ISHFT(IX,-17))
        IX=IEOR(IX,ISHFT(IX,5))
        K=IY/IQ
        IY=IA*(IY-K*IQ)-IR*K
        IF(IY<0) IY=IY+IM
        GGUBFS=AM*IOR(IAND(IM,IEOR(IX,IY)),1)
        END FUNCTION GGUBFS
