      program deconvolution 
*        (test)
*     *********************
*
      implicit none
*
      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		  MaxIter
*
      real*4 madrid(1)
*
      integer		  Nima,Npsf
      real*8  		  ImaC(maxdim)
      real*8		  PsfC(maxdim)
      real*8		  SigC(maxdim)
      real*8		  lambda
*
      common /vmr/madrid
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      call stspro('deconvolution')
*
      open(unit=1,status='old',file='decon.dat')
*
      read (1,*)name1
      read (1,*)name2
      read (1,*)name3
      read (1,*)name4
      read (1,*)lambda
      read(1,*)MaxIter
*
      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),lambda,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,lambda,
     * MaxIter)
*     *********************************************
*
      implicit none
*
      integer             NpIma,NpPsf
      real                Ima(NpIma)
      real                Psf(NpPsf)
      real                Mod(NpIma)
      real                Sig(NpIma)
*
      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
*
c     parameter(MaxIter=2000)
*
      integer             i,iter
      real*8              Khi2,Khi2s,sf,lambda,se,pf,pg,s
      real*8              p(maxdim)
      real*8		  Temp1(maxdim),Temp2(maxdim)
      real*8		  Temp3(maxdim)
      real*8		  Temp4(maxdim),Temp5(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
        p(i)=ImaC(i)
    5 continue
      do 6 iter=1,MaxIter
        call Convolve(p,Temp1)
        do 7 i=1,Nima
          Temp2(i)=Temp1(i)/SigC(i)**2
          Temp4(i)=ImaC(i)/SigC(i)**2
    7   continue
        call Convolve(Temp2,Temp3)
        call Convolve(Temp4,Temp5)
        sf=0.
        do 13 i=1,Nima
          sf=sf+p(i)
   13   continue
        do 8 i=1,Nima
          se=0.
          if(p(i).gt.0.)se=1.+log(p(i)/sf)
          p(i)=p(i)*(Temp5(i)-0.5*Lambda*se)/Temp3(i)
    8   continue
        s=0.d0
        do 18 i=1,Nima
          s=s+((ImaC(i)-Temp1(i))/SigC(i))**2
   18   continue
        Khi2=s/float(Nima)
        if(Khi2.lt.1.0)Lambda=1.05*Lambda
c        Khi2=Chi2(p)
        pf=int(iter/50.)
        pg=iter/50.
        if(abs(pf-pg).lt.0.01)write(*,101)iter,Khi2,Lambda
  101   format(' Iter:',I4,' Chi2:',f7.3,' Lambda: ',e10.3)
c        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 none
*
      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 none
*
      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
