!>
!> @file l2err.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine l2err ( iprfun, ftau, error )

!*************************************************************************
!
!! L2ERR computes the errors of an L2 approximation.
!
!  Discussion:
!
!    This routine computes various errors of the current L2-approximation,
!    whose piecewise polynomial representation is contained in common
!    block APPROX, to the given data contained in common block  data.
!
!    It prints out the average error errl1, the l2-error errl2, and the
!    maximum error errmax.
!
!  Reference:
!
!    Carl DeBoor,
!    A Practical Guide to Splines,
!    Springer Verlag.
!
!  Parameters:
!
!    Input, integer IPRFUN.  If iprfun= 1, the routine prints out
!    the value of the approximation as well as its error at
!    every data point.
!
!    Output, real ( kind = 8 ) FTAU(NTAU), contains the value of the computed
!    approximation at each value TAU(1:NTAU).
!
!    Output, error(1), ..., error(ntau),  with  error(i)=scale*(g-f)
!    at tau(i), all i. here, SCALE equals 1. in case
!    iprfun /= 1 , or the absolute error is greater than 100 some-
!    where. otherwise, SCALE is such that the maximum of
!    abs(error))  over all I lies between 10 and 100.  This
!    makes the printed output more illustrative.
!
  implicit none

  integer, parameter :: lpkmax = 100
  integer, parameter :: ntmax = 200
  integer, parameter :: ltkmax = 2000

  integer ntau

  real ( kind = 8 ) break
  real ( kind = 8 ) coef
  real ( kind = 8 ) err
  real ( kind = 8 ) errl1
  real ( kind = 8 ) errl2
  real ( kind = 8 ) errmax
  real ( kind = 8 ) error(ntau)
  real ( kind = 8 ) ftau(ntau)
  real ( kind = 8 ) gtau
  integer ie
  integer iprfun
  integer k
  integer l
  integer ll
  real ( kind = 8 ) scale
  real ( kind = 8 ) tau
  real ( kind = 8 ) totalw
  real ( kind = 8 ) weight

  COMMON /DATA/ tau(ntmax),gtau(ntmax),weight(ntmax),totalw,ntau
  common /approx/ break(lpkmax),coef(ltkmax),l,k

  errl1 = 0.0D+00
  errl2 = 0.0D+00
  errmax = 0.0D+00

  do ll = 1, ntau
    call ppvalu(break,coef,l,k,tau(ll),0,ftau(ll))
    error(ll) = gtau(ll)-ftau(ll)
    err = abs(error(ll))
    if ( errmax < err ) then
      errmax = err
    end if
    errl1 = errl1 + err * weight(ll)
    errl2 = errl2 + err**2 * weight(ll)
  end do

  errl1 = errl1 / totalw
  errl2 = sqrt ( errl2 / totalw )

  write ( *, * ) ' '
  write ( *, * ) ' Least square error =',errl2
  write ( *, * ) ' Average error     =',errl1
  write ( *, * ) ' Maximum error     =',errmax
  write ( *, * ) ' '

  if ( iprfun /= 1 ) then
    return
  end if
!
!  Scale error curve and print
!
  ie = 0
  scale = 1.0D+00

  if ( errmax < 10.0D+00 ) then

    do ie = 1, 9
      scale = scale * 10.0D+00
      if ( 10.0D+00 <= errmax * scale ) then
        exit
      end if
    end do

  end if

  error(1:ntau) = error(1:ntau) * scale

  write(*,60) ie, (ll,tau(ll),ftau(ll),error(ll),ll=1,ntau)

   60 format (///14x,'approximation and scaled error curve'/ &
       7x,'data point',7x,'approximation',3x,'deviation x 10**',i1/ &
       (i4, f16.8,f16.8,f17.6))

  return
end