!**********************************************************************!
!
! File: lfit.f90 (08-Jun-2001)
!                (08-Jun-2006)
!                (22-Apr-2008)
!
! Lineares Ausgleichsproblem mit Fitfunktion der Form
!
!            _m_
!           \
!    y(x) = /___ c(k)*phi_k(x)
!            k=1
!
! auf dem Intervall [0,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 (enthaelt in erster
!                Zeile Anzahl der Datenpunkte n)
!
! nx.............Anzahl der Datenpunkte fuer Darstellung der gefitteten
!                Funktion im Intervall [0,1]
! lfit_out.dat...Ausgabefile mit gefitteter Funktion
!
!**********************************************************************!

module lfit_m

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

  implicit none

  private

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

  public::gaussp,phi

contains

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

  function phi(k,x) result(phix)

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

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

    real(kind=double)::phix

! Basisfunktionen fuer Fit
    
    if(k==0) then
      phix=x*(1.0-x)
    else
      phix=x*(1.0-x)*(2.0*x-1.0)**k
    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(size(b))::iperm
    real(kind=double)::pivot,pmax
    real(kind=double),dimension(size(b),size(b)+1)::a

    n=size(a_in,dim=1)   ! 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

      m=maxloc(abs(a(iperm(i:n),i)),dim=1)+i-1
      pmax=abs(a(iperm(m),i))
      if(pmax<eps) then                        ! Matrix singulaer?
        write(unit=*,fmt="(a)") " 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)") " |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,nx
  real(kind=double)::chi2,dx,xfit,yfit
  real(kind=double),dimension(:),allocatable::b,c,sig,x,y
  real(kind=double),dimension(:,:),allocatable::a

! Anzahl Fitparameter & Punkte Ausgabefile

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

! Daten einlesen

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

  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

  call gaussp(a,b,c)

! Werte der Fitparamter

  do k=1,m
    write(unit=*,fmt="(a,i3,a,tr1,es12.5)") " c_k, k=",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="(a,es12.5)") " chi2_red=",chi2/(n-m)

! Fit auf File schreiben

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

  dx=1.0/nx

  do i=0,nx

    xfit=i*dx
    yfit=0.0
    do k=1,m
      yfit=yfit+c(k)*phi(k-1,xfit)
    end do

    write(unit=iout,fmt="(tr1,f7.5,tr1,es12.5)") xfit,yfit

  end do

  close(unit=iout)

end program lfit

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