!**********************************************************************!
!
! File: lfit.f90 (08-Jun-2001)
!                (03-Jun-2005)
!
! Lineares Ausgleichsproblem mit Fitfunktion der Form
!
!            _m_
!           \
!    y(x) = /___ c(k)*phi_k(x)
!            k=1
!
! wo
!
!   phi_k(x) = x*(1-x)*(2*x-1)**(k-1)
!
! Loesung des linearen Gleichungssystems fuer Fitparameter mit
! Gauss'scher Elimination und partieller Pivotierung
!
! m..............Anzahl der Fitparameter c_k
! lfit_in.dat....Eingabefile mit zu fittenden Daten
! lfit_out.dat...Ausgabefile mit gefitteter Funktion im Intervall [0,1]
!
!**********************************************************************!

module lfit_m

  implicit none

  public::gaussp,phi

  integer,parameter,public::double=selected_real_kind(15)
  real(kind=double),parameter,public::eps=1.0e-08

contains

!**********************************************************************!

  function phi(k,x) result(phix)

!**********************************************************************!

    integer,intent(in)::k
    real(kind=double),intent(in)::x
    real(kind=double)::phix

    if(k==0) then
      phix=x*(1.0-x)
    else
      phix=x*(1.0-x)*(2.0*x-1.0)**k   ! Basisfunktionen fuer Fit
    end if
    
    return
    
  end function phi

!**********************************************************************!

  subroutine gaussp(a_in,b,x)

!**********************************************************************!

    real(kind=double),dimension(:,:),intent(in)::a_in
    real(kind=double),dimension(:),intent(in)::b
    real(kind=double),dimension(:),intent(out)::x
    integer::i,itemp,j,k,m,n
    integer,dimension(1)::mvect
    integer,dimension(size(b))::iperm
    real(kind=double)::pivot,pmax
    real(kind=double),dimension(size(b),size(b)+1)::a

    n=size(b)   ! Dimension des Gleichungssystems

! Kopie der Koeffizientenmatrix a_in und der rechten Seite b
! in EINE n*(n+1) Matrix a (d.h. a_in wird nicht ueberschrieben)

    a(:,1:n)=a_in
    a(:,n+1)=b

! Permutationsvektor initialisieren

    iperm=(/(j,j=1,n)/)

    do i=1,n-1

! Maximales Pivotelement in Spalte i suchen

      mvect=maxloc(abs(a(iperm(i:n),i)))
      m=mvect(1)+i-1
      pmax=abs(a(iperm(m),i))
      if(pmax<=eps) then                       ! Matrix singulaer?
        write(unit=*,fmt="("" pmax<=eps"")")
        stop
      end if

! Zeilen (nur Indizes!) vertauschen

      if(m/=i) then
        itemp=iperm(i)
        iperm(i)=iperm(m)
        iperm(m)=itemp
      end if

! Elimination

      pivot=a(iperm(i),i)
      do k=i+1,n
       a(iperm(k),i+1:n+1)=a(iperm(k),i+1:n+1)- &
            a(iperm(i),i+1:n+1)*a(iperm(k),i)/pivot
      end do

    end do

! Ruecksubstitution

    if(abs(a(iperm(n),n))<=eps) then
      write(unit=*,fmt="("" |a(n,n)|<=eps"")")
      stop
    end if

    x(n)=a(iperm(n),n+1)/a(iperm(n),n)
    do i=n-1,1,-1
      x(i)=(a(iperm(i),n+1)-sum(a(iperm(i),i+1:n)*x(i+1:n)))/a(iperm(i),i)
    end do

    return

  end subroutine gaussp

end module lfit_m

!**********************************************************************!

program lfit

!**********************************************************************!

  use lfit_m

  implicit none

  integer,parameter::iin=1,iout=2

  integer::i,j,k,m,n
  real(kind=double)::chi2,xfit,yfit
  real(kind=double),dimension(:),allocatable::b,c,sig,x,y
  real(kind=double),dimension(:,:),allocatable::a

! Anzahl Fitparameter

  write(unit=*,fmt="("" m="")",advance="no")
  read(unit=*,fmt=*) m

! Daten einlesen

  open(unit=iin,file="lfit_in.dat",status="old",action="read", &
       form="formatted")

  read(unit=iin,fmt=*) n

  allocate(a(m,m),b(m),c(m),sig(n),x(n),y(n))

  do i=1,n
    read(unit=iin,fmt=*) x(i),y(i),sig(i)
  end do

  close(unit=iin)

! Koeffizienten des linearen Gleichungssystems berechnen

  a=0.0
  b=0.0
  do i=1,m
    do k=1,n
      b(i)=b(i)+phi(i-1,x(k))*y(k)/sig(k)**2
    end do
    do j=1,m
      do k=1,n
        a(i,j)=a(i,j)+phi(i-1,x(k))*phi(j-1,x(k))/sig(k)**2
      end do
    end do
  end do

! Gleichungssystem loesen

  c
all gaussp(a,b,c)

! Werte der Fitparamter

  do k=1,m
    write(unit=*,fmt="("" k,c_k: "",i3,"","",tr1,es12.5)") k,c(k)
  end do

! Reduziertes chi**2

  chi2=0.0
  do i=1,n
    yfit=0.0
    do k=1,m
      yfit=yfit+c(k)*phi(k-1,x(i))
    end do
    chi2=chi2+((yfit-y(i))/sig(i))**2
  end do
  write(unit=*,fmt="("" chi2_red="",es12.5)") chi2/(n-m)

! Fit auf File schreiben

  open(unit=iout,file="lfit_out.dat",status="replace",action="write", &
       form="formatted")

  do i=0,100
    xfit=0.01*i
    yfit=0.0
    do k=1,m
      yfit=yfit+c(k)*phi(k-1,xfit)
    end do
    write(unit=iout,fmt="(tr1,f5.3,tr1,es12.5)") xfit,yfit
  end do

  close(unit=iout)

end program lfit

!**********************************************************************!