*     **************
      program ajuste
*     **************

      implicit none

*     -------------------------
*     DECLARATION DES VARIABLES
*     -------------------------

*     Nombre maximum de points : max_data
*     Nombre maximum d'iterations : max_iter

      integer		max_data,max_iter

      parameter		(max_data=2000)
      parameter		(max_iter=100)

*     Nombre total de parametres : ma 

      integer		ma

      parameter		(ma=2)

      real*4		madrid(1)
      integer		naxis
      integer		npix(1)
      integer		pntr_in
      integer		no_in
      integer		pntr_sig
      integer		no_sig
      integer		pntr_out
      integer		no_out
      integer		acts
      integer		kun
      integer		knul
      integer		k_status
      integer		i_status
      integer		np

      real		p(ma)

      real*8		start(1)
      real*8		step(1)

      character*8	in_ima
      character*8	in_sig
      character*8	out_ima
      character*64	cunit
      character*72	ident

      common/vmr/madrid

      call stspro('ajuste')

*     --------------------------- 
*     Lecture des noms de fichier
*     --------------------------- 
 
      call stkrdc('IN_1',1,1,8,acts,in_ima,kun,knul,k_status)
      call stkrdc('IN_2',1,1,8,acts,in_sig,kun,knul,k_status)
      call stkrdc('OUT_A',1,1,8,acts,out_ima,kun,knul,k_status)

*     --------------------------- 
*     Lecture des valeurs initiales des parametres
*     --------------------------- 

      call stkrdr('IN_3',1,ma,np,p,kun,knul,k_status)

*     ----------------------------
*     Lecture du fichier des sigma
*     ----------------------------

      call stiget(in_sig,10,0,1,1,naxis,npix,start,step,
     +            ident,cunit,pntr_sig,no_sig,i_status)

*     --------------------------
*     Lecture du fichier donnees
*     --------------------------

      call stiget(in_ima,10,0,1,1,naxis,npix,start,step,
     +            ident,cunit,pntr_in,no_in,i_status)

*     -------------------------
*     Ajustement de la fonction
*     -------------------------

      call fit(madrid(pntr_in),madrid(pntr_sig),npix,start,step,p,np)

*     ----------------------------
*     Ecriture du fichier resultat
*     ---------------------------- 

      call stiput(out_ima,10,1,1,naxis,npix,start,step,
     +            ident,cunit,pntr_out,no_out,i_status)
      call write_res(madrid(pntr_out),npix,start,step,p)

*     ----------------
*     Fin du programme
*     ----------------

      call stfclo(no_in,i_status)
      call stfclo(no_sig,i_status)
      call stfclo(no_out,i_status)
      call stsepi

      stop
      end


*     ==============================================================
*                SOUS-ROUTINES DU PROGRAMME
*     ==============================================================

*     ===============================================
      subroutine fit(data,sigma,npix,start,step,p,np)
*     ===============================================

      implicit none

*     --------------------------------------
*     Nombre maximum de points : max_data
*     Nombre maximum d'iterations : max_iter
*     --------------------------------------
      integer		max_data,max_iter
      parameter		(max_data=2000)
      parameter		(max_iter=1000)

*     ---------------------------
*     Tolerance pour arret du fit
*     ---------------------------
      real		epsilon
      parameter		(epsilon=0.0001)

*     -------------------------------
*     Nombre total de parametres : ma 
*     -------------------------------
      integer		ma
      parameter		(ma=2)

*     ---------
*     variables 
*     ---------
      integer 		npix
      integer		i,j
      integer		ndata
      integer		np
      integer		mfit
      integer		lista(ma)
      integer		nca
      integer		iter
      integer		atest
      real		x(max_data)
      real		y(max_data)
      real		sig(max_data)
      real		data(npix)
      real		sigma(npix)
      real		p(ma)
      real		covar(ma,ma)
      real		alpha(ma,ma)
      real		chisq
      real		alambda
      real		chisq_old
      real		rel_change
      real		a(ma),beta(ma)
      real		q
      real		ochisq
      real*8		start(1)
      real*8		step(1)

      ochisq = 0.

      do i=1,max_data
         x(i) = 0.
         y(i) = 0.
         sig(i) = 0.
      end do
