c     PROGRAM OF DECONVOLUTION
c     ------------------------

c     VARIABLE DECLARATION
c     --------------------

      implicit none

      include 'param_dec.f'

      real*4 g(nbr1,nbg),sig2(nbr1,nbg),sig3(nbr1,nbg)
      real*4 s(nbr2,nbg),ss(2*xy2,nbg),trs(nbr2,nbg),strs(2*xy2,nbg)
      real*4 p(nbr2),lambda(nbr1)
      real*4 a(npic,nbg),b,c(npic,2),z(nbg,2),delta(nbg,2)
      real*4 nda(npic),ndc,ndp,ndz1,ndz2,ndd,k1
      real*4 ndamin(npic),ndcmin,ndpmin,ndz1min,ndz2min,nddmin,k3
      real*4 da(npic,nbg),dc(npic,2),dp(nbr2),dz(nbg,2),dd(nbg,2)
      real*4 rr(nbr2),srr(2*xy2),frf(nbr2)
      real*4 tmp(nbr2),tmp2(nbr2),tmp3(nbr1),tmp4(nbr1)
      real*4 stmp2(2*xy2),stmp(2*xy2)
      real*4 idlc,idfc,som,fwhm,expo,t1,t2,bg,qimin

      integer n1,n2,m,nxy1,nxy2,np,nax1,nax2,pas
      integer i,ii,i2,j,j2,k,k2,l
      integer iter,itmax,itmax2,idl,idf,module,cent,half
      integer ind,ind1,ind2,indice,posi
      
      character*20 name
      character*20 lambda_name,fond_name,idln,idfn
      character*20 entree,sortie,sortie2
      character*70 entete
      character*1 lettre,question,lettre2,lettre3
      character*8 fin
      character*42 ligne
        
      common /tailles/ n1,n2,np,m
      common /param/ b,nxy1,nxy2,pas
      common /gr1/ g,s,sig2,ss
      common /sig/ sig3
      common /trans/ trs,strs
      common /variables/ p,z,a,c,delta
      common /delta2/ da,dc,dp,dz,dd
      common /deltacons/ nda,ndc,ndp,ndz1,ndz2,ndd
      common /deltacons2/ ndamin,ndcmin,ndpmin,ndz1min,ndz2min,nddmin
      common /entropie/ rr,frf,srr
      common /lambd/ idl,idlc,idln
      common /positif/ posi
      common /gaussienne/ fwhm,bg

c     START OF THE PROGRAM
c     --------------------

c     READING OF FILE:  DECONV2.DAT:
c     - - - - - - - - - - - - - - -

      open(unit=4,file='deconv.txt')

      read(4,'(A70)') entete
      read(4,'(A42,I4)') ligne,m
      write(*,*) m
      read(4,'(A70)') entete
      read(4,'(A42,I4)') ligne,np
      read(4,'(A70)') entete
      read(4,'(A42,F16.10)') ligne,fwhm
      if (fwhm .ge. 2.) then
         b = 4.*log(2.)/(fwhm*fwhm)
      else
         write(*,*) 'FWHM < 2 --> erreur echantionnage'
         stop
      end if
      read(4,'(A70)') entete
      read(4,'(A42,F16.10)') ligne,fwhm
      if (fwhm .ge. 2.) then
         bg = 4.*log(2.)/(fwhm*fwhm)
      else
         write(*,*) 'FWHM < 2 --> erreur echantionnage'
         stop
      end if
      read(4,'(A70)') entete
      read(4,'(A42,A20)') ligne,entree
      read(4,'(A70)') entete
      read(4,'(A42,A20)') ligne,sortie
      read(4,'(A70)') entete
      read(4,'(A42,A20)') ligne,sortie2
      read(4,'(A70)') entete
      read(4,'(A42,I6)') ligne,itmax
      read(4,'(A70)') entete
      read(4,'(A42,F16.10)') ligne,qimin
      read(4,'(A70)') entete
      read(4,'(A42,I6)') ligne,itmax2
      read(4,'(A70)') entete
      read(4,'(A42,I4)') ligne,module
      read(4,'(A70)') entete
      read(4,'(A42,A1)') ligne,question
      if (question .eq. 'y') then
         posi = 1
      else
         posi = 0
      end if
      read(4,'(A70)') entete
      read(4,'(A42,A1)') ligne,question
      if (question.eq.'y') then
         idl = 1
         read(4,'(A42,A20)') ligne,idln
      else
         idl = 2
         read(4,'(A42,F16.10)') ligne,idlc
      endif
      read(4,'(A70)') entete
      read(4,'(A42,A1)') ligne,question
      if (question.eq.'y') then
         idf = 1
         read(4,'(A42,A20)') ligne,idfn
      else
         idf = 2
         read(4,'(A42,F16.10)') ligne,idfc
      endif
      close(unit=4)

c     READING OF DATA FROM INPUT FILE 
c     - - - - - - - - - - - - - - - -
      write(*,*) np,m
      open(unit=4,file=entree)
      do i=1,3
         read(4,'(A70)') entete
      end do

      read(4,'(A42,F16.10)') ligne,k1
      read(4,'(A42,F16.10)') ligne,k3
      read(4,'(A70)') entete
      read(4,'(A42,F16.10)') ligne,ndc
      read(4,'(A42,F16.10)') ligne,ndcmin
      read(4,'(A70)') entete
      read(4,'(A42,F16.10)') ligne,ndp
      read(4,'(A42,F16.10)') ligne,ndpmin
      read(4,'(A70)') entete
      read(4,'(A42,F16.10)') ligne,ndz1
      read(4,'(A42,F16.10)') ligne,ndz1min
      read(4,'(A42,F16.10)') ligne,ndz2
      read(4,'(A42,F16.10)') ligne,ndz2min
      read(4,'(A70)') entete
      read(4,'(A42,F16.10)') ligne,ndd
      read(4,'(A42,F16.10)') ligne,nddmin
      do i=1,5
         read(4,'(A70)') entete
      end do
      do k=1,np
         read(4,*) (a(k,j),j=1,m)
         read(4,*) c(k,1),c(k,2)
         write(*,*) k,c(k,1),c(k,2)
      end do
      do i=1,4
         read(4,'(A70)') entete
      end do
      do i=1,m
         read(4,*) z(i,1),z(i,2),delta(i,1),delta(i,2)
      end do
      close(unit=4)

c     LOADING OF IMAGES
c     - - - - - - - - -
      do ii=1,m
         lettre = char(48+int(ii/100))
         i2 = ii - int(ii/100)*100
         lettre2 = char(48+int(i2/10))
         lettre3 = char(48+mod(i2,10))
         fin = lettre // lettre2 // lettre3 // '.fits'

         name = 'g' // fin
         call openima(name,nax1,nax2,tmp3)
         expo = log(real(nax1))/log(2.)
         if (nax1 .ne. nax2) then
            write(*,*) 'Erreur: image g.fits non carrée'
            stop
         else if (int(expo) .ne. expo) then
            write(*,*) 'Erreur: g.fits pas une puissance de 2'
            stop
         end if
         n1 = nax1*nax1
         nxy1 = nax1

         name = 's' // fin
         call openima(name,nax1,nax2,tmp)
         expo = log(real(nax1))/log(2.)
         if (nax1 .ne. nax2) then
            write(*,*) 'Erreur: image s.fits non carrée'
            stop
         else if (nax1 .ne. xy2) then
            write(*,*) 'Erreur: XY2 not the same in images and param.f'
            write(*,*) 'use ini_dec.f'
         else if (int(expo) .ne. expo) then
            write(*,*) 'Erreur: s.fits pas une puissance de 2'
            stop
         end if
         n2 = nax1*nax1
         nxy2 = nax1
         pas = int(xy2/xy1)

         name = 'sig' // fin
         call openima(name,nax1,nax2,tmp4)
         expo = log(real(nax1))/log(2.)
         if (nax1 .ne. nax2) then
            write(*,*) 'Erreur: image sig.fits non carrée'
            stop
         else if (nax1 .ne. xy1) then
            write(*,*) 'image g.fits et sig.fits pas de meme taille or'
            write(*,*) 'Erreur: XY1 not the same in images and param.f'
            write(*,*) 'use ini_dec.f'
            stop
         else if (int(expo) .ne. expo) then
            write(*,*) 'Erreur: sig.fits pas une puissance de 2'
            stop
         end if

         do i=1,n1
            g(i,ii) = tmp3(i)
            sig2(i,ii) = tmp4(i)*tmp4(i)
         end do

         do j=1,nxy2
            do i=1,nxy2
               ind1 = i + (j-1)*nxy2
               ind2 = j + (i-1)*nxy2
               tmp2(ind2) = tmp(ind1)
            end do
         end do

         call rlft2(tmp,stmp,nxy2,1)
         call rlft2(tmp2,stmp2,nxy2,1)
         do i=1,2*nxy2
            ss(i,ii) = stmp(i)
            strs(i,ii) = stmp2(i)
         end do
         do i=1,n2
            s(i,ii) = tmp(i)
            trs(i,ii) = tmp2(i)
         end do
      end do

