      program deconvolution
*     *********************
*
*     Chi2 minimisation
*
      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		  acts
      integer  		  kun
      integer		  knul
      integer		  k_status
*
      real*4 madrid(1)
*
      integer		  Nima,Npsf
      real*4  		  ImaC(maxdim)
      real*4		  PsfC(maxdim)
      real*4		  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)
*
      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))
*
      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)
*     *********************************************
*
      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*4              ImaC(maxdim)
      real*4              PsfC(maxdim)
      real*4              SigC(maxdim)
      integer             Nima,Npsf
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      integer		  MaxIter
      real*4	          FTol
*
      parameter(MaxIter=100,FTol=1.0d-5)
*
      integer             i,j,ibig,iter
      real*4              Func,Chi2,Khi2
      real*4              p(maxdim),xi(maxdim,maxdim)
      real*4              del,t,Fret,fp,fptt
      real*4              pt(maxdim),ptt(maxdim),xit(maxdim)
*
      Nima=NpIma
      Npsf=NpPsf
      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
      Fret=Func(p)
      do 5 j=1,Nima
        pt(j)=p(j)
    5 continue
      do 21 j=1,Nima
        do 20 i=1,Nima
          xi(i,j)=0.d0
   20   continue
   21 continue
      do 22 i = 1,Nima
        xi(i,i)=1.0d-1*(1.0d-1+abs(p(i)))
   22 continue
      iter=0
    6 iter=iter+1
      Khi2=Chi2(p)
*
      write(*,101)iter,Fret,Khi2
  101 format(' Iter:',I3,' Fret:',e14.7,' Chi2:',f7.3)
      fp=Fret
      ibig=0
      del=0.d0
      do 8 i=1,Nima
        do 7 j=1,Nima
          xit(j)=xi(j,i)
    7   continue
        call LinMin(p,xit,Nima,Fret)
        if(abs(fp-Fret).gt.del)then
          del=abs(fp-Fret)
          ibig=i
        endif
    8 continue
      if(2.0d0*abs(fp-Fret).le.FTol*(abs(fp)+abs(Fret))) go to 11
      if(iter.eq.MaxIter) go to 11
      do 9 j=1,Nima
        ptt(j)=2.*p(j)-pt(j)
        xit(j)=p(j)-pt(j)
        pt(j)=p(j)
    9 continue
      fptt=Func(ptt)
      if(fptt.ge.fp) go to 6
      t=2.0d0*(fp-2.0d0*Fret+fptt)*(fp-Fret-del)**2-del*(fp-fptt)**2
      if(t.ge.0.d0) go to 6
      call LinMin(p,xit,Nima,Fret)
      do 10 j=1,Nima
        xi(j,ibig)=xit(j)
   10 continue
      go to 6
   11 do 12 i=1,Nima
        Mod(i)=p(i)
   12 continue
      return
      end
*
*     ******************************
      subroutine LinMin(p,xi,n,Fret)
*     ******************************
*
      implicit none
*
      integer n,maxdim
      real*4 Tol
*
      parameter(maxdim=64,Tol=1.0d-06)
*
      integer ncom,j
      real*4  pcom(maxdim),xicom(maxdim)
      real*4 ax,bx,xx,p(n),xi(n),Fret,fa,fb,fx,xmin,Brent
*
      common/F1Com/ncom,pcom,xicom
*
      ncom=n
      do 1 j=1,n
        pcom(j)=p(j)
        xicom(j)=xi(j)
    1 continue
      ax=0.0d0
      xx=1.0d0
      bx=2.0d0
      call MnBrak(ax,xx,bx,fa,fx,fb)
      Fret=Brent(ax,xx,bx,Tol,xmin)
      do 2 j=1,n
        xi(j)=xmin*xi(j)
        p(j)=p(j)+xi(j)
    2 continue
      return
      end
*
*     ************************
      real*4 function F1Dim(x)
*     ************************
*
      implicit none
*
      integer maxdim
*
      parameter(maxdim=64)
*
      integer ncom
      real*4  pcom(maxdim),xicom(maxdim)
*
      common/F1Com/ncom,pcom,xicom
*
      integer j
      real*4 x,xt(maxdim),Func
*
      do 1 j=1,ncom
        xt(j)=pcom(j)+x*xicom(j)
    1 continue
      F1Dim=Func(xt)
      return
      end
*
*     ***********************
      real*4 function Func(p)
*     ***********************
*
      implicit none
*
      integer             dim,maxdim
*
      parameter(dim=1,maxdim=64)
*
      integer             Nima,Npsf
      real*4		  ImaC(maxdim)
      real*4              PsfC(maxdim)
      real*4              SigC(maxdim)
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      integer             i
      real*4              p(maxdim),Conv(maxdim)
      real*4              s,sf
*
      call Convolve(p,Conv)
      s=0.d0
      sf=0.d0
      do 2 i=1,Nima
        sf=sf+p(i)
    2 continue
      do 1 i=1,Nima
        s=s+((Conv(i)-ImaC(i))/SigC(i))**2
    1 continue
      s=0.5d0*s
      Func=s
      return
      end
*
*     ************************************
      subroutine MnBrak(ax,bx,cx,fa,fb,fc)
*     ************************************
*
      implicit none
*
      real*4 gold,glimit,tiny,ax,bx,cx,fa,fb,fc,dum,r,q,u
      real*4 ulim,fu
      real*4 F1Dim