*     ------------------------------------------------
*     Valeurs de la variable independante x,
*             des donnees a ajuster y(x),
*             des incertitudes sur les donnees sig(x).
*     ------------------------------------------------
      do i=1,npix
        x(i)=start(1)+(i-1)*step(1)
        y(i)=data(i)
        sig(i)=sigma(i)
      end do

*     --------------------------------
*     Valeurs initiales des parametres
*     --------------------------------
      do i=1,np
        a(i)=p(i)
      end do

      if (np.lt.ma) then
         j=np+1
         do i=j,ma
            a(i)=0.
         end do
      end if

*     -------------------------------------
*     Nombre de points : ndata
*     Nombre de parametres a ajuster : mfit 
*     -------------------------------------
      ndata=npix
      mfit=ma

*     ---------------------------------------------------------
*     Numeros des parametres a ajuster : lista(i), i=1,...,mfit
*     ---------------------------------------------------------
      do i=1,mfit
         lista(i)=i
      end do

*     -------------------------------------------------------
*     Dimension des matrices covar(nca,nca) et alpha(nca,nca)
*     -------------------------------------------------------
      nca=ma
      do j=1,nca
         do i=1,nca
            covar(j,i) = 0.
            alpha(j,i) = 0.
         end do
      end do

*     --------------
*     Initialisation
*     --------------
      alambda=-0.001
      chisq = 0.
      atest = 0
      call mrqmin(x,y,sig,ndata,a,lista,mfit,covar,alpha,
     +            nca,chisq,ochisq,alambda,atest,beta)
      write(*,105)

*     -------------------------
*     Ajustement de la fonction
*     -------------------------
      do iter=1,max_iter
         chisq_old=chisq
         call mrqmin(x,y,sig,ndata,a,lista,mfit,covar,
     +                 alpha,nca,chisq,ochisq,alambda,atest,beta)
         rel_change=1.-chisq_old/chisq
         write(*,104)iter,alambda,chisq
         if (rel_change.lt.0..and.rel_change.gt.-epsilon.and.
     +      atest.eq.1.and.iter.gt.2) goto 4
      end do

*     --------------------------------------------------------
*     Arret de l'iteration, calcul de la matrice de covariance
*     --------------------------------------------------------
    4 alambda=0.
      call mrqmin(x,y,sig,ndata,a,lista,mfit,covar,
     +               alpha,nca,chisq,ochisq,alambda,atest,beta)

*     -------------------------------------
*     Accord entre le modele et les donnees
*     -------------------------------------
      call gammq(0.5*(ndata-2),0.5*chisq,q)


*     --------------------------------------------------
*     Valeurs des parametres ajustes et de la covariance
*     --------------------------------------------------
      write(*,100)
      do i=1,ma
         p(i)=a(i)
         write(*,101)i,p(i)
      end do
      write(*,100)
      do i=1,ma
         do j=1,ma
            write(*,102)i,j,covar(i,j)
         end do
      end do
      np=ma

*     ---------------------------
*     Valeur du chi-carre et de Q
*     ---------------------------
      write(*,100)
      write(*,103)chisq,q

*     -------
*     Formats
*     ------- 
  100 format(' ')
  101 format(' a(',I1,') = ',E14.8)
  102 format(' covar(',I1,',',I1,') = ',E12.4)
  103 format(' chi2 = ',F11.3,'   Q = ',E12.4)
  104 format(i5,2e10.3)
  105 format(' iter   lambda    chi2 ')

      return
      end

*     ============================================
      subroutine write_res(resu,npix,start,step,p)
*     ============================================

      implicit none