c     INITIALISATION OF DATA
c     - - - - - - - - - - - -

      do k=1,np
         nda(k) = a(k,1)*k1/(100.)
         ndamin(k) = a(k,1)*k3/(100.)
         do j=1,m
            da(k,j) = a(k,j)*k1/(100.)
            write(*,'(E12.6,A1,$)') a(k,j),' '
         end do
         dc(k,1) = ndc
         dc(k,2) = ndc
         write(*,*) c(k,1),c(k,2)       
      end do
      do j=1,m
         if (j .eq. 1) then
            dz(j,1) = 0.
            dz(j,2) = 0.
            dd(j,1) = 0.
            dd(j,2) = 0.
         else
            dz(j,1) = ndz1
            dz(j,2) = ndz2
            dd(j,1) = ndd
            dd(j,2) = ndd
         end if
         write(*,*) z(j,1),z(j,2)
         write(*,*) 'delta: ',delta(j,1),delta(j,2)
      end do
      write(*,*)

c     INITIALISATION OF TABLE LAMBDA
c     - - - - - - - - - - - - - - - -

      if (idl .eq. 1) then
         lambda_name = idln
         call openima(lambda_name,nax1,nax2,tmp)
         if (nax1 .eq. nxy1) then
            do ii=1,m 
               do j=1,n1
                  lambda(j) = tmp(j)
               end do
            end do
         else
            do j=1,nxy1
               do ii=1,nxy1
                  ind1 = nxy2*pas*(j-1)+pas*(ii-1)
                  indice = (j-1)*nxy1+ii
                  lambda(indice) = 0.
                  do l=1,pas
                     ind2 = (l-1)*nxy2
                     do k2=1,pas
                        lambda(indice) = lambda(indice) + 
     &                                   tmp(ind1+k2+ind2)  
                     end do
                  end do
                  lambda(indice) = lambda(indice)/real(pas*pas)
               end do
            end do
         end if
      else
         do j=1,n1
            lambda(j) = idlc
         end do
      end if

c     CALCULATION OF 1/(SIGMA AU CARRE) AND LAMBDA/(SIGMA AU CARRE):
c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      do ii=1,m
         do j=1,n1
            sig3(j,ii) = 1./sig2(j,ii)
            sig2(j,ii) = lambda(j)/sig2(j,ii)
         end do
      end do
         
c     READING OF INITIAL BACKGROUND IMAGE:
c     - - - - - - - - - - - - - - - - - -
      if (idf .eq. 1) then
         fond_name = idfn
         call openima(fond_name,nax1,nax2,p)
         if (posi .eq. 1) then
            do j=1,n2
               if (p(j) .lt. 0.) p(j) = 0.
            end do
         end if
      else
         do j=1,n2
            p(j) = idfc
         end do
      end if
      do j=1,n2
         dp(j) = ndp
      end do


c     CREATION OF TABLE R:
c     - - - - - - - - - - 
      som = 0.
      cent = (nxy2/2+1)
      do j=1,nxy2
         t2 = j-cent
         ind = (j-1)*nxy2
         do i=1,nxy2
            t1 = i-cent
            ind = ind+1
            tmp(ind) = exp(-b*(t1*t1+t2*t2))
            som = som + tmp(ind)
         end do
      end do

      half = int(nxy2/2)

      do j=1,nxy2
         j2 = j+half
         if (j2 .gt. nxy2) j2 = j2-nxy2
         ind = (j-1)*nxy2
         ind2 = (j2-1)*nxy2
         do i=1,nxy2
            i2 = i+half
            if (i2 .gt. nxy2) i2 = i2-nxy2
            rr(ind+i) = tmp(ind2+i2)/som
         end do
      end do
           
      call rlft2(rr,srr,nxy2,1)

c     START OF THE MINIMISATION :
c     - - - - - - - - - - - - - -
      open(unit=1,file=entree)
      open(unit=3,file=sortie2)
      do i=1,24
         read(1,'(A70)') entete
         write(3,'(A70)') entete
      end do
      close(unit=1)

      open(unit=2,file=sortie)
      write(2,'(A10)') 'Results : '
      write(2,'(A10)') '========= '
      write(2,*)

      call minimi(iter,itmax,module,qimin)
      write(2,'(A34,I4)') 'Number of iterations in minimi  : ',iter 
      call decompos()
      call residu()
      write(*,*)
      write(*,*) 'END OF DECONVOLUTION'

      close(unit=2)
      close(unit=3)

      end 

c     --------------------
c     SUBROUTINE: minimi
C     --------------------
c     MINIMISATION BY DECREASING STEPS.
c     --------------------

      subroutine minimi(iter,itmax,modu,qimin)

      implicit none

      include 'param_dec.f'

      real*4 eps
      parameter (eps=1.E-10)

      real*4 z(nbg,2),p(nbr2),delta(nbg,2),a(npic,nbg),c(npic,2)
      real*4 zN(nbg,2),pN(nbr2),deltaN(nbg,2),aN(npic,nbg),cN(npic,2)
      real*4 da(npic,nbg),dc(npic,2),dp(nbr2),dz(nbg,2),dd(nbg,2)
      real*4 daN(npic,nbg),dcN(npic,2),dpN(nbr2),dzN(nbg,2),ddN(nbg,2)
      real*4 nda(npic),ndc,ndp,ndz1,ndz2,ndd
      real*4 ndamin(npic),ndcmin,ndpmin,ndz1min,ndz2min,nddmin
      real*4 fp,qi,fpA,vecteurqi(10)
      real*4 xia(npic,nbg),xic(npic,2),xid(nbg,2),xiz(nbg,2),xip(nbr2)
      real*4 xiaA(npic,nbg),xicA(npic,2),xidA(nbg,2),xizA(nbg,2)
      real*4 xipA(nbr2),qimin
      integer posi,indi,condi,its,ind2,itmax,iter,modu
      integer n1,n2,np,m,ii,jj
      integer i,j,arret,arret2,k1,k10
      

      common /delta2/ da,dc,dp,dz,dd
      common /deltacons/ nda,ndc,ndp,ndz1,ndz2,ndd
      common /deltacons2/ ndamin,ndcmin,ndpmin,ndz1min,ndz2min,nddmin
      common /variables/ p,z,a,c,delta
      common /tailles/ n1,n2,np,m
      common /derivees/ xia,xic,xid,xiz,xip
      common /minimum/ fp,qi
      common /positif/ posi

      arret = 0
      arret2 = 1

c     COMPUTATION OF MINIMUM AND DERIVTIVES:
c     - - - - - - - - - - - - - - - - - - -
      call inifunc_dfunc(p,z,a,c,delta)
      vecteurqi(1) = qi
      do i=2,10
         vecteurqi(i) = 0.
      end do

c     CHECK OF THE SIGN OF THE DERIVATIVES OF F 
c     AND MODIFICATION OF STEP.
c     - - - - - - - - - - - - - - - - - - - - - 
      do i=1,n2
        dpN(i) = dp(i)
      end do
      do i=1,np
         do j=1,m
            daN(i,j) = da(i,j)
         end do
         dcN(i,1) = dc(i,1)
         dcN(i,2) = dc(i,2)
      end do
      do j=1,m
         do i=1,2
            dzN(j,i) = dz(j,i)
            ddN(j,i) = dd(j,i)
         end do
      end do

      indi = m+2
      condi = indi*np

      do its=1,itmax

         write(*,*) its,' M   min qi: ',fp,qi
         write(*,*)
         iter = its

c        MODIFICATION OF F
c        - - - - - - - - -
         do i=1,n2
            if (xip(i) .lt. 0.) then
                pN(i) = p(i) + dpN(i)
            else
               if (posi .eq. 0) then
                  pN(i) = p(i) - dpN(i)
               else
                  if (p(i)-dpN(i) .ge. 0) then
                      pN(i) = p(i) - dpN(i)
                  else
                      pN(i) = 0.
                  end if
              end if
            end if
         end do

c        MODIFICATION OF A, C1 AND C2
c        - - - - - - - - - - - - - - -
         do i=1,np
            write(*,'(''a,c :'',$)')
            do j=1,m
               if (xia(i,j) .lt. 0.) then
                  aN(i,j) = a(i,j) + daN(i,j)
               else
                  aN(i,j) = a(i,j) - daN(i,j)
               end if
               write(*,'(E12.6,A1,$)') aN(i,j),' '
            end do
            if(xic(i,1) .lt. 0.) then
               cN(i,1) = c(i,1) + dcN(i,1)
            else
               cN(i,1) = c(i,1) - dcN(i,1)
            end if
            if(xic(i,2) .lt. 0.) then
               cN(i,2) = c(i,2) + dcN(i,2)
            else
               cN(i,2) = c(i,2) - dcN(i,2)
            end if
            write(*,'(2(E12.6,A1))') cN(i,1),' ',cN(i,2),' '
         end do

