      program deconvolution 
*        (Gold's method)
*     *********************
*
      implicit undefined (a-z)
*
      integer             dim,maxdim
*
      parameter (dim=1,maxdim=64)
*
      character*8         name1,name2,name3,name4
      integer             npix1(dim)
      integer             npix2(dim)
      integer             npix3(dim)
      integer             npix4(dim)
      real*8  		  start1(dim)
      real*8 		  start2(dim)
      real*8 		  start3(dim)
      real*8 		  start4(dim)
      real*8  		  step1(dim)
      real*8 		  step2(dim)
      real*8  		  step3(dim)
      real*8  		  step4(dim)
      character*72        ident1,ident2,ident3,ident4
      character*64        cunit1,cunit2,cunit3,cunit4
      integer             pntr1,pntr2,pntr3,pntr4
      integer             stat1,stat2,stat3,stat4
      integer		  naxis1,naxis2,naxis3,naxis4
      integer             no1,no2,no3,no4
      integer		  acts
      integer  		  kun
      integer		  knul
      integer		  k_status
      integer		  MaxIter,np
*
      real*4 madrid(1)
*
      integer		  Nima,Npsf
      real*8  		  ImaC(maxdim)
      real*8		  PsfC(maxdim)
      real*8		  SigC(maxdim)
*
      common /vmr/madrid
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      call stspro('deconvolution')
*
*     Lecture des noms de fichiers
*
      call stkrdc('IN_A',1,1,8,acts,name1,kun,knul,k_status)
      call stkrdc('IN_B',1,1,8,acts,name2,kun,knul,k_status)
      call stkrdc('IN_C',1,1,8,acts,name3,kun,knul,k_status)
      call stkrdc('IN_D',1,1,8,acts,name4,kun,knul,k_status)
*
*     Lecture du nombre maximum d'iterations
*
      call stkrdi('MAXINT',1,1,np,MaxIter,kun,knul,k_status)
*
      call stiget(name1,10,0,1,dim,naxis1,npix1,start1,step1,
     +            ident1,cunit1,pntr1,no1,stat1)
*
      call stiget(name2,10,0,1,dim,naxis2,npix2,start2,step2,
     +            ident2,cunit2,pntr2,no2,stat2)
*
      call stiget(name3,10,2,1,dim,naxis3,npix3,start3,step3,
     +            ident3,cunit3,pntr3,no3,stat3)
*
      call stiget(name4,10,0,1,dim,naxis4,npix4,start4,step4,
     +            ident4,cunit4,pntr4,no4,stat4)
*
      call Decon(madrid(pntr1),madrid(pntr2),madrid(pntr3),
     +           madrid(pntr4),npix1(1),npix2(1),MaxIter)
*
      call stfclo(no1,stat1)
      call stfclo(no2,stat2)
      call stfclo(no3,stat3)
      call stfclo(no4,stat4)
*
      close(1)
*
      call stsepi
*
      end
*
*     *****************************************************
      subroutine Decon(Ima,Psf,Mod,Sig,NpIma,NpPsf,MaxIter)
*     *****************************************************
*
      implicit undefined (a-z)
*
      integer             NpIma,NpPsf
      real                Ima(NpIma)
      real                Psf(NpPsf)
      real                Mod(NpIma)
      real                Sig(NpIma)
      real		  f1,f2
*
      integer             dim,maxdim
*
      parameter(dim=1,maxdim=64)
*
      real*8              ImaC(maxdim)
      real*8              PsfC(maxdim)
      real*8              SigC(maxdim)
      integer             Nima,Npsf
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      integer		  MaxIter
*
      integer             i,iter
      real*8              Chi2,Khi2,Khi2s
      real*8              p(maxdim)
      real*8		  Temp1(maxdim),Temp2(maxdim)
*
      Nima=NpIma
      Npsf=NpPsf
      Khi2s=1.d30
      do 1 i=1,maxdim
        ImaC(i)=0.d0
        PsfC(i)=0.d0
        SigC(i)=0.d0
        p(i)=0.d0
    1   continue
      do 2 i=1,Nima
        ImaC(i)=Ima(i)
        p(i)=Mod(i)
        SigC(i)=Sig(i)
    2   continue
      do 3 i=1,Npsf
        PsfC(i)=Psf(i)
    3 continue
      if(ImaC(1).lt.0.d0)then
        ImaC(2)=ImaC(2)+ImaC(1)
        ImaC(1)=0.d0
      endif
      if(ImaC(Nima).lt.0.d0)then
        ImaC(Nima-1)=ImaC(Nima-1)+ImaC(Nima)
        ImaC(Nima)=0.d0
      endif
      do 4 i=2,Nima-1
        if(ImaC(i).lt.0.d0)then
          ImaC(i-1)=ImaC(i-1)+0.5*ImaC(i)
          ImaC(i+1)=ImaC(i+1)+0.5*ImaC(i)
          ImaC(i)=0.d0
        endif
    4 continue
      do 5 i=1,Nima
        if(ImaC(i).lt.0.d0)ImaC(i)=0.d0
    5 continue
*
      if(p(1).lt.0.d0)then
        p(2)=p(2)+p(1)
        p(1)=0.d0
      endif
      if(p(Nima).lt.0.d0)then
        p(Nima-1)=p(Nima-1)+p(Nima)
        p(Nima)=0.d0
      endif
      do 14 i=2,Nima-1
        if(p(i).lt.0.d0)then
          p(i-1)=p(i-1)+0.5*p(i)
          p(i+1)=p(i+1)+0.5*p(i)
          p(i)=0.d0
        endif
   14 continue
      do 15 i=1,Nima
        if(p(i).lt.0.d0)p(i)=0.d0
   15 continue
*
*     Iteration
*
      do 6 iter=1,MaxIter
        call Convolve(p,Temp1)
        do 7 i=1,Nima
          Temp2(i)=0.d0
          if(Temp1(i).gt.0.d0)Temp2(i)=ImaC(i)/Temp1(i)
    7   continue
        do 8 i=1,Nima
          p(i)=p(i)*Temp2(i)
    8   continue
        Khi2=Chi2(p)
*
        f1=iter/1.
        f2=iter/1
        if(f1.eq.f2)write(*,101)iter,Khi2
  101   format(' Iter:',I4,' Chi2:',f7.3)
        if(Khi2.lt.1.0d0)goto 9
c        if(Khi2.gt.Khi2s)goto 9
        Khi2s=Khi2
    6  continue
    9 do 12 i=1,Nima
        Mod(i)=p(i)
   12 continue
      return
      end
*
*     ****************************
      subroutine Convolve(Mod,Res)
*     ****************************
*
      implicit undefined (a-z)
*
      integer MaxDim
*
      parameter (maxdim=64)
*
      real*8  Mod(maxdim),Res(maxdim)
*
      real*8              ImaC(maxdim)
      real*8              PsfC(maxdim)
      real*8              SigC(maxdim)
      integer             Nima,Npsf
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      integer        i,ics,ics1,ip,imin,imax,ii
      real*8         SC
*
      ics=(Npsf-1)/2
      ics1=ics+1
      do 1 i=1,Nima
        Res(i)=Mod(i)
    1 continue
      do 15 i=1,Nima
        imin=max(1,i-ics)
        imax=min(i+ics,Nima)
        SC=0.d0
        ip=ics1-i
        do 13 ii=imin,imax
          SC=SC+Mod(ii)*PsfC(ii+ip)
   13   continue
        Res(i)=SC
   15 continue
      return
      end
*
*     *************************
      real*8 function Chi2(Mod)
*     *************************
*
      implicit undefined (a-z)
*
      integer MaxDim
*
      parameter (maxdim=64)
*
      real*8  Mod(maxdim),Res(maxdim)
*
      real*8              ImaC(maxdim)
      real*8              PsfC(maxdim)
      real*8              SigC(maxdim)
      integer             Nima,Npsf
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      integer i
      real*8  s
*
      call Convolve(Mod,Res)
      s=0.d0
      do 1 i=1,Nima
        s=s+((ImaC(i)-Res(i))/SigC(i))**2
    1 continue
      Chi2=s/float(Nima)
      return
      end