*     --------------------------------------
*     Nombre maximum de points : max_data
*     Nombre maximum de parametres : max_par
*     --------------------------------------
      integer		max_data
      parameter		(max_data=2000)

*     -------------------------------
*     Nombre total de parametres : ma 
*     -------------------------------
      integer		ma
      parameter		(ma=2)

      integer		npix
      integer		i
      real		x(max_data)
      real		resu(max_data)
      real		p(ma)
      real		dyda(ma)
      real		a(ma)
      real*8		start(1)
      real*8		step(1)

      do i=1,ma
         a(i)=p(i)
      end do

      do i=1,npix
         x(i)=start(1)+(i-1)*step(1)
      end do

      do i=1,npix
         call funcs(x(i),a,resu(i),dyda,ma)
      end do

      return
      end

*     ======================================================
      subroutine mrqmin(x,y,sig,ndata,a,lista,mfit,
     +       covar,alpha,nca,chisq,ochisq,alambda,itest,beta)
*     ======================================================

      implicit	none

*     -----------------------------------
*     Nombre maximum de points : max_data
*     -----------------------------------
      integer		max_data
      parameter		(max_data=2000)

*     -------------------------------
*     Nombre total de parametres : ma 
*     -------------------------------
      integer		ma
      parameter		(ma=2)

      integer 		lista(ma)
      integer		ndata,nca,itest
      integer		j,k,kk,mfit,ihit
      real		chisq,ochisq,alambda
      real 		x(max_data),y(max_data),sig(max_data),a(ma)
      real		covar(nca,nca),alpha(nca,nca),atry(ma)
      real              beta(ma),da(ma)

      if (alambda .lt. 0.) then
         kk=mfit+1
         do j=1,ma
            ihit=0
            do k=1,mfit
              if (lista(k) .eq. j) ihit=ihit+1
            end do
            if (ihit .eq. 0) then
               lista(kk)=j
               kk=kk+1
            else if (ihit .gt. 1) then
               pause 'Improper permutation in LISTA'
            endif
         end do
         if (kk .ne. (ma+1)) pause 'Improper permutation in LISTA'
         alambda=1.
         call mrqcof(x,y,sig,ndata,a,lista,mfit,alpha,beta,nca,chisq)
         ochisq=chisq
         do j=1,ma
            atry(j)=a(j)
         end do
      end if

      do j=1,mfit
         do k=1,mfit
            covar(j,k)=alpha(j,k)
         end do
         covar(j,j)=alpha(j,j)*(1.+alambda)
         da(j)=beta(j)
      end do

      call gaussj(covar,mfit,nca,da)

      if (alambda.eq.0.) then
         call covsrt(covar,nca,ma,lista,mfit)
         return
      end if

      do j=1,mfit
         atry(lista(j))=a(lista(j))+da(j)
      end do

      call mrqcof(x,y,sig,ndata,atry,lista,mfit,covar,da,nca,chisq)

      if (chisq.lt.ochisq) then
         alambda=0.333333*alambda
         itest=1
         ochisq=chisq
         do j=1,mfit
            do k=1,mfit
               alpha(j,k)=covar(j,k)
            end do
            beta(j)=da(j)
            a(lista(j))=atry(lista(j))
         end do
      else
         alambda=3.*alambda
         itest=2
         chisq=ochisq
      end if

      return
      end

*     =========================================================
      subroutine mrqcof(x,y,sig,ndata,a,lista,mfit,alpha,beta,
     +                  nca,chisq)
*     =========================================================

      implicit	none

*     -----------------------------------
*     Nombre maximum de points : max_data
*     -----------------------------------
      integer		max_data
      parameter		(max_data=2000)