c        MODIFICATION OF Z AND DELTA 
c           NOTE: THE PARAMETER OF IMAGE 1 ARE NOT MODIFIED
c        - - - - - - - - - - - - - - - - - - - - - - - - - -
         write(*,*) 'z:     ',z(1,1),z(1,2)
         zN(1,1) = z(1,1)
         zN(1,2) = z(1,2)
         write(*,*) 'delta: ',delta(1,1),delta(1,2)
         deltaN(1,1) = delta(1,1)
         deltaN(1,2) = delta(1,2)
         do i=2,m
            if (xiz(i,1) .lt. 0.) then
               zN(i,1) = z(i,1) + dzN(i,1)
            else
               zN(i,1) = z(i,1) - dzN(i,1)
            end if
            if (xiz(i,2) .lt. 0.) then
               zN(i,2) = z(i,2) + dzN(i,2)
            else
               zN(i,2) = z(i,2) - dzN(i,2)
            end if
            if (xid(i,1) .lt. 0.) then
               deltaN(i,1) = delta(i,1) + ddN(i,1)
            else
               deltaN(i,1) = delta(i,1) - ddN(i,1)
            end if
            if (xid(i,2) .lt. 0.) then
               deltaN(i,2) = delta(i,2) + ddN(i,2)
            else
               deltaN(i,2) = delta(i,2) - ddN(i,2)
            end if
            write(*,*) 'z:     ',zN(i,1),zN(i,2)
            write(*,*) 'delta: ',deltaN(i,1),deltaN(i,2)
         end do

         fpA = fp
         do i=1,n2
            xipA(i) = xip(i)
         end do
         do j=1,np
            do i=1,m
               xiaA(j,i) = xia(j,i)
            end do
            xicA(j,1) = xic(j,1)
            xicA(j,2) = xic(j,2)
         end do
         do j=1,m
            xizA(j,1) = xiz(j,1)
            xizA(j,2) = xiz(j,2)
            xidA(j,1) = xid(j,1)
            xidA(j,2) = xid(j,2)
         end do

         write(*,*) '---------------------------------------'

c        AUTOMATIC BACKUP AFTER X ITERATIONS:
c        - - - - - - - - - - - - - - - - - - 
         call inifunc_dfunc(pN,zN,aN,cN,deltaN)
         ii = mod(its+1,10)
         jj = ii+1
         if (ii.eq.0) ii = 10
         vecteurqi(ii) = qi
         if (abs(abs(vecteurqi(ii)) - 
     &        abs(vecteurqi(jj))) .le. qimin) then
            write(*,*)
            write(*,*) '===================================='
            write(*,*) '  LA DIFFERENCE ENTRE 10 KHI CARRE'
            write(*,*) '  EST INSIGNIFIANTE'
            write(*,*) '===================================='
            return
         end if
  
         if (mod(its,modu) .eq. 0) then
            write(2,*) 'Number of iterations in minimi  : ',iter 
            write(2,*) 
            call decompos()
            do j=1,9+(2+m)*np+m
               backspace(unit=2)
            end do
            do j=1,2*np+4+m
               backspace(unit=3)
            end do
         end if

c        TEST OF MINIMUM, IF THE MINIMUM INCREASES THEN STOP
c        - - - - - - - - - - - - - - - - - - - - - - - - - -
         if (fpA .gt. fp) then
            
            if (arret .ne. 0) arret = 0

c           FOND
c           - - - 
            do i=1,n2
               dp(i) = dpN(i)
               if (xipA(i)*xip(i) .lt. 0.)  then
                  dpN(i) = dp(i)/2.
               else
                  if (dp(i)*1.1 .le. ndp) then
                     dpN(i)= dp(i)*1.1
                  else
                     dpN(i) = dp(i)
                  end if
               end if
               if (arret2.eq.1 .and. dpN(i) .gt. ndpmin) arret2 = 0 
               p(i) = pN(i)
            end do

c           SOURCES
c           - - - - 
            do i=1,np
               do j=1,m
                 da(i,j) = daN(i,j)
                 if (xiaA(i,j)*xia(i,j) .lt. 0.) then
                     daN(i,j) = da(i,j)/2.
                  else
                     if (da(i,j)*1.1 .le. nda(i)) then
                        daN(i,j) = da(i,j)*1.1
                     else
                        daN(i,j) = da(i,j)
                     end if
                  end if
                  a(i,j) = aN(i,j)
                  if (arret2.eq.1 .and. daN(i,j).gt.ndamin(i)) arret2=0
               end do
               dc(i,1) = dcN(i,1)
               dc(i,2) = dcN(i,2)
               if (xicA(i,1)*xic(i,1) .lt. 0.) then
                  dcN(i,1) = dc(i,1)/2.
               else
                  if (dc(i,1)*1.1 .le. ndc) then
                     dcN(i,1) = dc(i,1)*1.1
                  else
                     dcN(i,1) = dc(i,1)
                  end if
               end if
               if (xicA(i,2)*xic(i,2) .lt. 0.) then
                  dcN(i,2) = dc(i,2)/2.
               else
                  if (dc(i,2)*1.1 .le. ndc) then
                     dcN(i,2) = dc(i,2)*1.1
                  else
                     dcN(i,2) = dc(i,2)
                  end if
               end if
               if (arret2.eq.1 .and. dcN(i,1) .gt. ndcmin) arret2 = 0
               if (arret2.eq.1 .and. dcN(i,2) .gt. ndcmin) arret2 = 0
               c(i,1) = cN(i,1)
               c(i,2) = cN(i,2)
            end do

c           CHECK OF THE SIGN OF THE DERIVATIVES OF Z AND DELTA 
c           AND MODIFICATION OF STEP.
c           - - - - - - - - - - - - - - - - - - - - - - - - - - 
            z(1,1) = zN(1,1)
            z(1,2) = zN(1,2)
            delta(1,1) = deltaN(1,1)
            delta(1,2) = deltaN(1,2)
            dz(1,1) = dzN(1,1)
            dz(1,2) = dzN(1,2)
            dd(1,1) = ddN(1,1)
            dd(1,2) = ddN(1,2)
            ind2 = n2+indi*np
            do j=2,m
               dz(j,1) = dzN(j,1)
               dz(j,2) = dzN(j,2)
               if (xizA(j,1)*xiz(j,1) .lt. 0.) then
                  dzN(j,1) = dz(j,1)/2.
               else 
                  if (dz(j,1)*1.1 .le. ndz1)  then
                     dzN(j,1) = dz(j,1)*1.1
                  else
                     dzN(j,1) = dz(j,1)
                  end if
               end if
               if (xizA(j,2)*xiz(j,2) .lt. 0.) then
                  dzN(j,2) = dz(j,2)/2.
               else
                  if (dz(j,2)*1.1 .le. ndz2)  then
                     dzN(j,2) = dz(j,2)*1.1
                  else
                     dzN(j,2) = dz(j,2)
                  end if
               end if
               if (arret2.eq.1 .and. dzN(j,1) .gt. ndz1min) arret2 = 0
               if (arret2.eq.1 .and. dzN(j,2) .gt. ndz2min) arret2 = 0
               z(j,1) = zN(j,1)
               z(j,2) = zN(j,2)
               dd(j,1) = ddN(j,1)
               dd(j,2) = ddN(j,2)
               if (xidA(j,1)*xid(j,1) .lt. 0.) then
                  ddN(j,1) = dd(j,1)/2.
               else
                 if (dd(j,1)*1.4 .le. ndd)  then
                    ddN(j,1) = dd(j,1)*1.4
                 else
                    ddN(j,1) = dd(j,1)
                 end if
               end if
               if (xidA(j,2)*xid(j,2) .lt. 0.) then
                  ddN(j,2) = dd(j,2)/2.
               else
                 if (dd(j,2)*1.4 .le. ndd)   then
                    ddN(j,2) = dd(j,2)*1.4
                 else
                    ddN(j,2) = dd(j,2)
                 end if
               end if
               if (arret2.eq.1 .and. ddN(j,1) .gt. nddmin) arret2 = 0
               if (arret2.eq.1 .and. ddN(j,2) .gt. nddmin) arret2 = 0
               delta(j,1) = deltaN(j,1)
               delta(j,2) = deltaN(j,2)
           end do

           if (arret2 .eq. 1) then
               write(*,*)
               write(*,*) '================================='
               write(*,*) '  TOUT LES PAS SONT TROP PETITS'
               write(*,*) '================================='
               return
           end if
           arret2 = 1

         else

            arret = arret + 1
            write(*,*)
            write(*,*) '========================'
            write(*,*) '  LE MINIMUM AUGMENTE'
            write(*,*) '========================'
            write(*,*) 'AUGMENTATION N°',arret
            write(*,*)

            if (arret .eq. 10) then
               write(*,*)
               write(*,*) '================================='
               write(*,*) '  LE MINIMUM A AUGMENTE 10 fois'
               write(*,*) '================================='
               return
            end if

c           FOND
c           - - - 
            do i=1,n2
               if (xipA(i)*xip(i) .lt. 0.)  then
                  dp(i) = dp(i)/2.
               end if
            end do

c           SOURCES
c           - - - - 
            do i=1,np
               do j=1,m
                  if (xiaA(i,j)*xia(i,j) .lt. 0.) then
                     da(i,j) = da(i,j)/2.
                  end if
               end do
               if (xicA(i,1)*xic(i,1) .lt. 0.) then
                  dc(i,1) = dc(i,1)/2.
               end if
               if (xicA(i,2)*xic(i,2) .lt. 0.) then
                  dc(i,2) = dc(i,2)/2.
               end if
            end do