*
      parameter (gold=1.618034d0,glimit=100.0d0,tiny=1.0d-06)
*
      fa=F1Dim(ax)
      fb=F1Dim(bx)
      if(fb.gt.fa)then
        dum=ax
        ax=bx
        bx=dum
        dum=fb
        fb=fa
        fa=dum
      endif
      cx=bx+gold*(bx-ax)
      fc=F1Dim(cx)
    1 if(fb.ge.fc)then
        r=(bx-ax)*(fb-fc)
        q=(bx-cx)*(fb-fa)
        u=bx-((bx-cx)*q-(bx-ax)*r)/(2.*sign(max(abs(q-r),tiny),q-r))
        ulim=bx+glimit*(cx-bx)
        if((bx-u)*(u-cx).gt.0.)then
          fu=F1Dim(u)
          if(fu.lt.fc)then
            ax=bx
            fa=fb
            bx=u
            fb=fu
            go to 1
          else if(fu.gt.fb)then
            cx=u
            fc=fu
            go to 1
          endif
          u=cx+gold*(cx-bx)
          fu=F1Dim(u)
        else if((cx-u)*(u-ulim).gt.0.)then
          fu=F1Dim(u)
          if(fu.lt.fc)then
            bx=cx
            cx=u
            u=cx+gold*(cx-bx)
            fb=fc
            fc=fu
            fu=F1Dim(u)
          endif
        else if((u-ulim)*(ulim-cx).ge.0.)then
          u=ulim
          fu=F1Dim(u)
        else
          u=cx+gold*(cx-bx)
          fu=F1Dim(u)
        endif
        ax=bx
        bx=cx
        cx=u
        fa=fb
        fb=fc
        fc=fu
        go to 1
      endif
      return
      end
*
*     ****************************************
      real*4 function Brent(ax,bx,cx,tol,xmin)
*     ****************************************
*
      implicit none
*
      integer itmax,iter
      real*4 cgold,zeps,ax,bx,cx,tol,xmin,a,b,v,w,x,e,fx,fv,fw
      real*4 xm,tol1,tol2,r,q,p,etemp,d,u,fu
      real*4 F1Dim
*
      parameter (itmax=200,cgold=.381966d0,zeps=1.0d-05)
*
      a=min(ax,cx)
      b=max(ax,cx)
      v=bx
      w=v
      x=v
      e=0.
      fx=F1Dim(x)
      fv=fx
      fw=fx
      do 11 iter=1,itmax
        xm=0.5*(a+b)
        tol1=tol*abs(x)+zeps
        tol2=2.*tol1
        if(abs(x-xm).le.(tol2-0.5*(b-a)))goto 3
        if(abs(e).gt.tol1)then
          r=(x-w)*(fx-fv)
          q=(x-v)*(fx-fw)
          p=(x-v)*q-(x-w)*r
          q=2.*(q-r)
          if(q.gt.0.) p=-p
          q=abs(q)
          etemp=e
          e=d
          if(abs(p).ge.abs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x))
     +      goto 1
          d=p/q
          u=x+d
          if(u-a.lt.tol2.or.b-u.lt.tol2)d=sign(tol1,xm-x)
          goto 2
        endif
    1 if(x.ge.xm)then
        e=a-x
      else
        e=b-x
      endif
      d=cgold*e
    2 if(abs(d).ge.tol1)then
        u=x+d
      else
        u=x+sign(tol1,d)
      endif
      fu=F1Dim(u)
      if(fu.le.fx)then
        if(u.ge.x)then
          a=x
        else
          b=x
        endif
        v=w
        fv=fw
        w=x
        fw=fx
        x=u
        fx=fu
      else
        if(u.lt.x)then
          a=u
        else
          b=u
        endif
        if(fu.le.fw.or.w.eq.x)then
          v=w
          fv=fw
          w=u
          fw=fu
        else if(fu.le.fv.or.v.eq.x.or.v.eq.w)then
          v=u
          fv=fu
        endif
      endif
   11 continue
c      pause 'Brent exceed maximum iterations'
    3 xmin=x
      Brent=fx
      return
      end
*
*     **********************
      real*4 function Erf(x)
*     **********************
*
      implicit none
*
      real*4      t,s,x
*
      s=1.d0
      if(x.lt.5)then
        t=1.0d0/(1.0d0+0.47046*x)
        s=t*(0.3480242+t*(-0.0958798+0.7478556*t))
        s=1.0d0-s*Exp(-x*x)
      endif
      Erf=s
*
      return
      end
*
*     ****************************
      subroutine Convolve(Mod,Res)
*     ****************************
*
      implicit none
*
      integer MaxDim
*
      parameter (maxdim=64)
*
      real*4  Mod(maxdim),Res(maxdim)
*
      real*4              ImaC(maxdim)
      real*4              PsfC(maxdim)
      real*4              SigC(maxdim)
      integer             Nima,Npsf
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      integer        i,ics,ics1,ip,imin,imax,ii
      real*4         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*4 function Chi2(Mod)
*     *************************
*
      implicit none
*
      integer MaxDim
*
      parameter (maxdim=64)
*
      real*4  Mod(maxdim),Res(maxdim)
*
      real*4              ImaC(maxdim)
      real*4              PsfC(maxdim)
      real*4              SigC(maxdim)
      integer             Nima,Npsf
*
      Common/CF/ImaC,PsfC,SigC,Nima,Npsf
*
      integer i
      real*4  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