*     -------------------------------
*     Nombre total de parametres : ma 
*     -------------------------------
      integer		ma
      parameter		(ma=2)

      integer		lista(ma)
      integer		mfit,nca,ndata
      integer		i,j,k
      real		chisq,ymod,sig2i,dy,wt
      real 		x(max_data),y(max_data),sig(max_data),a(ma)
      real		alpha(nca,nca),beta(ma),dyda(ma)

      ymod = 0.
      do j=1,mfit
         do k=1,j
            alpha(j,k)=0.
         end do
         beta(j)=0.
      end do

      do i=1,ma
         dyda(i) = 0.
      end do

      chisq=0.
      do i=1,ndata
         call funcs(x(i),a,ymod,dyda,ma)
         sig2i=1./(sig(i)*sig(i))
         dy=y(i)-ymod
         do j=1,mfit
            wt=dyda(lista(j))*sig2i
            do k=1,j
               alpha(j,k)=alpha(j,k)+wt*dyda(lista(k))
            end do
            beta(j)=beta(j)+dy*wt
         end do
         chisq=chisq+dy*dy*sig2i
      end do

      do j=2,mfit
         do k=1,j-1
            alpha(k,j)=alpha(j,k)
         end do
      end do

      return
      end


*     ================================
      subroutine gaussj(a,n,np,b)
*     ================================

      implicit none

      integer		nmax
      parameter		(nmax=40)

      integer		n
      integer		np
      integer		ipiv(nmax)
      integer		indxr(nmax)
      integer		indxc(nmax)
      integer		i,j,irow,icol,k,l,ll
      real		a(np,np)
      real		b(np)
      real		big,dum,pivinv

      do j=1,n
         ipiv(j)=0
      end do

      do i=1,n
         big=0.
         do j=1,n
            if (ipiv(j) .ne. 1) then
               do k=1,n
                  if (ipiv(k) .eq. 0) then
                     if (abs(a(j,k)) .ge. big) then
                        big=abs(a(j,k))
                        irow=j
                        icol=k
                     end if
                  else if (ipiv(k) .gt. 1) then
                     pause 'Singular matrix'
                  end if
               end do
            end if
         end do
         ipiv(icol)=ipiv(icol)+1

         if (irow.ne.icol) then
            do l=1,n
               dum=a(irow,l)
               a(irow,l)=a(icol,l)
               a(icol,l)=dum
            end do
               dum=b(irow)
               b(irow)=b(icol)
               b(icol)=dum
         end if
         indxr(i)=irow
         indxc(i)=icol
         if (a(icol,icol) .eq. 0.) pause 'Singular matrix'
         pivinv=1./a(icol,icol)
         a(icol,icol)=1.
         do l=1,n
            a(icol,l)=a(icol,l)*pivinv
         end do
            b(icol)=b(icol)*pivinv
         do ll=1,n
            if (ll .ne. icol) then
               dum=a(ll,icol)
               a(ll,icol)=0.
               do l=1,n
                  a(ll,l)=a(ll,l)-a(icol,l)*dum
               end do
                  b(ll)=b(ll)-b(icol)*dum
            end if
         end do

      end do

      do l=n,1,-1
         if (indxr(l) .ne. indxc(l)) then
            do k=1,n
               dum=a(k,indxr(l))
               a(k,indxr(l))=a(k,indxc(l))
               a(k,indxc(l))=dum
            end do
         end if
      end do

      return
      end


*     ===========================================
      subroutine covsrt(covar,ncvm,ma,lista,mfit)
*     ===========================================

      implicit	none

      integer		ncvm,ma,mfit
      integer		lista(ma)
      integer		i,j
      real		swap  
      real 		covar(ncvm,ncvm)

      do j=1,ma-1
         do i=j+1,ma
            covar(i,j)=0.
         end do
      end do

      do i=1,mfit-1
         do j=i+1,mfit
            if (lista(j) .gt. lista(i)) then
               covar(lista(j),lista(i))=covar(i,j)
            else
               covar(lista(i),lista(j))=covar(i,j)
            end if
         end do
      end do

      swap=covar(1,1)
      do j=1,ma
         covar(1,j)=covar(j,j)
         covar(j,j)=0.
      end do

      covar(lista(1),lista(1))=swap
      do j=2,mfit
         covar(lista(j),lista(j))=covar(1,j)
      end do

      do j=2,ma
         do i=1,j-1
            covar(i,j)=covar(j,i)
         end do
      end do

      return
      end