c           Z ET DECALLAGE
c           - - - - - - - -
            do j=2,m
               if (xizA(j,1)*xiz(j,1) .lt. 0.) then
                  dz(j,1) = dz(j,1)/2.
               end if
               if (xizA(j,2)*xiz(j,2) .lt. 0.) then
                  dz(j,2) = dz(j,2)/2.
               end if
               if (xidA(j,1)*xid(j,1) .lt. 0.) then
                  dd(j,1) = dd(j,1)/2.
               end if
               if (xidA(j,2)*xid(j,2) .lt. 0.) then
                  dd(j,2) = dd(j,2)/2.
               end if
           end do

         end if
      end do

      write(*,*)
      write(*,*) '============================='
      write(*,*) '  NOMBRE ITERATIONS EXEDEES'
      write(*,*) '============================='

      return
      end


c     -------------------------
c     SUBROUTINE: inifunc_dfunc
c     -------------------------
c     ROUTINE WITH COMPUTE THE MINIMUM OF THE FUNCTION AND SOME TABLES 
c     WHICH WILL BE USED IN DFUNC
c     ---------------------

      subroutine inifunc_dfunc(p,z,a,c,delta)

      implicit none

      include 'param_dec.f'

      integer nb4,nb
      parameter(nb4=nbr2+4*(xy2+1))
      parameter (nb=nbr2+(nbg+2)*npic+(nbg-1)*4)

      real*4 p(nbr2),z(nbg,2),delta(nbg,2),a(npic,nbg),c(npic,2),b
      real*4 g(nbr1,nbg),sig2(nbr1,nbg),s(nbr2,nbg),ss(2*xy2,nbg)
      real*4 sig3(nbr1,nbg)
      real*4 far(nbr2,nbg),far2(nbr2,nbg)
      real*4 sf(nbr2),sf2(nbr2),ssf(2*xy2),ssf2(2*xy2)
      real*4 sfbis(nbr1,nbg),sfbis2
      real*4 tmp(nbr2),stmp(2*xy2),tmp2(nbr2),stmp2(2*xy2)
      real*4 tmp3(nbr2),tmp4(nbr2),tmp5(nbr2)
      real*4 stmp3(2*xy2),stmp4(2*xy2),stmp5(2*xy2)
      real*4 rr(nbr2),srr(2*xy2),frf(nbr2)
      real*4 w(9),inter(nb4),shfz(nb4,nbg)
      real*4 tamp1,tamp2,tm,delx,dely
      real*4 fp,qi,sum2,sum1,sum5,tamp,t1,t2,fwhm,bg
      real*4 som1,som2,cx,cy,ebc,tamp3,srcbis,srcbis2,srcbis3
      real*4 xia(npic,nbg),xic(npic,2),xid(nbg,2),xiz(nbg,2),xip(nbr2)
      real*4 trs(nbr2,nbg),strs(2*xy2,nbg)
      integer n1,n2,np,m,nxy1,nxy2,pas,debutx,debuty,finx,finy
      integer ind,ind1,ind2,ind3,indice,indice2
      integer i,ii,j,k,l,deb,indi,ind4,k2,dif

      common /gr1/ g,s,sig2,ss
      common /param/ b,nxy1,nxy2,pas
      common /gr5/ shfz,sfbis
      common /entropie/ rr,frf,srr
      common /sig/ sig3
      common /poids/ w
      common /tailles/ n1,n2,np,m
      common /minimum/ fp,qi
      common /trans/ trs,strs
      common /derivees/ xia,xic,xid,xiz,xip
      common /gaussienne/ fwhm,bg

c     CALCULATION OF : FAR(i,j) = z(i,1)*F(j) + z(i,2)
c                      FAR2(i,j) = SOMMEsurK (A(k,i)*R(C(k)-delta(i)))
c                      WHERE i = IMAGE NUMBER
c                            j = INDEX OF THE COMPUTED PARAMETER
c                            k = NUMBER OF POINT SOURCES
c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

      indi = m+2

      do i=1,n2
         tmp2(i) = p(i)
         xip(i) = 0.
      end do
      do ii=1,m
         do i=1,n2
            far(i,ii) = z(ii,1)*p(i)+z(ii,2)
            far2(i,ii) = 0.
         end do
      end do

c     CALCULATION OF R CONVOLVED WITH F, USED IN THE SMOOTHING TERM
c     AND CALCULATION OF THIS TERM:
c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      call rlft2(tmp2,stmp2,nxy2,1)
      call produit(tmp2,rr,stmp2,srr,xy2)
      call rlft2(tmp2,stmp2,nxy2,-1)

      sum2 = 0.
      do i=1,n2
         frf(i) = p(i)-tmp2(i)
         sum2 = sum2 + frf(i)**2
      end do

      do k=1,np
         xic(k,1) = 0.
         xic(k,2) = 0.
         do ii=1,m
            t1 = c(k,1)-delta(ii,1)
            t2 = c(k,2)-delta(ii,2)
            xia(k,ii) = 0.
            do j=1,nxy2
               tamp2 = j-t2
               ind2 = (j-1)*nxy2
               do i=1,nxy2
                  ind = ind2+i
                  tamp1 = i-t1
                  tm = exp(-b*(tamp1*tamp1+tamp2*tamp2))
                  far2(ind,ii) = far2(ind,ii) + a(k,ii)*tm
               end do
            end do
         end do
      end do


c     CALCULATION OF S CONVOLVED WITH FAR AND S CONVOLVED WITH FAR2,
c     AND CALCULATION OF THE MINIMUM
c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      sum1 = 0.
      qi = 0.
      sum5 = 0.
      do ii=1,m
         xid(ii,1) = 0.
         xid(ii,2) = 0.
         xiz(ii,1) = 0.
         xiz(ii,2) = 0.

c        CALCULATION OF W (WEIGHT OF PIXELS) FOR THE SHIFT BETWEEN IMAGES
c        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         if (delta(ii,1) .lt. 0.) then
            delx = 1+delta(ii,1)
         else
            delx = delta(ii,1)
         end if
         if (delta(ii,2) .lt. 0.) then
            dely = 1+delta(ii,2)
         else
            dely = delta(ii,2)
         end if
         if (pas .eq. 1) then
c           WHEN ONE STAYS IN BIG PIXELS
c           - - - - - - - - - - - - - - 
            w(1) = (1-delx)*(1-dely)
            w(2) = delx*(1-dely)
            w(3) = dely*(1-delx)
            w(4) = delx*dely
         else
c           WHEN ONE CHANGES TO SMALL PIXELS
c           - - - - - - - - - - - - - - -
            w(1) = (1-delx)*(1-dely)/4
            w(2) = (1-dely)/4
            w(3) = delx*(1-dely)/4
            w(4) = (1-delx)/4
            w(5) = 0.25
            w(6) = delx/4
            w(7) = dely*(1-delx)/4
            w(8) = dely/4
            w(9) = delx*dely/4
         end if
         if (delta(ii,1) .ge. 0. .and. delta(ii,2) .ge. 0.) then
            deb = nxy2+3
         else if (delta(ii,1) .ge. 0. .and. delta(ii,2) .lt. 0.) then
            deb = 1
         else if (delta(ii,1) .lt. 0. .and. delta(ii,2) .ge. 0.) then
            deb = nxy2+2
         else
            deb = 0
         end if

         do i=1,n2
            sf(i) = far(i,ii)
            sf2(i) = far2(i,ii)
            tmp(i) = s(i,ii)
            tmp2(i) = 0.
            tmp3(i) = 0.
         end do   
         do i=1,xy2*2
            stmp(i) = ss(i,ii)
         end do

         call rlft2(sf,ssf,nxy2,1)
         call rlft2(sf2,ssf2,nxy2,1)
         call produit(sf,tmp,ssf,stmp,nxy2)
         call produit(sf2,tmp,ssf2,stmp,nxy2)
         call rlft2(sf,ssf,nxy2,-1)
         call rlft2(sf2,ssf2,nxy2,-1)