*     =======================
      subroutine gammq(a,x,q)
*     =======================

*     ------------------------------------------------------
*     Calcul de la fonction gamma incomplete Q(a,x)=1-P(a,x)
*     ------------------------------------------------------

      implicit	none

      real		a
      real		x
      real		q
      real		gamser,gammcf
      real		gln

      if ((x .lt. 0.) .or. (a .le. 0.)) pause
      if (x .lt. a+1.) then
         call gser(gamser,a,x,gln)
         q=1.-gamser
      else
         call gcf(gammcf,a,x,gln)
         q=gammcf
      end if

      return
      end

*     ===============================
      subroutine gser(gamser,a,x,gln)
*     ===============================

      implicit	none

      integer		itmax,n
      real		a
      real		x
      real		gamser
      real		gln
      real		ap,sum,del
      real		eps
      real              gammln

      parameter(itmax=100,eps=3.E-7)

      gln=gammln(a)
      if (x .le. 0.) then
         if (x .lt. 0.) pause
         gamser=0.
         return
      endif

      ap=a
      sum=1./a
      del=sum
      do n=1,itmax
         ap=ap+1.
         del=del*x/ap
         sum=sum+del
         if (abs(del) .lt. abs(sum)*eps) goto 2
      end do

      pause 'A trop grand, ITMAX trop petit'
    2 gamser=sum*exp(-x+a*log(x)-gln)

      return
      end

*     ==============================
      subroutine gcf(gammcf,a,x,gln)
*     ==============================

      implicit	none

      integer		itmax,n
      real		a
      real		x
      real		gammcf
      real		gln
      real		a0,a1,b0,b1,fac,an,ana,anf,gold,g
      real		eps
      real              gammln

      parameter(itmax=100,eps=3.E-7)

      gln=gammln(a)
      gold=0.
      a0=1.
      a1=x
      b0=0.
      b1=1.
      fac=1.
      do n=1,itmax
         an=float(n)
         ana=an-a
         a0=(a1+a0*ana)*fac
         b0=(b1+b0*ana)*fac
         anf=an*fac
         a1=x*a0+anf*a1
         b1=x*b0+anf*b1
         if (a1 .ne. 0.) then
            fac=1./a1
            g=b1*fac
            if (abs((g-gold)/g) .lt. eps) goto 2
            gold=g
         endif
      end do

      pause 'A trop grand, ITMAX trop petit'
    2 gammcf=exp(-x+a*alog(x)-gln)*g

      return
      end

*     ===================
      function gammln(xx)
*     ===================

      implicit	none

      integer		j
      real		xx,gammln
      real*8 		cof(6),stp,half,one,fpf,x,tmp,ser

      data cof,stp/76.18009173D0,-86.50532033D0,24.01409822D0,
     +     -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
      data half,one,fpf/0.5D0,1.0D0,5.5D0/

      x=xx-one
      tmp=x+fpf
      tmp=(x+half)*log(tmp)-tmp
      ser=one

      do j=1,6
         x=x+one
         ser=ser+cof(j)/x
      end do

      gammln=tmp+log(stp*ser)

      return
      end

*     ===============================
      subroutine funcs(x,a,y,dyda,na)
*     ===============================

      implicit	none

      integer		na
      real	 	x,y,a(na),dyda(na)

*     ---------------
*     fit exponentiel
*     ---------------

      y=a(1)*exp(-a(2)*x)
      dyda(1)=exp(-a(2)*x)
      dyda(2)=-a(1)*x*exp(-a(2)*x)     

      return
      end