c        SF = (S)*(Z1*F+Z2) AND SF2 = (S)*(SOMME(A*R(C-DELTE))

c        SF MUST BE SHIFTED
c        THE TABLE IS MODIFIED AND THE STARTING POINT IS COMPUTED
c        YOU CHANGE FROM BIG PIXELS TO SMALL PIXELS TO TAKE INTO
c        ACCOUNT THE WEIGHTS OF PIXELS.
c        - - - - - - - - - - - - - - - - - - - - - - - - - - -
c        SFBIS2 = SF SHIFTED AND PUT IN BIG PIXELS.
c        SF2 IS ALREADY SHIFTED, ON ONLY CHANGE FROM SMALL TO BIG PIXELS
c        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

         do j=1,nxy2
            inter(j*(nxy2+2)+1) = sf(j*nxy2)
            do i=1,nxy2
               indice = (j-1)*nxy2+i
               indice2 = j*(nxy2+2)+1+i
               inter(indice2) = sf(indice)
            end do
            inter((j+1)*(nxy2+2)) = sf((j-1)*nxy2+1)
         end do
         do i=1,nxy2+2
            inter(i) = inter((nxy2+2)*nxy2+i)
            inter((nxy2+2)*(nxy2+1)+i) = inter(nxy2+2+i)
         end do

         do i=1,n2+4*(nxy2+1)
            shfz(i,ii) = inter(i)
         end do

         if (delta(ii,1) .ge. 0. .and. delta(ii,2) .ge. 0.) then
            deb = nxy2+3
         else if (delta(ii,1) .ge. 0. .and. delta(ii,2) .lt. 0.) then
            deb = 1
         else if (delta(ii,1) .lt. 0. .and. delta(ii,2) .ge. 0.) then
            deb = nxy2+2
         else
            deb = 0
         end if 

         do j=1,nxy1
            do i=1,nxy1
               ind1 = deb+(nxy2+2)*pas*(j-1)+pas*(i-1)
               indice = (j-1)*nxy1+i
               sfbis2 = 0.
               sfbis(indice,ii) = 0.
               do l=1,pas+1
                  ind2 = (l-1)*(nxy2+2)
                  ind3 = (pas+1)*(l-1)
                  do k=1,pas+1
                    sfbis2 = sfbis2 + 
     &                       w(ind3+k)*inter(ind1+k+ind2) 
                  end do
               end do
               ind1 = nxy2*pas*(j-1)+pas*(i-1)
               do l=1,pas
                  ind2 = (l-1)*nxy2
                  do k=1,pas
                    sfbis(indice,ii) = sfbis(indice,ii) + 
     &                              sf2(ind1+k+ind2) 
                  end do
               end do
               sfbis(indice,ii) = sfbis(indice,ii)/real(pas**2) +
     &                            sfbis2
               tamp = (sfbis(indice,ii)-g(indice,ii))**2 
               qi = qi + tamp*sig3(indice,ii)
               sum1 = sum1 + tamp*sig2(indice,ii)
            end do
         end do

c        SFBIS = SFBIS2 + SF2 PUT IN BIG PIXELS
c        - - - - - - - - - - - - - - - - - - - 

         do k=1,np
            cx = c(k,1)-delta(ii,1)
            cy = c(k,2)-delta(ii,2)
            do j=1,nxy2
               tamp2 = j-cy
               ind2 = (j-1)*nxy2
               do i=1,nxy2
                  ind = ind2+i
                  tamp1 = i-cx
                  ebc = exp(-b*((tamp1*tamp1)+(tamp2*tamp2)))
                  tmp2(ind) = tmp2(ind) + ebc*tamp1*a(k,ii)
                  tmp3(ind) = tmp3(ind) + ebc*tamp2*a(k,ii)
               end do
            end do

            cx = cx/2.
            cy = cy/2.
            debutx = anint(cx)-fwhm
            debuty = anint(cy)-fwhm
            finx = anint(cx)+fwhm
            finy = anint(cy)+fwhm
            if (debutx .le. 0) then
               dif = 1-debutx
               debutx = 1
               finx = finx-dif
            end if
            if (debuty .le. 0) then
               dif = 1-debuty
               debuty=1
               finy = finy-dif
            end if
            if (finx .gt. nxy1) then
               dif = finx-nxy1
               debutx = debutx+dif
               finx=nxy1
            end if
            if (finy .gt. nxy1) then
               dif = finy-nxy1
               debuty = debuty+dif
               finy = nxy1
            end if
            do j=debuty,finy
               tamp2 = j-cy
               do i=debutx,finx
                  tamp1 = i-cx
                  indice = (j-1)*nxy1+i
                  ebc = exp(-bg*((tamp1*tamp1)+(tamp2*tamp2)))
                  tamp3 = (sfbis(indice,ii)-g(indice,ii))
     &                     *sig2(indice,ii)
                  xia(k,ii) = xia(k,ii) + ebc*tamp3
                  xic(k,1) = xic(k,1) + a(k,ii)*ebc*tamp3*tamp1
                  xic(k,2) = xic(k,2) + a(k,ii)*ebc*tamp3*tamp2
               end do
            end do
         end do

c        CALCULATION OF DERIVATIVES OF Z AND DELTA
c        - - - - - - - - - - - - - - - - - - - - - 
         if (ii .ne. 1) then
            do i=1,n2
               tmp4(i) = p(i)
            end do
            call rlft2(tmp4,stmp4,nxy2,1)
            call rlft2(tmp2,stmp2,nxy2,1)
            call rlft2(tmp3,stmp3,nxy2,1)
            call produit(tmp4,tmp,stmp4,stmp,nxy2)
            call produit(tmp2,tmp,stmp2,stmp,nxy2)
            call produit(tmp3,tmp,stmp3,stmp,nxy2)
            call rlft2(tmp4,stmp4,nxy2,-1)
            call rlft2(tmp2,stmp2,nxy2,-1)
            call rlft2(tmp3,stmp3,nxy2,-1)

            do j=1,nxy2
               inter(j*(nxy2+2)+1) = tmp4(j*nxy2)
               do i=1,nxy2
                  indice = (j-1)*nxy2+i
                  indice2 = j*(nxy2+2)+1+i
                  inter(indice2) = tmp4(indice)
               end do
               inter((j+1)*(nxy2+2)) = tmp4((j-1)*nxy2+1)
            end do
            do i=1,nxy2+2
               inter(i) = inter((nxy2+2)*nxy2+i)
               inter((nxy2+2)*(nxy2+1)+i) = inter(nxy2+2+i)
            end do
 
            do j=1,nxy1
               do i=1,nxy1
                  ind1 = deb+(nxy2+2)*pas*(j-1)+pas*(i-1)
                  ind4 = nxy2*pas*(j-1)+pas*(i-1)
                  indice = (j-1)*nxy1+i
                  srcbis = 0.
                  srcbis2 = 0.
                  srcbis3 = 0.
                  tamp2 = (sfbis(indice,ii)-g(indice,ii))
     &                    *sig2(indice,ii)
                  do l=1,pas+1
                     ind2 = (l-1)*(nxy2+2)
                     ind3 = (pas+1)*(l-1)
                     do k2=1,pas+1
                        srcbis = srcbis + 
     &                           w(ind3+k2)*inter(ind1+k2+ind2)  
                     end do
                  end do
                  do l=1,pas
                     ind2 = (l-1)*nxy2
                     do k2=1,pas
                        srcbis2 = srcbis2 + tmp2(ind4+k2+ind2)  
                        srcbis3 = srcbis3 + tmp3(ind4+k2+ind2)  
                        tmp5(ind4+k2+ind2) = tamp2
                     end do
                  end do
                  srcbis2 = srcbis2/real(pas**2)
                  srcbis3 = srcbis3/real(pas**2)
                  if (pas .eq. 1) then
                     som1 = shfz(ind1+1,ii)*(dely-1) +
     &                      shfz(ind1+2,ii)*(1-dely) -
     &                      shfz(ind1+(nxy2+2)+1,ii)*dely +
     &                      shfz(ind1+(nxy2+2)+2,ii)*dely
                     som2 = shfz(ind1+1,ii)*(delx-1) -
     &                      shfz(ind1+2,ii)*delx +
     &                      shfz(ind1+(nxy2+2)+1,ii)*(1-delx) +
     &                      shfz(ind1+(nxy2+2)+2,ii)*delx
                  else
                     som1 = (shfz(ind1+1,ii)*(dely-1) +
     &                      shfz(ind1+3,ii)*(1-dely) -     
     &                      shfz(ind1+(nxy2+2)+1,ii) +
     &                      shfz(ind1+(nxy2+2)+3,ii) -
     &                      shfz(ind1+2*(nxy2+2)+1,ii)*dely +
     &                      shfz(ind1+2*(nxy2+2)+3,ii)*dely)*0.25 
                     som2 = (shfz(ind1+1,ii)*(delx-1) -
     &                      shfz(ind1+2,ii) -
     &                      shfz(ind1+3,ii)*delx +
     &                      shfz(ind1+2*(nxy2+2)+1,ii)*(1-delx)+
     &                      shfz(ind1+2*(nxy2+2)+2,ii) +
     &                      shfz(ind1+2*(nxy2+2)+3,ii)*delx)*0.25
                  end if
                  xiz(ii,1) = xiz(ii,1) + 2*tamp2*srcbis
                  xiz(ii,2) = xiz(ii,2) + 2.*tamp2
                  xid(ii,1) = xid(ii,1)-4.*b*tamp2*srcbis2
     &                 + 2.*som1*tamp2
                  xid(ii,2) = xid(ii,2) - 4.*b*tamp2*srcbis3
     &                 + 2.*som2*tamp2
               end do
            end do
         else 
            do j=1,nxy1
               do i=1,nxy1
                  ind1 = nxy2*pas*(j-1)+pas*(i-1)
                  indice = (j-1)*nxy1+i
                  tamp2 = (sfbis(indice,ii)-g(indice,ii))*
     &                    sig2(indice,ii)
                  do l=1,pas
                     ind2 = (l-1)*nxy2
                     do k=1,pas
                        tmp5(ind1+k+ind2) = tamp2
                     end do
                  end do
               end do
            end do
         end if  

c        CALCULTATION OF DERIVATIVES OF F  
c        - - - - - - - - - - - - - - - - -

         call rlft2(tmp5,stmp5,nxy2,1)
         
         do i=1,n2
            tmp3(i) = trs(i,ii)
         end do
         do i=1,2*nxy2
            stmp3(i) = strs(i,ii)
         end do

         call produit(tmp3,tmp5,stmp3,stmp5,nxy2)
         call rlft2(tmp3,stmp3,nxy2,-1)

         tamp1 = z(ii,1)*2.
         do i=1,n2
            xip(i) =  xip(i) + tamp1*tmp3(i) 
         end do
      end do
 
      fp = sum1 + sum2 + sum5

c     FP = THE MINIMUM VALUE
c     - - - - - - - - - - -
c      do i=1,np
c        write(*,'(''xi a,c: '',$)') 
c        do j=1,m 
c           write(*,'(E12.6,A1,$)') xia(i,j),' '
c        end do
c        write(*,'(2(E12.6,A1))') xic(i,1),' ',xic(i,2),' '
c      end do
c      do i=2,m
c       write(*,*) 'xi z1, z2 : ',xiz(i,1),xiz(i,2)
c       write(*,*) 'xi delta1, delta2 : ',xid(i,1),xid(i,2)
       
c      end do

c     CALCULTATION OF THE SMOOTHING TERM FOR THE DERIVATIVES OF F
c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      do i=1,n2
         tmp2(i) = frf(i)
      end do

      call rlft2(tmp2,stmp2,nxy2,1)
      call produit(tmp2,rr,stmp2,srr,nxy2)
      call rlft2(tmp2,stmp2,nxy2,-1)

      do i=1,n2
         xip(i) = xip(i) + 2.*(frf(i)-tmp2(i))
      end do

c      write(*,*) 'xi f   : ',xip(1),xip(n2/2),xip(n2)
c      write(*,*)

      return
      end


c     ----------------------
c     SUBROUTINE: decompos
c     ----------------------
c     THIS ROUTINE IS USED FOR THE OUTPUT OF THE RESULTS.
c     ----------------------

      subroutine decompos()

      implicit none

      include 'param_dec.f'

      real*4 z(nbg,2),p(nbr2),a(npic,nbg),c(npic,2),b,delta(nbg,2)
      real*4 g(nbr1,nbg),sig2(nbr1,nbg),s(nbr2,nbg),ss(2*xy2,nbg)
      real*4 sig3(nbr1,nbg)
      real*4 sfbis(nbr1,nbg)
      real*4 shfz(nbr2+4*(xy2+1),nbg),inter(nbr2+4*(xy2+1)),w2(4)
      real*4 tmp2(nbr2),stmp2(2*xy2),tmp(nbr2),stmp(2*xy2)
      real*4 tamp1,tamp2,t1,t2,tamp,idlc,fp,qi,delx,dely
      integer n1,n2,np,m,pas,nxy1,nxy2
      integer i,ii,i2,j,k,l,ind1,ind2,indice,indi,idl,deb,indice2
      character*20 results,idln
      character*1 lettre1,lettre3,lettre2
      character*8 fin

      common /param/ b,nxy1,nxy2,pas
      common /gr5/ shfz,sfbis
      common /gr1/ g,s,sig2,ss
      common /lambd/ idl,idlc,idln
      common/sig/ sig3
      common /variables/ p,z,a,c,delta
      common /tailles/ n1,n2,np,m
      common /minimum/ fp,qi

      indi = m+2

      do ii=1,m
         do j=1,nxy2
            inter(j*(nxy2+2)+1) = p(j*nxy2)
            do i=1,nxy2
               indice = (j-1)*nxy2+i
               indice2 = j*(nxy2+2)+1+i
               inter(indice2) = p(indice)
            end do
            inter((j+1)*(nxy2+2)) = p((j-1)*nxy2+1)
         end do
         do i=1,nxy2+2
            inter(i) = inter((nxy2+2)*nxy2+i)
            inter((nxy2+2)*(nxy2+1)+i) = inter(nxy2+2+i)
         end do

         if (delta(ii,1) .lt. 0.) then
            delx = 1+delta(ii,1)
         else
            delx = delta(ii,1)
         end if
         if (delta(ii,2) .lt. 0.) then
            dely = 1+delta(ii,2)
         else
            dely = delta(ii,2)
         end if
         if (delta(ii,1) .ge. 0. .and. delta(ii,2) .ge. 0.) then
            deb = nxy2+3
         else if (delta(ii,1) .ge. 0. .and. delta(ii,2) .lt. 0.) then
            deb = 1
         else if (delta(ii,1) .lt. 0. .and. delta(ii,2) .ge. 0.) then
            deb = nxy2+2
         else
            deb = 0
         end if

         w2(1) = (1-delx)*(1-dely)
         w2(2) = delx*(1-dely)
         w2(3) = dely*(1-delx)
         w2(4) = delx*dely

         lettre1 = char(48+int(ii/100))
         i2 = ii - int(ii/100)*100
         lettre2 = char(48+int(i2/10))
         lettre3 = char(48+mod(i2,10))
         fin = lettre1 // lettre2 // lettre3 // '.fits'
         do j=1,nxy2
            do i=1,nxy2
               tamp = 0.
               indice = i+(j-1)*nxy2
               indice2 = deb+(j-1)*(nxy2+2)+i
               do k=1,np
                  t1 = i-(c(k,1)-delta(ii,1))
                  t2 = j-(c(k,2)-delta(ii,2))
                  tamp = tamp + a(k,ii)*exp(-b*((t1**2)+(t2**2)))
               end do
               tamp2 = w2(1)*inter(indice2) + w2(2)*inter(indice2+1) +
     &                 w2(3)*inter(indice2+nxy2+2) +
     &                 w2(4)*inter(indice2+nxy2+3)
               tmp2(indice) = z(ii,1)*tamp2+z(ii,2) + tamp
            end do
         end do

c        TE IMAGE "DEC" IS THE SOLUTION OF THE DECONVOLUTION
c        THE NUMBER OF "DEC" IMAGES = THE NUMBER OF IMAGES AT ENTRANCE
c        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
         results = 'dec' // fin 
         call makeima(results,nxy2,tmp2)
   
         call rlft2(tmp2,stmp2,nxy2,1)
         call produit(tmp2,tmp,stmp2,stmp,nxy2)
         call rlft2(tmp2,stmp2,nxy2,-1)

         do j=1,nxy1
            do i=1,nxy1
               ind1 = nxy2*pas*(j-1)+pas*(i-1)
               indice = (j-1)*nxy1+i
               tamp1 = (sfbis(indice,ii)-g(indice,ii))*
     &                 sqrt(sig3(indice,ii))

               do l=1,pas
                  ind2 = (l-1)*nxy2
                  do k=1,pas
                    tmp2(ind1+k+ind2) = tamp1
                  end do
               end do
            end do
         end do

c        THE IMAGE "RESI" IS THE IMAGE OF RESIDUALS
c        THE NUMBER OF "RESI" IMAGES = THE NUMBER OF IMAGES AT ENTRANCE.
c        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         results = ' '
         results = 'resi' // fin
         call makeima(results,nxy2,tmp2)
         do i=1,n2
            tmp2(i) = tmp2(i)*tmp2(i)
            tmp(i) = s(i,ii)
         end do

         do i=1,2*nxy2
            stmp(i) = ss(i,ii)
         end do
         call rlft2(tmp2,stmp2,nxy2,1)
         call produit(tmp2,tmp,stmp2,stmp,nxy2)
         call rlft2(tmp2,stmp2,nxy2,-1)

c        THE IMAGE "RESI_SM" IS THE IMAGE OF SMOOTHED RESIDUALS
c        THE NUMBER OF "RESI_SM" IMAGES = THE NUMBER OF "XIXI" IMAGES
c        - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         results = ' '
         results = 'resi_sm' // fin
         call makeima(results,nxy2,tmp2)
      end do

c     BACKUP OF RESULTS IN 2 OUTPUT FILES
c     - - - - - - - - - - - - - - - - - - 

      write(2,'(A12,F20.6)') 'minimum   : ',fp
      write(2,'(A12,F16.8)') 'chi carre : ',qi
      write(2,*)
      write(2,'(A46)') '* Value of intencity and coordinates x and y :'
      write(2,*)
      do k=1,np
         write(2,'(A19,I4)') '  - Source number: ',k 
         do i=1,m
            write(2,'(A10,3(F15.6,A1))') '          ',a(k,i),
     &                ' ',
     &               c(k,1)-delta(i,1),' ',
     &               c(k,2)-delta(i,2),' '
            write(3,'(F15.6,A1,$)') a(k,i),' '
        end do
        write(3,*)  
        write(3,'(2F12.6)') c(k,1),c(k,2)
        write(2,*)
      end do
      write(3,'(A55)') '------------------------------------------------
     &-------'
      write(3,'(A55)') '|Initial values for z1, z2 delta1 et delta2     
     &      |'
      write(3,'(A55)') '|nb * (4 values) nb = number of images          
     &      |'
      write(3,'(A55)') '------------------------------------------------
     &-------'
      write(2,'(A37)') '* Values of z1, z2 delta1 and delta2:'
      do i=1,m
         write(2,'(A10,4(F12.6,A1))') '          ',z(i,1),' ',
     &   z(i,2),' ',delta(i,1),' ',delta(i,2),' '
         write(3,'(4F12.6)') z(i,1),z(i,2),delta(i,1),delta(i,2)
      end do

      results = 'back' // '.fits'
      call makeima(results,nxy2,p)


      return
      end

c     ----------------------
c     SUBROUTINE: residu
c     ----------------------

      subroutine residu()

      implicit none

      include 'param_dec.f'

      integer nb
      parameter (nb=(xy2+2)*(xy2+2))

      real*4 g(nbr1,nbg),sig2(nbr1,nbg),s(nbr2,nbg),ss(2*xy2,nbg)
      real*4 z(nbg,2),p(nbr2),a(npic,nbg),c(npic,2),b,delta(nbg,2)
      real*4 w(4),delx,dely,sfbis(nbr1,nbg),shfz(nbr2+4*(xy2+1),nbg)
      real*4 inter(nb),sig3(nbr1,nbg)
      real*4 inter2(nb),inter3(nb),fg(nbr1)
      real*4 resim(nbr2),resismm(nbr2),fonddec(nbr2),sfonddec(2*xy2)
      real*4 tmp(nbr2),tmp2(nbr2),tmp3(nbr2),stmp3(2*xy2),stmp2(2*xy2)
      integer nxy1,nxy2,pas,n1,n2,np,m
      integer ii,i,j,indice,indice2,deb,ind1,ind3,ind2,k,l,i2
      character*20 results
      character*1 question,lettre1,lettre2,lettre3
      character*8 fin

      common /param/ b,nxy1,nxy2,pas
      common /gr1/ g,s,sig2,ss
      common /sig/ sig3
      common /tailles/ n1,n2,np,m
      common /variables/ p,z,a,c,delta
      common /gr5/ shfz,sfbis

      do i=1,n2
         resim(i) = 0.
         resismm(i) = 0.
      end do

c     demande de mise à zero des valeurs négatives
c     - - - - - - - - - - - - - - - - - - - - - - -
      question = 'a'
      do while (question.ne.'y' .and. question .ne.'n')
         write(*,'(''Voulez-vous mettre à zéro les valeurs '',$)')
         write(*,'(''négatives du fond (y ou n) : '',$)')
         read(5,*) question
      end do
      if (question .eq. 'y') then
         do i=1,n2
            if (p(i) .lt. 0.) p(i) = 0.
         end do
      end if

c     inter 3 = image de f avec un bord de 1 pixel supplémentaire
c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      do j=1,nxy2
         inter3(j*(nxy2+2)+1) = p(j*nxy2)
         do i=1,nxy2
            indice = (j-1)*nxy2+i
            indice2 = j*(nxy2+2)+1+i
            inter3(indice2) = p(indice)
         end do
         inter3((j+1)*(nxy2+2)) = p((j-1)*nxy2+1)
      end do
      do i=1,nxy2+2
         inter3(i) = inter3((nxy2+2)*nxy2+i)
         inter3((nxy2+2)*(nxy2+1)+i) = inter3(nxy2+2+i)
      end do
      results = 'inter3' // '.fits'
      call makeima(results,nxy2+2,inter3)


      do ii=1,m

         do j=1,nxy1
            do i=1,nxy1
               ind1 = nxy2*pas*(j-1)+pas*(i-1)
               indice = (j-1)*nxy1+i
               do l=1,pas
                  ind2 = (l-1)*nxy2
                  do k=1,pas
                     tmp(ind1+k+ind2) = (sfbis(indice,ii)-g(indice,ii))*
     &                 sqrt(sig3(indice,ii))
                  end do
               end do
            end do
         end do
 
         do i=1,n2
            tmp2(i) = tmp(i)*tmp(i)
            tmp3(i) = s(i,ii)
            fonddec(i) = 0.
         end do

         do i=1,2*nxy2
            stmp3(i) = ss(i,ii)
         end do
         call rlft2(tmp2,stmp2,nxy2,1)
         call produit(tmp2,tmp3,stmp2,stmp3,nxy2)
         call rlft2(tmp2,stmp2,nxy2,-1)


         if (delta(ii,1) .lt. 0.) then
            delx = 1+delta(ii,1)
         else
            delx = delta(ii,1)
         end if
         if (delta(ii,2) .lt. 0.) then
            dely = 1+delta(ii,2)
         else
            dely = delta(ii,2)
         end if

         w(1) = (1-delx)*(1-dely)
         w(2) = delx*(1-dely)
         w(3) = dely*(1-delx)
         w(4) = delx*dely

         do j=1,nxy2
            inter(j*(nxy2+2)+1) = tmp(j*nxy2)
            inter2(j*(nxy2+2)+1) = tmp2(j*nxy2)
            do i=1,nxy2
               indice = (j-1)*nxy2+i
               indice2 = j*(nxy2+2)+1+i
               inter(indice2) = tmp(indice)
               inter2(indice2) = tmp2(indice)
            end do
            inter((j+1)*(nxy2+2)) = tmp((j-1)*nxy2+1)
            inter2((j+1)*(nxy2+2)) = tmp2((j-1)*nxy2+1)
         end do
         do i=1,nxy2+2
            inter(i) = inter((nxy2+2)*nxy2+i)
            inter2(i) = inter2((nxy2+2)*nxy2+i)
            inter((nxy2+2)*(nxy2+1)+i) = inter(nxy2+2+i)
            inter2((nxy2+2)*(nxy2+1)+i) = inter2(nxy2+2+i)
         end do

         results = 'inter' // '.fits'
         call makeima(results,nxy2+2,inter)
         results = 'inter2' // '.fits'
         call makeima(results,nxy2+2,inter2)

         if (delta(ii,1) .ge. 0. .and. delta(ii,2) .ge. 0.) then
            deb = nxy2+3
         else if (delta(ii,1) .ge. 0. .and. delta(ii,2) .lt. 0.) then
            deb = 1
         else if (delta(ii,1) .lt. 0. .and. delta(ii,2) .ge. 0.) then
            deb = nxy2+2
         else
            deb = 0
         end if 

         do j=1,nxy2
            do i=1,nxy2
               ind1 = deb+(nxy2+2)*(j-1)+(i-1)
               indice = (j-1)*nxy2+i
               do l=1,2
                  ind2 = (l-1)*(nxy2+2)
                  ind3 = 2*(l-1)
                  do k=1,2
                    resim(indice) = resim(indice) + 
     &                  w(ind3+k)*inter(ind1+k+ind2) 
                    resismm(indice) = resismm(indice) + 
     &                  w(ind3+k)*inter2(ind1+k+ind2) 
                    fonddec(indice) = fonddec(indice) + 
     &                  w(ind3+k)*inter3(ind1+k+ind2) 
                  end do
               end do
            end do
         end do

c        convolution avec s du fond interpoler
c        - - - - - - - - - - - - - - - - - - - 
         call rlft2(fonddec,sfonddec,nxy2,1)
         call produit(fonddec,tmp3,sfonddec,stmp3,nxy2)
         call rlft2(fonddec,sfonddec,nxy2,-1)

c        mise en gros pixel du fond interpolé et convoluté
c        - - - - - - - - - - - - - - - - - - - - - - - - - 
         do j=1,nxy1
            do i=1,nxy1
               ind1 = nxy2*pas*(j-1)+pas*(i-1)
               indice = (j-1)*nxy1+i
               fg(indice) = 0.
               do l=1,pas
                  ind2 = (l-1)*nxy2
                  do k=1,pas
                     fg(indice) = fg(indice) + fonddec(ind1+k+ind2) 
                  end do
               end do
               fg(indice) = fg(indice)/real(pas*pas)
            end do
         end do

         lettre1 = char(48+int(ii/100))
         i2 = ii - int(ii/100)*100
         lettre2 = char(48+int(i2/10))
         lettre3 = char(48+mod(i2,10))
         fin = lettre1 // lettre2 // lettre3 // '.fits'
         results = 'fg' // fin
         call makeima(results,nxy1,fg)
      end do

      results = 'resi_moyenne' // '.fits'
      call makeima(results,nxy2,resim)
      results = 'resi_sm_moyenne' // '.fits'
      call makeima(results,nxy2,resismm)

      return
      end


c     ----------------------------
c     SUBROUTINE: nouvelle fourier
c     ----------------------------

      subroutine rlft2(data,speq,nn1,isign)
 
      integer isign,nn1
      reaL*4 data(nn1*nn1)
      integer i1,i2,j1,j2,nn(2),i,ind1,ind2,nb
      double precision theta,wi,wpi,wpr,wr,wtemp
      real*4 c1,c2,h1r,h1i,h2r,h2i
      real*4 speq(nn1*2)

      c1 = 0.5
      c2 = -0.5*isign
      theta = 6.28318530717959d0/dble(isign*nn1)
      wpr = -2.0d0*sin(0.5d0*theta)**2
      wpi = sin(theta)
      nn(1) = nn1/2
      nn(2) = nn1
      nb = nn(1)*nn(2)
      if (isign .eq. 1) then
         call fourn(data,nn,nb,2,isign)
        do i=1,nn1
            speq(i*2-1) = data((i-1)*nn1+1)
            speq(i*2) = data((i-1)*nn1+2)
        end do
      end if
      wr = 1.0d0
      wi = 0.0d0
      do i2=1,nn1/2+1,2
         j2 = nn1-i2+2
         do i1=1,nn1
            j1 = 1
            if (i1 .ne. 1) j1 = nn1-i1+2
            if (i2 .eq. 1) then
               h1r = c1*(data((i1-1)*nn1+1)+speq(j1*2-1))
               h1i = c1*(data((i1-1)*nn1+2)-speq(j1*2))
               h2i = c2*(data((i1-1)*nn1+1)-speq(j1*2-1))
               h2r = -c2*(data((i1-1)*nn1+2)+speq(j1*2))
               data((i1-1)*nn1+1) = h1r+h2r
               data((i1-1)*nn1+2) = h1i+h2i
               speq(j1*2-1) = h1r-h2r
               speq(j1*2) = -h1i+h2i
            else
               ind1 = (i1-1)*nn1+i2
               ind2 = (j1-1)*nn1+j2
               h1r = c1*(data(ind1)+data(ind2))
               h1i = c1*(data(ind1+1)-data(ind2+1))
               h2i = c2*(data(ind1)-data(ind2))
               h2r = -c2*(data(ind1+1)+data(ind2+1))
               data(ind1) = h1r+wr*h2r-wi*h2i
               data(ind1+1) = h1i+wr*h2i+wi*h2r
               data(ind2) = h1r-wr*h2r+wi*h2i
               data(ind2+1) = -h1i+wr*h2i+wi*h2r
           end if
         end do
         wtemp = wr
         wr = wr*wpr-wi*wpi+wr
         wi = wi*wpr+wtemp*wpi+wi
      end do
      if (isign .eq. -1) then
         call fourn(data,nn,nb,2,isign)
      end if

      return
      end 

c     ---------------------
c     SUBROUTINE: fourier
c     ---------------------
c     ROUTINE WHICH CALCULATES THE FOURIER TRANSFORM OR INVERSE FOURIER
c     TRANSFORM
c     ---------------------

      subroutine fourn(data,nn,nb,ndim,isign)

c      real*8 wr,wi,wpr,wpi,wtemp,theta
c      real*4 data(nb*2)
c      integer nn(ndim)
      integer nn(ndim),isign,ndim,nb
      double precision wr,wi,wpr,wpi,wtemp,theta
      real*4 data(nb*2),tempi,tempr
      integer i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3
      integer k1,k2,n,nprev,nrem,ntot


      ntot = 1
      do idim=1,ndim
         ntot = ntot*nn(idim)
      end do    
      nprev = 1
      do idim=1,ndim
         n = nn(idim)
         nrem = ntot/(n*nprev)
         ip1 = 2*nprev
         ip2 = ip1*n
         ip3 = ip2*nrem
         i2rev = 1
         do i2=1,ip2,ip1
            if (i2 .lt. i2rev) then
               do i1=i2,i2+ip1-2,2
                  do i3=i1,ip3,ip2
                     i3rev = i2rev+i3-i2
                     tempr = data(i3)
                     tempi = data(i3+1)
                     data(i3) = data(i3rev)
                     data(i3+1) = data(i3rev+1)
                     data(i3rev) = tempr
                     data(i3rev+1) = tempi
                 end do
               end do
            end if
            ibit = ip2/2
30          if ((ibit .ge. ip1) .and. (i2rev .gt. ibit)) then
               i2rev = i2rev-ibit
               ibit = ibit/2
               go to 30
            end if
            i2rev = i2rev+ibit 
        end do
        ifp1 = ip1
31      if (ifp1 .lt. ip2) then
           ifp2 = 2*ifp1
           theta = isign*6.28318530717959D0/(ifp2/ip1)
           wpr = -2.D0*sin(0.5D0*theta)**2
           wpi = sin(theta)
           wr = 1.D0
           wi = 0.D0
           do i3=1,ifp1,ip1
              do i1=i3,i3+ip1-2,2
                do i2=i1,ip3,ifp2
                   k1 = i2
                   k2 = k1+ifp1
                   tempr = sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
                   tempi = sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
                   data(k2) = data(k1)-tempr
                   data(k2+1) = data(k1+1)-tempi
                   data(k1) = data(k1)+tempr
                   data(k1+1) = data(k1+1)+tempi
                end do
              end do
              wtemp = wr
              wr = wr*wpr-wi*wpi+wr
              wi = wi*wpr+wtemp*wpi+wi
           end do
           ifp1 = ifp2
           go to 31
        end if
        nprev = n*nprev
      end do
      
      return
      end

c     ---------------------
c     SUBROUTINE: produit
c     ---------------------
c     ROUTINE WHICH COMPUTES THE PRODUCT OF 2 COMPLEX TABLES
c     ---------------------

      subroutine produit(f1,f2,s1,s2,n)

      integer n,j
      real*4 f1(n*n),f2(n*n),s1(2*n),s2(2*n),sc,r,i
  

      sc = 2./(n*n)
      do j=1,n*n,2
         r = f1(j)*f2(j) - f1(j+1)*f2(j+1)
         i = f1(j)*f2(j+1) + f1(j+1)*f2(j)
         f1(j) = r*sc
         f1(j+1) = i*sc
      end do
      do j=1,2*n,2
         r = s1(j)*s2(j) - s1(j+1)*s2(j+1)
         i = s1(j)*s2(j+1) + s1(j+1)*s2(j)
         s1(j) = r*sc
         s1(j+1) = i*sc
      end do

      return
      end

c     ---------------------
c     SUBROUTINE: openima
c     ---------------------

      subroutine openima(name,xpix,ypix,data)

      implicit none

      integer xpix,ypix
      integer status,unit,readwrite,blocksize,group
      integer nfound_in,naxes(2),npix,fpix

      real*4 data(xpix*ypix)
      real*4 null

      logical anynull

      character name*20

      status=0     
      call ftgiou(unit,status)
      readwrite=0
      blocksize=2880
      call ftopen(unit,name,readwrite,blocksize,status)
      call ftgknj(unit,'NAXIS',1,2,naxes,nfound_in,status)
      npix=naxes(1)*naxes(2)
      xpix=naxes(1)
      ypix=naxes(2)
      group=1
      fpix=1
      null=0.
      call ftgpve(unit,group,fpix,npix,null,data,anynull,status)
      call ftclos(unit,status)
      call ftfiou(unit,status)   
      
      if (status.gt.0) call printerror(status)

      return
      end

c     ---------------------
c     SUBROUTINE: makeima
c     ---------------------

      subroutine makeima(name,xpix,data)
 
      implicit none

      integer status,unit,blocksize
      integer xpix,fpixel,nelement
      integer bitpix
      integer naxes(2)
      integer group
      real*4 data(xpix*xpix)
      logical extend,simple
      character name*20

      status=0
      call delfil(name,status)
      call ftgiou(unit,status)
      blocksize=2880
      call ftinit(unit,name,blocksize,status)
      simple=.true.
      extend=.true.
      bitpix=-32
      naxes(1)=xpix
      naxes(2)=xpix
      call ftphpr(unit,simple,bitpix,2,naxes,0,1,extend,status)
      if (status.gt.0) then
         call printerror(status)
         stop
      end if
      group = 1
      fpixel = 1
      nelement = xpix*xpix
      call ftppre(unit,group,fpixel,nelement,data,status)
      call ftclos(unit,status)
      call ftfiou(unit,status)
      if (status.gt.0) then
         call printerror(status)
         stop
      end if

      return
      end

c     ------------------------
c     SUBROUTINE: printerror
c     ------------------------

      subroutine printerror(status)


      integer status
      character errtext*30,errmessage*80

      if (status .le. 0) return

      call ftgerr(status,errtext)
      print *,'FITSIO Error Status =',status,': ',errtext

 1    call ftgmsg(errmessage)
      do while (errmessage .ne. ' ')
          print *,errmessage
 2        call ftgmsg(errmessage)
      end do

      return
      end

c     --------------------
c     SUBROUTINE: delfil
c     --------------------

      subroutine delfil(filename,status)


      integer status,unit,blocksize
      character*(*) filename
 
      if (status .gt. 0)return

 1    call ftgiou(unit,status)

 2    call ftopen(unit,filename,1,blocksize,status)

      if (status .eq. 0)then
 3        call ftdelt(unit,status)
      else if (status .eq. 103)then
          status=0
 4        call ftcmsg
      else
          status=0
 5        call ftcmsg
          call ftdelt(unit,status)
      end if

 6    call ftfiou(unit, status)
      
      return
      end
