%\documentclass[a4paper,10pt]{article}
\documentclass[a4paper,english]{article}

\usepackage[utf8]{inputenc}
\usepackage{a4wide}
\usepackage{eurosym}
\usepackage{url}
\usepackage[colorlinks=true,hyperfigures=true]{hyperref}
%\usepackage{aeguill}

\usepackage{graphicx}
\usepackage{babel}
\makeatother


%opening
\title{FVN Documentation}
\author{William Daniau}


\begin{document}

\maketitle

%\begin{abstract}

%\end{abstract}
\tableofcontents

\section{Whatis fvn,license}
\subsection{Whatis fvn}
fvn is a Fortran95 mathematical library with several modules. It provides various usefull subroutine covering linear algebra, numerical integration, least square polynomial, spline interpolation, zero finding, special functions etc.

Most of the work for linear algebra is done by interfacing Lapack \url{http://www.netlib.org/lapack} which means that Lapack and Blas \url{http://www.netlib.org/blas} must be available on your system for linking fvn. If you use an AMD microprocessor, the good idea is to use ACML ( AMD Core Math Library \url{http://developer.amd.com/acml.jsp} which contains an optimized Blas/Lapack.

fvn include some integrated libraries : integration tasks uses a slightly modified version of Quadpack \url{http://www.netlib.org/quadpack}, the fnlib library \url{http://www.netlib.org/fn} is used for special functions and sparse system resolution uses SuiteSparse \url{http://www.cise.ufl.edu/research/sparse/SuiteSparse/}.

This library has been initially written for the use of the ``Acoustic and microsonic'' group leaded by Sylvain Ballandras in the Time and Frequency Department of institute Femto-ST \url{http://www.femto-st.fr/}.

\subsection{License}
    Your use or distribution of fvn or any modified version of
    fvn implies that you agree to this License.

    This library 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.

    This library 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
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
    USA

    Permission is hereby granted to use or copy this program under the
    terms of the GNU LGPL, provided that the Copyright, this License,
    and the Availability of the original version is retained on all copies.
    User documentation of any code that uses this code or any modified
    version of this code must cite the Copyright, this License, the
    Availability note, and "Used by permission." Permission to modify
    the code and to distribute modified code is granted, provided the
    Copyright, this License, and the Availability note are retained,
    and a notice that the code was modified is included.

\subsubsection*{Authors and Contributors}
\begin{itemize}
 \item William Daniau (william.daniau@femto-st.fr)
 \item Sylvain Ballandras (sylvain.ballandras@femto-st.fr)
 \item Christian Waterkeyn (christian.waterkeyn@epcos.com)
\end{itemize}

\section{Naming scheme and convention}
The naming scheme of the routines is as follow :
\begin{verbatim}
    fvn_*_name()
\end{verbatim} 
where * can be s,d,c or z. 
\begin{itemize}
 \item s is for single precision real (real,real*4,real(4),real(kind=4))
 \item d for double precision real (double precision,real*8,real(8),real(kind=8))
 \item c for single precision complex (complex,complex*8,complex(4),complex(kind=4))
 \item z for double precision complex (double complex,complex*16,complex(8),complex(kind=8))
\end{itemize}
In the following description of subroutines parameters, input parameters are followed by (in), output parameters by (out) and parameters which are used as input and modified by the subroutine are followed by (inout).

For each routine, there is a generic interface (simply remove \verb'_*' in the name), so using the specific routine is not mandatory.


There's a general module called ``\verb'fvn''' that include all fvn submodules. So whatever part of the library is used in the program a ``\verb'use fvn''' will be sufficient instead of specifying the specific module.

\section{kind specification}
In \verb'fvn' the following definitions are made (this is done in module \verb'fvn_common') to ease portability:

\begin{verbatim}
integer, parameter :: ip_kind = kind(1)
integer, parameter :: sp_kind = kind(1.0E0)
integer, parameter :: dp_kind = kind(1.0D0)
\end{verbatim}

Although not mandatory, it is a good idea to use these definitions when programming with fvn, that is for example : 
\begin{verbatim}
real(kind=sp_kind) :: x 
real(kind=dp_kind) :: y
complex(kind=dp_kind) :: z
\end{verbatim}
instead of
\begin{verbatim}
real :: x
double precision :: y
complex(8) :: z
\end{verbatim}


\section{Linear algebra}
The linear algebra routines of fvn are an interface to lapack, which make it easier to use.
\subsection{Matrix inversion}
\begin{verbatim}
Module : use fvn_linear
call fvn_matinv(d,a,inva,status)
\end{verbatim}
\begin{itemize}
 \item d (in) is an integer equal to the matrix rank
 \item a (in) is a real or complex matrix. It will remain untouched.
 \item inva (out) is a real or complex matrix which contain the inverse of a at the end of the routine
 \item status (out) is an optional integer equal to zero if something went wrong
\end{itemize}

\subsubsection*{Example}
\begin{verbatim}
program inv
 use fvn_linear
 implicit none

 real(8),dimension(3,3) :: a,inva

 call random_number(a)
 a=a*100

 call fvn_matinv(3,a,inva)
 write (*,*) a
 write (*,*)
 write (*,*) inva
 write (*,*)
 write (*,*) matmul(a,inva)
end program
\end{verbatim}



\subsection{Matrix determinants}
\begin{verbatim}
Module : use fvn_linear
det=fvn_det(d,a,status)
\end{verbatim}
\begin{itemize}
 \item d (in) is an integer equal to the matrix rank
 \item a (in) is a real or complex matrix. It will remain untouched.
 \item status (out) is an optional integer equal to zero if something went wrong
\end{itemize}

\subsubsection*{Example}
\begin{verbatim}
program det
 use fvn_linear
 implicit none

 real(8),dimension(3,3) :: a
 real(8) :: deta
 integer :: status

 call random_number(a)
 a=a*100

 deta=fvn_det(3,a,status)
 write (*,*) a
 write (*,*)
 write (*,*) "Det = ",deta
end program

\end{verbatim}



\subsection{Matrix condition}
\begin{verbatim}
Module : use fvn_linear
call fvn_matcon(d,a,rcond,status)
\end{verbatim}
\begin{itemize}
 \item d (in) is an integer equal to the matrix rank
 \item a (in) is a real or complex matrix. It will remain untouched.
 \item rcond (out) is a real of same kind as matrix a, it will contain the reciprocal condition number of the matrix
 \item status (out) is an optional integer equal to zero if something went wrong
\end{itemize}

The reciprocal condition number is evaluated using the 1-norm and is define as in equation \ref{rconddef}
\begin{equation}
 R = \frac{1}{norm(A)*norm(invA)}
 \label{rconddef}
\end{equation}

The 1-norm itself is defined as the maximum value of the columns absolute values (modulus for complex) sum as in equation \ref{l1norm}
\begin{equation}
 L1 = max_j ( \sum_i{\mid A(i,j)\mid}  )
 \label{l1norm}
\end{equation}

\subsubsection*{Example}
\begin{verbatim}
program cond
 use fvn_linear
 implicit none

 real(8),dimension(3,3) :: a
 real(8) :: rcond
 integer :: status

 call random_number(a)
 a=a*100

 call fvn_d_matcon(3,a,rcond,status)
 write (*,*) a
 write (*,*)
 write (*,*) "Cond = ",rcond
end program

\end{verbatim}


\subsection{Eigenvalues/Eigenvectors}
\begin{verbatim}
Module : use fvn_linear
call fvn_matev(d,a,evala,eveca,status,sortval)
\end{verbatim}
\begin{itemize}
 \item d (in) is an integer equal to the matrix rank
 \item a (in) is a real or complex matrix. It will remain untouched.
 \item evala (out) is a complex array of same kind as a. It contains the eigenvalues of matrix a
 \item eveca (out) is a complex matrix of same kind as a. Its columns are the eigenvectors of matrix a : eveca(:,j)=jth eigenvector associated with eigenvalue evala(j).
 \item status (out) is an optional integer equal to zero if something went wrong
 \item sortval (in) is an optional logical, if it is true the eigenvalues (and eigenvectors) are sorted in decreasing order of eigenvalues's modulus.
\end{itemize}


\subsubsection*{Example}
\begin{verbatim}
program eigen
 use fvn_linear
 implicit none

 real(8),dimension(3,3) :: a
 complex(8),dimension(3) :: evala
 complex(8),dimension(3,3) :: eveca
 integer :: status,i,j

 call random_number(a)
 a=a*100

 call fvn_matev(3,a,evala,eveca,status)
 write (*,*) a
 write (*,*)
 do i=1,3
    write(*,*) "Eigenvalue ",i,evala(i)
    write(*,*) "Associated Eigenvector :"
    do j=1,3
        write(*,*) eveca(j,i)
    end do
    write(*,*)
 end do

end program

\end{verbatim}


\subsection{Sparse matrix}
By interfacing Tim Davis's SuiteSparse from university of Florida \url{http://www.cise.ufl.edu/research/sparse/SuiteSparse/} which is a reference for this kind of problems, fvn provides simple subroutines for solving linear sparse systems and calculating determinants.

\subsubsection{Sparse solving}
The provided routines solves the equation $Ax=B$ where A is sparse and given in its triplet form.

Note that interface has changed from previous versions, splitting complex inputs into real and imaginary part and using 0-based indices, to fit umfpack routines inputs and avoid duplicating datas inside the routines.

\begin{verbatim}
Module : fvn_sparse
call fvn_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det)
call fvn_sparse_solve(n,nz,Tx,Tz,Ti,Tj,Bx,Bz,x,status,det)
\end{verbatim}
\begin{itemize}
 \item For this family of subroutines the two letters (zl,zi,dl,di) of the specific interface name decribe the arguments's type. z is for complex(8), d for real(8), l for integer(8) and i for integer(4)
 \item \texttt{n} (in) is an integer equal to the matrix rank
 \item \texttt{nz} (in) is an integer equal to the number of non-zero elements
 \item \texttt{T(nz) or Tx(nz),Tz(nz)} (in) is a real array (or two real arrays for real/imaginary in case of a complex system) containing the non-zero elements
 \item \texttt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, this has to be 0-based as in C.
 \item \texttt{B(n) or Bx(n),Bz(n)} (in) is a real array or two real arrays for real/imaginary in case of a complex system) containing the second member of the equation.
 \item \texttt{x(n)} (out) is a complex/real array containing the solution
 \item \texttt{status} (out) is an integer which contain non-zero is something went wrong
 \item \texttt{det} (out), is an optional real(8) array of dimension 2 for dl and di specific interface (real systems) and dimension 3 for zl and zi interface (complex systems)
\end{itemize}


\subsubsection*{Example}
\begin{verbatim}
program test_sparse

 use fvn_sparse
 implicit none

 integer(8), parameter :: nz=12
 integer(8), parameter :: n=5
 complex(8),dimension(nz) :: A
 real(8), dimension(nz) :: Ax,Az
 integer(8),dimension(nz) :: Ti,Tj
 complex(8),dimension(n) :: B,x
 real(8), dimension(n) :: Bx,Bz
 integer(8) :: status

 A = (/ (2.,0.),(3.,0.),(3.,0.),(-1.,0.),(4.,0.),(4.,0.),(-3.,0.),&
          (1.,0.),(2.,0.),(2.,0.),(6.,0.),(1.,0.) /)
 B = (/ (8.,0.), (45.,0.), (-3.,0.), (3.,0.), (19.,0.)/)
 Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
 Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)
 Ax=real(A)
 Az=aimag(A)
 Bx=real(B)
 Bz=aimag(B)

 ! 1-based to 0-based translation
 Ti=Ti-1
 Tj=Tj-1

 !specific routine that will be used here
 !call fvn_zl_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status)
 call fvn_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status)
 write(*,*) x

end program


program test_sparse

use fvn_sparse
implicit none

integer(4), parameter :: nz=12
integer(4), parameter :: n=5
real(8),dimension(nz) :: A
integer(4),dimension(nz) :: Ti,Tj
real(8),dimension(n) :: B,x
integer(4) :: status

A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
B = (/ 8., 45., -3., 3., 19./)
Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)

! 1-based to 0-based translation
Ti=Ti-1
Tj=Tj-1
 
!specific routine that will be used here
!call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
write(*,*) x

end program

\end{verbatim}

If optional parameter \texttt{det} is given, the routine will also calculates the matrix determinant and returns it on a mantissa + exponent form, that is the actual determinant will be $det(1).10^{det(2)}$ for real problems and $(det(1)+i.det(2)).10^{det(3)}$ for complex problems. This is given in this form as the determinant can be considerably higher/lower than the biggest/lowest usable double precision real. There's an example of how to use this in following paragraph.

\subsubsection{Sparse determinant}

The provided subroutines calculates the determinant of a matrix given in its triplet form.
\begin{verbatim}
Module : fvn_sparse
call fvn_sparse_det(n,nz,T,Ti,Tj,det,status)
call fvn_sparse_det(n,nz,Tx,Tz,Ti,Tj,det,status)
\end{verbatim}

\begin{itemize}
 \item For this family of subroutines the two letters (zl,zi,dl,di) of the specific interface name decribe the arguments's type. z is for complex(8), d for real(8), l for integer(8) and i for integer(4)
 \item \texttt{n} (in) is an integer equal to the matrix rank
 \item \texttt{nz} (in) is an integer equal to the number of non-zero elements
 \item \texttt{T(nz) or Tx(nz),Tz(nz)} (in) is a real array (or two real arrays for real/imaginary in case of a complex system) containing the non-zero elements
 \item \texttt{Ti(nz)},\texttt{Tj(nz)} (in) are the indexes of the corresponding element of \texttt{T} in the original matrix, it has to be 0-based as in C.
 \item \texttt{det} (out), a real(8) array of dimension 2 for dl and di specific interface (real systems) and dimension 3 for zl and zi interface (complex systems)
\item \texttt{status} (out) is an integer which contain non-zero is something went wrong
\end{itemize}

The matrix determinant is returned on a mantissa + exponent form, that is the actual determinant will be $det(1).10^{det(2)}$ for real problems and $(det(1)+i.det(2)).10^{det(3)}$ for complex problems. This is given in this form as the determinant can be considerably higher/lower than the biggest/lowest usable double precision real.

Here are the possibly returned errors in \texttt{status} parameter :
\begin{itemize}
 \item \texttt{0} : no errors
 \item \texttt{-1}: out of memory
 \item \texttt{1} : singular matrix
 \item \texttt{2} : determinant underflow, the ``natural'' form of the determinant $det(1).10^{det(2)}$ or $(det(1)+i.det(2)).10^{det(3)}$ will underflow.
 \item \texttt{3} : determinant overflow, the ``natural'' form of the determinant (as above) will overflow
\end{itemize}


And here's an example using this

\begin{verbatim}
program test_sparse
use fvn
implicit none
integer(kind=sp_kind), parameter :: nz=12
integer(kind=sp_kind), parameter :: n=5
complex(kind=dp_kind),dimension(nz) :: A
real(kind=dp_kind), dimension(nz) :: Ax,Az
complex(kind=dp_kind),dimension(n,n) :: As
integer(kind=sp_kind),dimension(nz) :: Ti,Tj
complex(kind=dp_kind),dimension(n) :: B,x
real(kind=dp_kind), dimension(n) :: Bx,Bz
integer(kind=sp_kind) :: status,i
real(kind=dp_kind),dimension(3) :: det
character(len=80) :: fmcmplx

fmcmplx='(5("(",f8.5,",",f8.5,")  "))'

! Description of the matrix in triplet form
A = (/ (2.,-1.),(3.,2.),(3.,1.),(-1.,5.),(4.,-7.),(4.,0.),(-3.,-4.),(1.,3.),(2.,0.),(2.,-2.),(6.,4.),(1.,0.) /)
B = (/ (8.,3.), (45.,1.), (-3.,-2.), (3.,0.), (19.,2.) /)
Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)

Ax=real(A)
Az=aimag(A)
Bx=real(B)
Bz=aimag(B)

! Reconstruction of the matrix in standard form
As=0.
do i=1,nz
    As(Ti(i),Tj(i))=A(i)
end do

write(*,*) "Matrix in standard representation :"
do i=1,5
    write(*,fmcmplx) As(i,:)
end do
write(*,*)
write(*,*) "Standard determinant : ",fvn_det(5,As)
write(*,*)
write(*,*) "Right hand side :"
write(*,fmcmplx) B

! can use either specific interface, fvn_zi_sparse_det
! either generic one fvn_sparse_det
call fvn_zi_sparse_det(n,nz,Ax,Az,Ti,Tj,det,status)
write(*,*)
write(*,*) "Sparse Det = ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
! can use either specific interface fvn_zi_sparse_solve
! either generic one fvn_sparse_solve
! parameter det is optional
call fvn_zi_sparse_solve(n,nz,Ax,Az,Ti,Tj,Bx,Bz,x,status,det)
write(*,*)
write(*,*) "Sparse Det as solve option= ",cmplx(det(1),det(2),kind=dp_kind)*10**det(3)
write(*,*)
write(*,*) "Solution :"
write(*,fmcmplx) x
write(*,*)
write(*,*) "Product matrix Solution :"
write(*,fmcmplx) matmul(As,x)
end program

\end{verbatim}


\subsection{Identity matrix}
\begin{verbatim}
 Module : use fvn_linear
 I=fvn_*_ident(n)     (*=s,d,c,z)
\end{verbatim}
\begin{itemize}
 \item n (in) is an integer equal to the matrix rank
\end{itemize}

This function return the identity matrix of rank n, in the specified type. No generic interface for this one.

\subsection{Operators}
fvn defines some linear operators similar to those defined in IMSL\textregistered (\url{http://www.vni.com/products/imsl/}), that can be used for matrix operations.
\begin{verbatim}
 Module : use fvn_linear
\end{verbatim}

\subsubsection{Unary operators}
\paragraph{.i.}
This operator gives the inverse matrix of the argument which must be a square matrix. The status of the operation can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).
\[
 \textrm{b=.i.a}~\Longleftrightarrow~b=a^{-1}
\]

\paragraph{.t.}
This operator gives the transpose matrix of the argument.
\[
 \textrm{b=.t.a}~\Longleftrightarrow~b= {}^t \! a
\]


\paragraph{.h.}
This operator gives the conjugate transpose matrix of the argument (also called Hermitian transpose or adjoint matrix).
\[
 \textrm{b=.h.a}~\Longleftrightarrow~b=a^*=\overline{({}^t\!a)}={}^t \overline{a}
\]


\subsubsection{Binary operators}
\paragraph{.x.}
This operator gives the matrix product of the two operands.
\[
 \textrm{c=a.x.b}~\Longleftrightarrow~c=ab
\]


\paragraph{.ix.}
This operator gives the matrix product of the inverse of the first operand and the second one. The status of the inversion can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).
\[
 \textrm{c=a.ix.b}~\Longleftrightarrow~c=a^{-1}b
\]

\paragraph{.xi.}
This operator gives the matrix product of the first operand and the inverse of the second one. The status of the inversion can be found in the module variable fvn\_status (fourth parameter fvn\_matinv).
\[
 \textrm{c=a.xi.b}~\Longleftrightarrow~c=ab^{-1}
\]

\paragraph{.tx.}
This operator gives the matrix product of the transpose matrix of the first operand and the second one.
\[
 \textrm{c=a.tx.b}~\Longleftrightarrow~c={}^t\!ab
\]


\paragraph{.xt.}
This operator gives the matrix product of the first operand and the transpose matrix of the second one.
\[
 \textrm{c=a.xt.b}~\Longleftrightarrow~c=a{}^t\!b
\]

\paragraph{.hx.}
This operator gives the matrix product of the conjugate transpose matrix of the first operand and the second one.
\[
 \textrm{c=a.hx.b}~\Longleftrightarrow~c=a^*b
\]

\paragraph{.xh.}
This operator gives the matrix product of thee first operand and the conjugate transpose matrix of the second one.
\[
 \textrm{c=a.xh.b}~\Longleftrightarrow~c=ab^*
\]


\section{Interpolation}

\subsection{Quadratic Interpolation}
fvn provide function for interpolating values of a tabulated function of 1, 2 or 3 variables, for both single and double precision.
\subsubsection{One variable function}
\begin{verbatim}
 Module : use fvn_interpol
 value=fvn_quad_interpol(x,n,xdata,ydata)
\end{verbatim}
\begin{itemize}
 \item x is the real where we want to evaluate the function
 \item n is the number of tabulated values
 \item xdata(n) contains the tabulated coordinates
 \item ydata(n) contains the tabulated function values ydata(i)=y(xdata(i))
\end{itemize}
xdata must be strictly increasingly ordered.
x must be within the range of xdata to actually perform an interpolation, otherwise the resulting value is an extrapolation
\paragraph*{Example}
\begin{verbatim}
program inter1d

use fvn_interpol
implicit none

integer(kind=4),parameter :: ndata=33
integer(kind=4) :: i,nout
real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)
real(kind=8) ::tv

intrinsic sin

f(x)=sin(x)

xdata(1)=0.
fdata(1)=f(xdata(1))
h=1./32.
do i=2,ndata
      xdata(i)=xdata(i-1)+h
      fdata(i)=f(xdata(i))
end do
call random_seed()
call random_number(x)

q=fvn_d_quad_interpol(x,ndata,xdata,fdata)

tv=f(x)
write(*,*) "x ",x
write(*,*) "Calculated (real) value :",tv
write(*,*) "fvn interpolation :",q
write(*,*) "Relative fvn error :",abs((q-tv)/tv)

end program

\end{verbatim}


\subsubsection{Two variables function}
\begin{verbatim}
Module : use fvn_interpol
value=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
\end{verbatim}
\begin{itemize}
 \item x,y are the real coordinates where we want to evaluate the function
 \item nx is the number of tabulated values along x axis
 \item xdata(nx) contains the tabulated x
 \item ny is the number of tabulated values along y axis
 \item ydata(ny) contains the tabulated y
 \item zdata(nx,ny) contains the tabulated function values zdata(i,j)=z(xdata(i),ydata(j))
\end{itemize}
xdata and ydata must be strictly increasingly ordered.
(x,y) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation

\paragraph*{Example}

\begin{verbatim}
program inter2d
use fvn_interpol
implicit none

integer(kind=4),parameter  :: nx=21,ny=42
integer(kind=4) :: i,j
real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
real(kind=8) :: tv

intrinsic dble,sin

f(x,y)=sin(x+2.*y)
do i=1,nx
      xdata(i)=dble(i-1)/dble(nx-1)
end do
do i=1,ny
      ydata(i)=dble(i-1)/dble(ny-1)
end do
do i=1,nx
      do j=1,ny
            fdata(i,j)=f(xdata(i),ydata(j))
      end do
end do
call random_seed()
call random_number(x)
call random_number(y)

q=fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
tv=f(x,y)

write(*,*) "x y",x,y
write(*,*) "Calculated (real) value :",tv
write(*,*) "fvn interpolation :",q
write(*,*) "Relative fvn error :",abs((q-tv)/tv)

end program

\end{verbatim}



\subsubsection{Three variables function}
\begin{verbatim}
 Module : use fvn_interpol
value=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
\end{verbatim}
\begin{itemize}
 \item x,y,z are the real coordinates where we want to evaluate the function
 \item nx is the number of tabulated values along x axis
 \item xdata(nx) contains the tabulated x
 \item ny is the number of tabulated values along y axis
 \item ydata(ny) contains the tabulated y
 \item nz is the number of tabulated values along z axis
 \item zdata(ny) contains the tabulated z
 \item tdata(nx,ny,nz) contains the tabulated function values tdata(i,j,k)=t(xdata(i),ydata(j),zdata(k))
\end{itemize}
xdata, ydata and zdata must be strictly increasingly ordered.
(x,y,z) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation

\paragraph*{Example}
\begin{verbatim}
program inter3d
use fvn_interpol

implicit none

integer(kind=4),parameter  :: nx=21,ny=42,nz=18
integer(kind=4) :: i,j,k
real(kind=8) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz)
real(kind=8) :: tv

intrinsic dble,sin

f(x,y,z)=sin(x+2.*y+3.*z)
do i=1,nx
      xdata(i)=2.*(dble(i-1)/dble(nx-1))
end do
do i=1,ny
      ydata(i)=2.*(dble(i-1)/dble(ny-1))
end do
do i=1,nz
      zdata(i)=2.*(dble(i-1)/dble(nz-1))
end do
do i=1,nx
      do j=1,ny
            do k=1,nz
                  fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
            end do
      end do
end do
call random_seed()
call random_number(x)
call random_number(y)
call random_number(z)

q=fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
tv=f(x,y,z)

write(*,*) "x y z",x,y,z
write(*,*) "Calculated (real) value :",tv
write(*,*) "fvn interpolation :",q
write(*,*) "Relative fvn error :",abs((q-tv)/tv)

end program

\end{verbatim}

\subsubsection{Utility procedure}
fvn provides a simple utility procedure to locate the interval in which a value is located in an increasingly ordered array.
\begin{verbatim}
Module : use fvn_interpol
call fvn_find_interval(x,i,xdata,n)
\end{verbatim}
\begin{itemize}
 \item x (in) the real value to locate
 \item i (out) the resulting indice
 \item xdata(n) (in) increasingly ordered array
 \item n (in) size of the array
\end{itemize}
The resulting integer i is as : $xdata(i) <= x < xdata(i+1)$. If $x < xdata(1)$ then $i=0$ is returned. If $x > xdata(n)$ then $i=n$ is returned. Finally if $x=xdata(n)$ then $i=n-1$ is returned.



\subsection{Akima spline}
fvn provides Akima spline interpolation and evaluation for both single and double precision real.
\subsubsection{Interpolation}
\begin{verbatim}
Module : use fvn_interpol
call fvn_akima(n,x,y,br,co)
\end{verbatim}
\begin{itemize}
 \item n (in) is an integer equal to the number of points
 \item x(n) (in) ,y(n) (in) are the known couples of coordinates
 \item br (out) on output contains a copy of x
 \item co(4,n) (out) is a real matrix containing the 4 coefficients of the Akima interpolation spline for a given interval.
\end{itemize}

\subsubsection{Evaluation}
\begin{verbatim}
Module : use fvn_interpol
y=fvn_spline_eval(x,n,br,co)
\end{verbatim}
\begin{itemize}
 \item x (in) is the point where we want to evaluate
 \item n (in) is the number of known points and br(n) (in), co(4,n) (in) \\
are the outputs of fvn\_x\_akima(n,x,y,br,co) 
\end{itemize}

\subsubsection{Example}
In the following example we will use Akima splines to interpolate a sinus function with 30 points between -10 and 10. We then use the evaluation function to calculate the coordinates of 1000 points between -11 and 11, and write a 3 columns file containing : x, calculated sin(x), interpolation evaluation of sin(x).

One can see that the interpolation is very efficient even with only 30 points. Of course as soon as we leave the -10 to 10 interval, the values are extrapolated and thus can lead to very inacurrate values.

\begin{verbatim}
program akima
 use fvn_interpol
 implicit none

 integer :: nbpoints,nppoints,i
 real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
 real(8),dimension(:,:),allocatable :: coeff_fvn_d
 real(8) :: xstep_d,xp_d,ty_d,fvn_y_d

 open(2,file='fvn_akima_double.dat')
 open(3,file='fvn_akima_breakpoints_double.dat')
 nbpoints=30
 allocate(x_d(nbpoints))
 allocate(y_d(nbpoints))
 allocate(breakpoints_d(nbpoints))
 allocate(coeff_fvn_d(4,nbpoints))

 xstep_d=20./dfloat(nbpoints)
 do i=1,nbpoints
    x_d(i)=-10.+dfloat(i)*xstep_d
    y_d(i)=dsin(x_d(i))
    write(3,44) (x_d(i),y_d(i))
 end do
 close(3)

 call fvn_d_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)

 nppoints=1000 
 xstep_d=22./dfloat(nppoints)
 do i=1,nppoints
    xp_d=-11.+dfloat(i)*xstep_d
    ty_d=dsin(xp_d)
    fvn_y_d=fvn_d_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)
    write(2,44) (xp_d,ty_d,fvn_y_d)
 end do

 close(2)

44      FORMAT(4(1X,1PE22.14))

end program

\end{verbatim}
Results are plotted on figure \ref{akima}

\begin{figure}
 \begin{center}
 \includegraphics[width=0.9\textwidth]{akima.pdf}
 % akima.pdf: 504x720 pixel, 72dpi, 17.78x25.40 cm, bb=0 0 504 720
 \caption{Akima Spline Interpolation}
 \label{akima}
\end{center}

\end{figure}



\section{Least square polynomial}
fvn provide a function to find a least square polynomial of a given degree, for real in single or double precision. It is performed using Lapack subroutine sgels (dgels), which solve this problem.

\begin{verbatim}
Module : use fvn_linear
call fvn_lspoly(np,x,y,deg,coeff,status)
\end{verbatim}
\begin{itemize}
 \item np (in) is an integer equal to the number of points
 \item x(np) (in),y(np) (in) are the known coordinates
 \item deg (in) is an integer equal to the degree of the desired polynomial, it must be lower than np.
 \item coeff(deg+1) (out) on output contains the polynomial coefficients
 \item status (out) is an integer containing 0 if a problem occured.
\end{itemize}

\subsection*{Example}
Here's a simple example : we've got 13 measurement points and we want to find the least square degree 3 polynomial for these points :
\begin{verbatim}
 program lsp
 use fvn_linear
 implicit none

 integer,parameter :: npoints=13,deg=3
 integer :: status,i
 real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
 real(kind=8) :: coeff(deg+1)

 xm = (/ -3.8,-2.7,-2.2,-1.9,-1.1,-0.7,0.5,1.7,2.,2.8,3.2,3.8,4. /)
 ym = (/ -3.1,-2.,-0.9,0.8,1.8,0.4,2.1,1.8,3.2,2.8,3.9,5.2,7.5  /)

 open(2,file='fvn_lsp_double_mesure.dat')
 open(3,file='fvn_lsp_double_poly.dat')

 do i=1,npoints
    write(2,44) xm(i),ym(i)
 end do
 close(2)


 call fvn_d_lspoly(npoints,xm,ym,deg,coeff,status)

 xstep=(xm(npoints)-xm(1))/1000.
 do i=1,1000
    xc=xm(1)+(i-1)*xstep
    yc=poly(xc,coeff)
    write(3,44) xc,yc
 end do
 close(3)

44      FORMAT(4(1X,1PE22.14))

contains
function poly(x,coeff)
    implicit none
    real(8) :: x
    real(8) :: coeff(deg+1)
    real(8) :: poly
    integer :: i

    poly=0.

    do i=1,deg+1
        poly=poly+coeff(i)*x**(i-1)
    end do

end function
end program
\end{verbatim}
The results are plotted on figure \ref{lsp} .

\begin{figure}[ht]
 \begin{center}
 \includegraphics[width=0.9\textwidth]{lsp.pdf}
 \caption{Least Square Polynomial}
 \label{lsp}
 \end{center}
\end{figure}

\section{General least square fitting}
fvn provides a routine performing a general least square fitting using Levenberg-Marquardt algorithm. It uses routine \verb'lmdif' from \verb'MINPACK' (\url{http://www.netlib.org/minpack}). The used version has been converted to fortran90 by Alan Miller \verb'amiller @ bigpond.net.au'.
\newline
The purpose of \verb'fvn_lm' is to minimize the sum of the squares of m nonlinear
functions in n variables by a modification of the Levenberg-Marquardt
algorithm.  This is done by using the more general least-squares
solver lmdif.  The user must provide a subroutine which calculates the
functions.  The jacobian is then calculated by a forward-difference
approximation.

\begin{verbatim}
Module : use fvn_misc
call fvn_lm(fcn,m,n,a,info,tol)
\end{verbatim}
\begin{itemize}
\item fcn is the user-supplied subroutine which calculates
        the functions.  fcn must follow the following interface that must
        be declared in the calling subroutine :
\begin{verbatim}
         interface
             subroutine fcn(m,n,a,fvec,iflag)
                 use fvn_common
                 integer(ip_kind), intent(in) :: m
                 integer(ip_kind), intent(in) :: n
                 real(dp_kind), dimension(:), intent(in) :: a
                 real(dp_kind), dimension(:), intent(inout) :: fvec
                 integer(ip_kind), intent(inout) :: iflag
             end subroutine
         end interface
 
\end{verbatim}

          This is the function which calculate the differences for which which square sum 
          will be minimized outputing this difference in vector fvec.
           Parameters of fcn are as follows :
\begin{itemize}
 \item m : positive integer input variable set to the number of functions
               (number of measurement points)
 \item n : positive integer input variable set to the number of variables
               (number of parameters in the function to fit)
 \item a : vector of length n containing parameters for which fcn should
               perform the calculation
 \item fvec : vector of length m containing the resulting evaluation
 \item iflag : integer normally not used, can be used to exit the
               the algorithm by setting it to a negative value
\end{itemize}

\item m : positive integer input variable set to the number of functions
               (number of measurement points)
\item n : positive integer input variable set to the number of variables
               (number of parameters in the function to fit)
\item a(n) : vector of length n, on input must contains an initial guess (or unity vector)
           and on output the solution vector
\item info : is an output positive integer
\begin{itemize}
 \item info = 0  improper input parameters.
 \item info = 1  algorithm estimates that the relative error
                 in the sum of squares is at most tol.
 \item info = 2  algorithm estimates that the relative error
                between x and the solution is at most tol.
 \item info = 3  conditions for info = 1 and info = 2 both hold.
 \item info = 4  fvec is orthogonal to the columns of the
                jacobian to machine precision.
 \item info = 5  number of calls to fcn has reached or exceeded 200*(n+1).
 \item info = 6  tol is too small. no further reduction in
                the sum of squares is possible.
 \item info = 7  tol is too small.  No further improvement in
                the approximate solution x is possible.
\end{itemize}
         
\item tol : is an optional positive value. Termination occurs when the
      algorithm estimates either that the relative error in the sum of
      squares is at most tol or that the relative error between x and the
      solution is at most tol. If not provided default value is :
               sqrt(epsilon(0.0d0))
\end{itemize}

\subsection*{Example}
Here's a simple example solving the same problem as for the least square polynomial example but using \verb'fvn_lm' instead :
\begin{verbatim}
module excursion
real(8), dimension(:), pointer :: xm => NULL()
real(8), dimension(:), pointer :: ym => NULL()
end module


program gls
use fvn_misc
use excursion
implicit none

    interface
        subroutine zef(m, n, x, fvec, iflag)
            integer(4), intent(in) :: m,n
            real(8), dimension(:), intent(in) :: x
            real(8), dimension(:), intent(inout) :: fvec
            integer(4), intent(inout) :: iflag
        end subroutine
    end interface


 integer,parameter :: npoints=13,deg=3
 integer :: status,i
 real(kind=dp_kind) :: xstep,xc,yc
 real(kind=dp_kind) :: coeff(deg+1)
 real(kind=dp_kind) :: fvec(npoints)
 integer(4) :: iwa(deg+1)
 integer(4) :: info
 real(8) :: tol
 
 allocate(xm(npoints),ym(npoints))
 
 xm = (/ -3.8,-2.7,-2.2,-1.9,-1.1,-0.7,0.5,1.7,2.,2.8,3.2,3.8,4. /)
 ym = (/ -3.1,-2.,-0.9,0.8,1.8,0.4,2.1,1.8,3.2,2.8,3.9,5.2,7.5 /)
 open(2,file='fvn_lsp_double_mesure.dat')
 open(3,file='fvn_lm.dat')
 do i=1,npoints
    write(2,44) xm(i),ym(i)
 end do
 close(2)
 
 coeff=1.
 call fvn_lm(zef,npoints,deg+1,coeff,info)
 write(*,*) "info : ",info
 
 xstep=(xm(npoints)-xm(1))/1000.
 do i=1,1000
    xc=xm(1)+(i-1)*xstep
    yc=poly(xc,coeff)
    write(3,44) xc,yc
 end do
 close(3)
write(*,*) "All done, plot results with gnuplot using command :"
write(*,*) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lm.dat' u 1:2 w l"
44      FORMAT(4(1X,1PE22.14))
contains
function poly(x,coeff)
    use fvn_common
    implicit none
    real(kind=dp_kind) :: x
    real(kind=dp_kind) :: coeff(deg+1)
    real(kind=dp_kind) :: poly
    integer :: i
    poly=0.
    do i=1,deg+1
         poly=poly+coeff(i)*x**(i-1)
    end do
end function
end program

subroutine zef(m, n, x, fvec, iflag)
    use excursion
    implicit none
    integer(4), intent(in) :: m,n
    real(8), dimension(:), intent(in) :: x
    real(8), dimension(:), intent(inout) :: fvec
    integer(4), intent(inout) :: iflag
    
    integer(4) :: i
        
    do i=1,m
        fvec(i)=ym(i)-(x(1)+x(2)*xm(i)+x(3)*xm(i)**2+x(4)*xm(i)**3)
    enddo
end subroutine
\end{verbatim}
Note the need to use a supplementary module \verb'excursion', to allow the subroutine \verb'zef' to have access to the measurement points.

\section{Zero finding}
fvn provides a routine for finding zeros of a complex function using Muller algorithm (only for double complex type). It is based on a version provided on the web by Hans D Mittelmann \url{http://plato.asu.edu/ftp/other\_software/muller.f}.

\begin{verbatim}
 Module : use fvn_misc
 call fvn_muller(f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
\end{verbatim}
\begin{itemize}
 \item f (in) is the complex function (kind=8) for which we search zeros
 \item eps (in) is a real(8) corresponding to the first stopping criterion : let fp(z)=f(z)/p where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) and z(1),...,z(k-1) are previously found roots.  if ((cdabs(f(z)).le.eps) .and. (cdabs(fp(z)).le.eps)), then z is accepted as a root.
 \item eps1 (in) is a real(8) corresponding to the second stopping criterion : a root is accepted if two successive approximations to a given root agree within eps1. Note that if either or both of the stopping criteria are fulfilled, the root is accepted.
 \item kn (in) is an integer equal to the number of known roots, which must be stored in x(1),...,x(kn), prior to entry in the subroutine.
 \item nguess (in) is the number of initial guesses provided. These guesses must be stored in x(kn+1),...,x(kn+nguess). nguess must be set equal to zero if no guesses are provided.
 \item n (in) is an integer equal to the number of new roots to be found.
 \item x (inout) is a complex(8) vector of length kn+n. x(1),...,x(kn) on input must contain any known roots.  x(kn+1),..., x(kn+n) on input may, on user option, contain initial guesses for the n new roots which are to be computed. If the user does not provide an initial guess, zero is used. On output, x(kn+1),...,x(kn+n) contain the approximate roots found by the subroutine.
 \item itmax (in) is an integer equal to the maximum allowable number of iterations per root.
 \item infer (out) is an integer vector of size kn+n. On output infer(j) contains the number of iterations used in finding the j-th root when convergence was achieved.  If convergence was not obtained in itmax iterations, infer(j) will be greater than itmax
 \item ier (out) is an integer used as an error parameter. ier = 33 indicates failure to converge within itmax iterations for at least one of the (n) new roots.
\end{itemize}
This subroutine always returns the last approximation for root j in x(j). if the convergence criterion is satisfied, then infer(j) is less than or equal to itmax. if the convergence criterion is not satisified, then infer(j) is set to either itmax+1 or itmax+k, with k greater than 1. infer(j) = itmax+1 indicates that muller did not obtain convergence in the allowed number of iterations. in this case, the user may wish to set itmax to a larger value. infer(j) = itmax+k means that convergence was obtained (on iteration k) for the deflated function fp(z) = f(z)/((z-z(1)...(z-z(j-1))) but failed for f(z). in this case, better initial guesses might help or, it might be necessary to relax the convergence criterion.

\subsection*{Example}
Example to find the ten roots of $x^{10}-1$
\begin{verbatim}
 program muller
 use fvn_misc
 implicit none

 integer :: i,info
 complex(8),dimension(10) :: roots
 integer,dimension(10) :: infer
 complex(8), external :: f

 call fvn_z_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)

 write(*,*) "Error code :",info
 do i=1,10
    write(*,*) roots(i),infer(i)
 enddo
 end program

 function f(x)
    complex(8) :: x,f
    f=x**10-1
 end function

\end{verbatim}



\section{Numerical integration}
Using an integrated slightly modified version of quadpack \url{http://www.netlib.org/quadpack}, fvn provide adaptative numerical integration (Gauss Kronrod) of real functions of 1 and 2 variables. fvn also provide a function to calculate Gauss-Legendre abscissas and weight, and a simple non adaptative integration subroutine. All routines exists only in fvn for double precision real.

\subsection{Gauss Legendre Abscissas and Weigth}
This subroutine was inspired by Numerical Recipes routine gauleg.
\begin{verbatim}
Module : use fvn_integ
call fvn_gauss_legendre(n,qx,qw)
\end{verbatim}
\begin{itemize}
 \item n (in) is an integer equal to the number of Gauss Legendre points
 \item qx (out) is a real(8) vector of length n containing the abscissas.
 \item qw (out) is a real(8) vector of length n containing the weigths.
\end{itemize}
This subroutine computes n Gauss-Legendre abscissas and weigths

\subsection{Gauss Legendre Numerical Integration}
\begin{verbatim}
Module : use fvn_integ
call fvn_gl_integ(f,a,b,n,res)
\end{verbatim}
\begin{itemize}
 \item f (in) is a real(8) function to integrate
 \item a (in) and b (in) are real(8) respectively lower and higher bound of integration
 \item n (in) is an integer equal to the number of Gauss Legendre points to use
 \item res (out) is a real(8) containing the result
\end{itemize}
This function is a simple Gauss Legendre integration subroutine, which evaluate the integral of function f as in equation \ref{intsple} using n Gauss-Legendre pairs.

\subsection{Gauss Kronrod Adaptative Integration}
This kind of numerical integration is an iterative procedure which try to achieve a given precision.
\subsubsection{Numerical integration of a one variable function}
\begin{verbatim}
Module : use fvn_integ
call fvn_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
\end{verbatim}
This routine evaluate the integral of function f as in equation \ref{intsple}
\begin{itemize}
 \item f (in) is an external real(8) function of one variable
 \item a (in) and b (in) are real(8) respectively lower an higher bound of integration
 \item epsabs (in) and epsrel (in) are real(8) respectively desired absolute and relative error
 \item key (in) is an integer between 1 and 6 correspondind to the Gauss-Kronrod rule to use :
    \begin{itemize}
        \item 1 : 7 - 15 points
        \item 2 : 10 - 21 points
        \item 3 : 15 - 31 points
        \item 4 : 20 - 41 points
        \item 5 : 25 - 51 points
        \item 6 : 30 - 61 points
    \end{itemize}
 \item res (out) is a real(8) containing the estimation of the integration.
 \item abserr (out) is a real(8) equal to the estimated absolute error
 \item ier (out) is an integer used as an error flag
    \begin{itemize}
        \item 0 : no error
        \item 1 : maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yield no improvement it is advised to analyze the integrand in order to determine the integration difficulaties. If the position of a local difficulty can be determined (i.e.singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. If possible, an appropriate special-purpose integrator should be used which is designed for handling the type of difficulty involved.
        \item 2 : the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved.
        \item 3 : extremely bad integrand behaviour occurs at some points of the integration interval.
        \item 6 : the input is invalid, because (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) or limit.lt.1 or lenw.lt.limit*4. result, abserr, neval, last are set to zero. Except when lenw is invalid, iwork(1), work(limit*2+1) and work(limit*3+1) are set to zero, work(1) is set to a and work(limit+1) to b.
    \end{itemize}
 \item limit (in) is an optional integer equal to maximum number of subintervals in the partition of the given integration interval (a,b). If the parameter is not present a default value of 500 will be used.
\end{itemize}

\begin{equation}
 \int_a^b f(x)~dx
 \label{intsple}
\end{equation}




\subsubsection{Numerical integration of a two variable function}
\begin{verbatim}
Module : use fvn_integ
call fvn_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)
\end{verbatim}
This function evaluate the integral of a function f(x,y) as defined in equation \ref{intdble}. The parameters of same name as in the previous paragraph have exactly the same function and behaviour thus only what differs is decribed here
\begin{itemize}
 \item a (in) and b (in) are real(8) corresponding respectively to lower and higher bound of integration for the x variable.
 \item g(x) (in) and h(x) (in) are external functions describing the lower and higher bound of integration for the y variable as a function of x.
\end{itemize}

\begin{equation}
 \int_a^b \int_{g(x)}^{h(x)} f(x,y)~dy~dx
 \label{intdble}
\end{equation}

\subsubsection*{Example}
\begin{verbatim}
program integ
 use fvn_integ
 implicit none

 real(8), external :: f1,f2,g,h
 real(8) :: a,b,epsabs,epsrel,abserr,res
 integer :: key,ier

 a=0.
 b=1.
 epsabs=1d-8
 epsrel=1d-8
 key=2
 call fvn_d_integ_1_gk(f1,a,b,epsabs,epsrel,key,res,abserr,ier,500)
 write(*,*) "Integration of x*x between 0 and 1 : "
 write(*,*) res

 call fvn_d_integ_2_gk(f2,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,500)
 write(*,*) "Integration of x*y between 0 and 1 on both x and y : "
 write(*,*) res
 

end program

function f1(x)
 implicit none
    real(8) :: x,f1
    f1=x*x
end function

function f2(x,y)
 implicit none
    real(8) :: x,y,f2
    f2=x*y
end function

function g(x)
 implicit none
    real(8) :: x,g
    g=0.
end function

function h(x)
 implicit none
    real(8) :: x,h
    h=1.
end function
\end{verbatim}


\section{Special functions}
Specials functions are available in fvn by using an implementation of fnlib \url{http://www.netlib.org/fn} with some additions. This can be used separatly from the rest of fvn by using the module \verb'fvn_fnlib' and linking the library \verb'libfvn_fnlib.a' . The module provides a generic interfaces to all the routines. Specific names of the routines are given in the description.

\begin{verbatim}
 Module : use fvn_fnlib
\end{verbatim}

\paragraph{Important Note}
Due to the addition of fnlib to fvn, some functions that were in fvn and are redondant are now removed from fvn, so update your code now and replace them with the fnlib version. These are listed here after :
\begin{itemize}
 \item \verb'fvn_z_acos' replaced by \verb'acos'
 \item \verb'fvn_z_asin' replaced by \verb'asin'
 \item \verb'fvn_d_asinh' replaced by \verb'asinh'
 \item \verb'fvn_d_acosh' replaced by \verb'acosh'
 \item \verb'fvn_s_csevl' replaced by \verb'csevl'
 \item \verb'fvn_d_csevl' replaced by \verb'csevl'
 \item \verb'fvn_d_factorial' replaced by \verb'fac'
 \item \verb'fvn_d_lngamma' replaced by \verb'alngam'
\end{itemize}


\subsection{Elementary functions}
\subsubsection{carg}
\begin{verbatim}
 carg(z)
\end{verbatim}
\begin{itemize}
 \item z (in) is a complex
\end{itemize}
This function evaluates the argument of the complex z. That is $\theta$ for $z=\rho e^{i\theta}$.

Specific interfaces : \verb'carg,zarg'


\subsubsection{cbrt}
\begin{verbatim}
 cbrt(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the cubic root of the argument x.

Specific interfaces : \verb'cbrt,dcbrt,ccbrt,zcbrt'


\subsubsection{exprl}
\begin{verbatim}
 exprl(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates ${e^x-1}\over x$.

Specific interfaces : \verb'exprel,dexprl,cexprl,zexprl'

\subsubsection{log10}
\begin{verbatim}
 log10(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function is an extension of the intrinsic function log10 to complex arguments.

Specific interfaces : \verb'clog10,zlog10'


\subsubsection{alnrel}
\begin{verbatim}
 alnrel(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates $ln(1+x)$.

Specific interfaces : \verb'alnrel,dlnrel,clnrel,zlnrel'


\subsection{Trigonometry}
\subsubsection{tan}
\begin{verbatim}
 tan(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the tangent of the argument. It is an extension of the intrinsic function tan to complex arguments.

Specific interfaces : \verb'ctan,ztan'


\subsubsection{cot}
\begin{verbatim}
 cot(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluate the cotangent of the argument.

Specific interfaces : \verb'cot,dcot,ccot,zcot'

\subsubsection{sindg}
\begin{verbatim}
 sindg(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluate the sinus of the argument expressed in degrees.

Specific interfaces : \verb'sindg,dsindg'


\subsubsection{cosdg}
\begin{verbatim}
 cosdg(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluate the cosinus of the argument expressed in degrees.

Specific interfaces : \verb'cosdg,dcosdg'


\subsubsection{asin}
\begin{verbatim}
 asin(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the arc sine of the argument. It is an extension of the intrinsic function asin to complex arguments.

Specific interfaces : \verb'casin,zasin'

\subsubsection{acos}
\begin{verbatim}
 acos(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the arc cosine of the argument. It is an extension of the intrinsic function acos to complex arguments.

Specific interfaces : \verb'cacos,zacos'


\subsubsection{atan}
\begin{verbatim}
 atan(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the arc tangent of the argument. It is an extension of the intrinsic function atan to complex arguments.

Specific interfaces : \verb'catan,zatan'

\subsubsection{atan2}
\begin{verbatim}
 atan2(x,y)
\end{verbatim}
\begin{itemize}
 \item x,y are real or complex
\end{itemize}
This function evaluates the arc tangent of $x \over y$. It is an extension of the intrinsic function atan2 to complex arguments.

Specific interfaces : \verb'catan2,zatan2'

\subsubsection{sinh}
\begin{verbatim}
 sinh(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the hyperbolic sine of the argument. It is an extension of the intrinsic function sinh to complex arguments.

Specific interfaces : \verb'csinh,zsinh'


\subsubsection{cosh}
\begin{verbatim}
 cosh(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the hyperbolic cosine of the argument. It is an extension of the intrinsic function cosh to complex arguments.

Specific interfaces : \verb'ccosh,zcosh'

\subsubsection{tanh}
\begin{verbatim}
 tanh(x)
\end{verbatim}
This function evaluates the hyperbolic tangent of the argument. It is an extension of the intrinsic function tanh to complex arguments.

Specific interfaces : \verb'ctanh,ztanh'

\subsubsection{asinh}
\begin{verbatim}
 asinh(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the arc hyperbolic sine of the argument.

Specific interfaces : \verb'asinh,dasinh,casinh,zasinh'

\subsubsection{acosh}
\begin{verbatim}
 acosh(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the arc hyperbolic cosine of the argument.

Specific interfaces : \verb'acosh,dacosh,cacosh,zacosh'

\subsubsection{atanh}
\begin{verbatim}
 atanh(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the arc hyperbolic tangent of the argument.

Specific interfaces : \verb'atanh,datanh,catanh,zatanh'

\subsection{Exponential Integral and related}
\subsubsection{ei}
\begin{verbatim}
 ei(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the exponential integral for argument greater then 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{ei} for $x \neq 0$.
\begin{equation}
\label{ei}
 ei(x)= - \int _{-x} ^\infty {e^{-t}\over t}dt
\end{equation}

Specific interfaces : \verb'ei,dei'


\subsubsection{e1}
\begin{verbatim}
 e1(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
For a real argument, this function evaluates the exponential integral for argument greater than 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{e1} for $x \neq 0$.
\begin{equation}
\label{e1}
 e1(x)= \int _{x} ^\infty {e^{-t}\over t}dt
\end{equation}

For a complex argument, the notation in equation \ref{e1cplx} is used (Abramowitz and Stegun, p.228 \url{http://www.math.ucla.edu/~cbm/aands/page_228.htm}):
\begin{equation}
\label{e1cplx}
 e1(z) = \int _{z} ^\infty {e^{-t}\over t}dt \textrm{~with~} \left|  arg(z) \right| < \pi
\end{equation}
For positive values of real part of $z$, this can be written as in equation \ref{e1cplx_pos} :
\begin{equation}
\label{e1cplx_pos}
 e1(z)= \int _{1} ^\infty {e^{-tz}\over t}dt \textrm{~with~} Re(z) > 0
\end{equation}

 

Specific interfaces : \verb'e1,de1,ze1'

\subsubsection{ali}
\begin{verbatim}
 ali(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the logarithm integral. it is define by equation \ref{ali} for $x > 0$ and $x \neq 1$.
\begin{equation}
 \label{ali}
 ali(x)= - \int _0 ^x {dt \over ln(x)}
\end{equation}

Specific interfaces : \verb'ali,dli'

\subsubsection{si}
\begin{verbatim}
 si(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the sine integral defined by equation \ref{si}.
\begin{equation}
 \label{si}
 si(x)= \int _0 ^x {sin(t) \over t }dt
\end{equation}

Specific interfaces : \verb'si,dsi'


\subsubsection{ci}
\begin{verbatim}
 ci(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the cosine integral defined by equation \ref{ci} where $\gamma \approx 0.57721566$ represent Euler's constant.
\begin{equation}
 \label{ci}
 ci(x)= \gamma + ln(x) + \int _0 ^x {{1-cos(t)} \over t} dt
\end{equation}

Specific interfaces : \verb'ci,dci'

\subsubsection{cin}
\begin{verbatim}
 cin(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the cosine integral alternate definition given by equation \ref{cin}.
\begin{equation}
 \label{cin}
 cin(x)= \int _0 ^x {{1-cos(t)} \over t} dt
\end{equation}

Specific interface : \verb'cin,dcin'

\subsubsection{shi}
\begin{equation}
 shi(x)
\end{equation}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the hyperbolic sine integral defined by equation \ref{shi}.
\begin{equation}
 \label{shi}
 shi(x) = \int _0 ^x {sinh(t) \over t}dt
\end{equation}

Specific interfaces : \verb'shi,dshi'

\subsubsection{chi}
\begin{verbatim}
 chi(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the hyperbolic cosine integral defined by equation \ref{chi} where $\gamma \approx 0.57721566$ represent Euler's constant.
\begin{equation}
 \label{chi}
 chi(x)= \gamma + ln(x) + \int _0 ^x {{cosh(t) -1} \over t}dt
\end{equation}

Specific interfaces : chi,dchi


\subsubsection{cinh}
\begin{verbatim}
 cinh(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the hyperbolic cosine integral alternate definition given by equation \ref{cinh}.
\begin{equation}
 \label{cinh}
 cinh(x) = \int _0 ^x {{cosh(t) -1} \over t}dt
\end{equation}

Specific interfaces : cinh,dcinh


\subsection{Gamma function and related}
\subsubsection{fac}
\begin{verbatim}
 fac(n)
 dfac(n)
\end{verbatim}
\begin{itemize}
 \item n is an integer
\end{itemize}
This function return $n!$ as a real(4) or real(8) for dfac. There's no generic interface for this one.

Specific interfaces : \verb'fac,dfac'

\subsubsection{binom}
\begin{verbatim}
 binom(n,m)
 dbinom(n,m)
\end{verbatim}
\begin{itemize}
 \item n,m are integers
\end{itemize}
This function return the binomial coefficient defined by equation \ref{binom} with $n \geq m \geq 0$. binom returns a real(4), dbinom a real(8). There's no generic interface for this one.
\begin{equation}
 \label{binom}
 binom(n,m) = C_n^m = {{n!} \over {m!(n-m)!}}
\end{equation}

Specific interfaces : \verb'binom,dbinom'


\subsubsection{gamma}
\begin{verbatim}
 gamma(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates $ \Gamma (x) $ defined by equation \ref{gamma}.
\begin{equation}
 \label{gamma}
 \Gamma (x) = \int _0 ^{\infty} t^{x-1}e^{-t}dt
\end{equation}
Note that $n!=\Gamma (n+1)$.

Specific interfaces :\verb'gamma,dgamma,cgamma,zgamm'

\subsubsection{gamr}
\begin{verbatim}
 gamr(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the reciprocal gamma function $gamr(x)= {1 \over \Gamma(x)}$


\subsubsection{alngam}
\begin{verbatim}
 alngam(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates $ln(|\Gamma(x)|)$

Specific interfaces : \verb'alngam,dlngam,clngam,zlngam'


\subsubsection{algams}
\begin{verbatim}
 call algams(x,algam,sgngam)
\end{verbatim}
\begin{itemize}
 \item x (in) is a real
 \item algam (out) is a real
 \item sgngam (out) is a real
\end{itemize}
This subroutine evaluates the logarithm of the absolute value of gamma and the sign of gamma. 
$algam=ln(|\Gamma(x)|)$ and $sgngam=1.0$ or $-1.0$ according to the sign of $\Gamma(x)$.

Specific interfaces : \verb'algams,dlgams'

\subsubsection{gami}
\begin{verbatim}
 gami(a,x)
\end{verbatim}
\begin{itemize}
 \item x is a positive real
 \item a is a strictly positive real
\end{itemize}
This function evaluates the incomplete gamma function defined by equation \ref{gami}.
\begin{equation}
 \label{gami}
 gami(a,x)=\gamma(a,x)=\int _0 ^x t^{a-1} e^{-t}dt
\end{equation}

Specific interfaces : \verb'gami,dgami'

\subsubsection{gamic}
\begin{verbatim}
 gamic(a,x)
\end{verbatim}
\begin{itemize}
 \item x is a positive real
 \item a is a real
\end{itemize}
This function evaluates the complementary incomplete gamma function defined by equation \ref{gamic}.
\begin{equation}
 \label{gamic}
 gamic(a,x)=\Gamma(a,x)=\int _x ^\infty t^{a-1} e^{-t}dt
\end{equation}

Specific interfaces : \verb'gamic,dgamic'

\subsubsection{gamit}
\begin{verbatim}
 gamit(a,x)
\end{verbatim}
\begin{itemize}
 \item x is a positive real
 \item a is a real
\end{itemize}
This function evaluates the Tricomi's incomplete gamma function defined by equation \ref{gamit}.
\begin{equation}
 \label{gamit}
 gamit(a,x)=\gamma^* (a,x)= {{x^{-a}\gamma(a,x)}\over \Gamma(a)}
\end{equation}

Specific interfaces : \verb'gamit,dgamit'


\subsubsection{psi}
\begin{verbatim}
 psi(x)
\end{verbatim}
\begin{itemize}
 \item x is a real or complex
\end{itemize}
This function evaluates the psi function which is the logarithm derivative of the gamma function as defined in equation \ref{psi}.
\begin{equation}
 \label{psi}
  psi(x)= \psi(x) = {d\over dx} ln(\Gamma(x))
\end{equation}
x must not be zero or a negative integer.

Specific interfaces : \verb'psi,dpsi,cpsi,zpsi'


\subsubsection{poch}
\begin{verbatim}
 poch(a,x)
\end{verbatim}
\begin{itemize}
 \item x is a real
 \item a is a real
\end{itemize}
This function evaluates a generalization of Pochhammer's symbol.

Pochhammer's symbol for n a positive integer is given by equation \ref{poch_int}
\begin{equation}
 \label{poch_int}
(a)_n = a(a-1)(a-2)...(a-n+1)
\end{equation}

The generalization of Pochhammer's symbol is given by equation \ref{poch}
\begin{equation}
 \label{poch}
 poch(a,x)= (a)_x = {\Gamma(a+x) \over \Gamma(a)}
\end{equation}

Specific interfaces : \verb'poch,dpoch'


\subsubsection{poch1}
\begin{verbatim}
 poch1(a,x)
\end{verbatim}
\begin{itemize}
 \item x is a real
 \item a is a real
\end{itemize}
This function is defined by equation \ref{poch1}. It is usefull for certains situations, especially when x is small.

\begin{equation}
 \label{poch1}
 poch1(a,x)={{(a)_x-1} \over x}
\end{equation}

Specific interfaces : \verb'poch1,dpoch1'

\subsubsection{beta}
\begin{verbatim}
 beta(a,b)
\end{verbatim}
\begin{itemize}
 \item a,b are real positive or complex
\end{itemize}
This function evaluates $\beta$ function defined by equation \ref{beta}.
\begin{equation}
 \label{beta}
beta(a,b)=\beta(a,b)={  {\Gamma(a) \Gamma(b)} \over \Gamma(a+b) }
\end{equation}

Specific interfaces : \verb'beta,dbeta,cbeta,zbeta'


\subsubsection{albeta}
\begin{verbatim}
 albeta(a,b)
\end{verbatim}
\begin{itemize}
 \item a,b are real positive or complex
\end{itemize}
This function evaluates the natural logarithm of beta function : $ln(\beta(a,b))$

Specific interfaces : \verb'albeta,dlbeta,clbeta,zlbeta'

\subsubsection{betai}
\begin{verbatim}
 betai(x,pin,qin)
\end{verbatim}
\begin{itemize}
 \item x is a real in [0,1]
 \item pin and qin are strictly positive real
\end{itemize}
This function evaluates the incomplete beta function ratio, that is the probability that a random variable from a beta distribution having parameters pin and qin will be less than or equal to x. It is defined by equation \ref{betai}.

\begin{equation}
 \label{betai}
betai(x,pin,qin)=I_x(pin,qin)={1 \over \beta(pin,qin)} \int _0 ^x t^{pin-1}(1-t)^{qin-1}dt
\end{equation}

Specific interfaces : \verb'betai,dbetai'

\subsection{Error function and related}
\subsubsection{erf}
\begin{verbatim}
 erf(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the error function defined by equation \ref{erf}.
\begin{equation}
 \label{erf}
 erf(x)={2\over \sqrt{ \pi}} \int _0 ^x e^{-t^2}dt
\end{equation}

Specific interfaces : \verb'erf,derf'

\subsubsection{erfc}
\begin{verbatim}
 erfc(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the complimentary error function defined by equation \ref{erfc}.
\begin{equation}
 \label{erfc}
 erfc(x)={2\over \sqrt{ \pi}} \int _x ^\infty e^{-t^2}dt
\end{equation}

Specific interfaces : \verb'erfc,derfc'


\subsubsection{daws}
\begin{verbatim}
 daws(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates Dawson's function defined by equation \ref{daws}.
\begin{equation}
 \label{daws}
daws(x)=e^{-x^2} \int _0 ^x e^{t^2}dt
\end{equation}

Specific interfaces : \verb'daws,ddaws'

\subsection{Bessel functions and related}
\subsubsection{bsj0}
\begin{verbatim}
 bsj0(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates Bessel function of the first kind of order 0 defined by equation \ref{bsj0}.
\begin{equation}
 \label{bsj0}
 bsj0(x)=J_0(x)= {1 \over \pi} \int _0 ^\pi cos(x sin(\theta)) d\theta
\end{equation}

Specific interfaces : \verb'besj0,dbesj0'

\subsubsection{bsj1}
\begin{verbatim}
 bsj1(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates Bessel function of the first kind of order 1 defined by equation \ref{bsj1}.
\begin{equation}
 \label{bsj1}
bsj1(x)=J_1(x)={1 \over \pi} \int _0 ^\pi cos(x sin(\theta)-\theta) d\theta
\end{equation}

Specific interfaces : \verb'besj1,dbesj1'

\subsubsection{bsjn}
\begin{verbatim}
 bsjn(n,x,factor,big)
\end{verbatim}
\begin{itemize}
 \item n is an integer
 \item x is a real
 \item factor is an optional integer
 \item big is an optional real
\end{itemize}
This function evaluates Bessel function of the first kind of order n (plotted in figure \ref{besseljnfamily}). These functions satisfy the recurrent relation \ref {bsjn}.
\begin{equation}
 \label{bsjn}
J_{n+1}(x)={{2n}\over x} J_n(x) - J_{n-1}(x)
\end{equation}
This relation is directly used in upward direction to compute $J_n(x)$ for $x>n$. However it is unstable for $x<n$, therefore a Miller's Algorithm is used. The principle of this method is to use the reccurent relation downward from an arbitrary higher than n order with an arbitrary seed and then normalize the solution with \ref{bsjnnorm}
\begin{equation}
 \label{bsjnnorm}
 1=J_0+2J_2+2J_4+2J_6+...
\end{equation}
The optional parameters \verb'factor' and \verb'big' can be used to modify the behaviour of the algorithm. \verb'factor' is used in determining the arbitrary starting order ( an even integer near $n+\sqrt{factor~n}$), the default $factor$ value is 40 for single precision and 150 for double precision. \verb'big' is a real determining the threshold for which anti-overflow counter measures has to be taken, default value is $1.10^{10}$

By convenience, the routine accept $n=0$ and $n=1$, in that cases a call to $bsj0(x)$ or $bsj1(x)$ is actually performed.

\begin{figure}
 \begin{center}
 \includegraphics[width=0.9\textwidth]{besselj-mono.pdf}
 \caption{Bessel $J_n$ functions family}
 \label{besseljnfamily}
 \end{center}
\end{figure}


Specific interfaces : \verb'besjn,dbesjn'


\subsubsection{besrj}
\begin{verbatim}
 call besrj(x,n,b)
\end{verbatim}
\begin{itemize}
 \item x (in) is a real
 \item n (in) is an integer
 \item b (out) is a real array of dimension n
\end{itemize}
This subroutine evaluates Bessel function of the first kind of order \verb'0' to \verb'n-1' for argument x and return the result in array b, which then contain \verb'b(1)='$J_0(x)$,\verb'b(2)='$J_1(x)$,...,\verb'b(n)='$J_{n-1}(x)$.

The algorithm is different from the one used in \verb'bsjn', the choice between the two depends on accuracy and timing considerations, there are test programs (\verb'test_besrj' and \verb'test_bestime') in the \verb'fvn_test' directory which can help choosing the good one. 

Specific interfaces : \verb'besrj,dbesrj'

\subsubsection{bsy0}
\begin{verbatim}
 bsy0(x)
\end{verbatim}
\begin{itemize}
 \item x is a strictly positive real
\end{itemize}
This function evaluates the Bessel function of the second kind of order 0 defined by equation \ref{bsy0}
\begin{equation}
 \label{bsy0}
 bsy0(x)=Y_0(x)={1 \over \pi} \int _0 ^\pi sin(x sin(\theta))d\theta -{2 \over \pi} \int _0 ^\infty e^{-x sinh(t)}dt
\end{equation}

Specific interfaces : \verb'besy0,dbesy0'

\subsubsection{bsy1}
\begin{verbatim}
 bsy1(x)
\end{verbatim}
\begin{itemize}
 \item x is a strictly positive real
\end{itemize}
This function evaluates the Bessel function of the second kind of order 1 defined by equation \ref{bsy1}.
\begin{equation}
 \label{bsy1}
 bsy1(x)=Y_1(x)=-{1 \over \pi} \int _0 ^\pi sin (\theta - x sin(\theta)) d\theta 
- {1 \over \pi} \int _0 ^\infty (e^t -e^{-t})e^{-x sinh(t)}dt
\end{equation}

Specific interfaces : \verb'besy1,dbesy1'


\subsubsection{bsyn}
\begin{verbatim}
 bsyn(n,x)
\end{verbatim}
\begin{itemize}
 \item n is an integer
 \item x is a strictly positive real
\end{itemize}
This function evaluates the Bessel function of the second kind of order n (plotted in figure \ref{besselynfamily}). These functions satisfy the recurrent relation \ref {bsyn}.
\begin{equation}
 \label{bsyn}
Y_{n+1}(x)={{2n}\over x} Y_n(x) - Y_{n-1}(x)
\end{equation}
This recurrent relation is directly used in the upward direction to compute $Y_n(x)$.

By convenience, the routine accept $n=0$ and $n=1$, in that cases a call to $bsy0(x)$ or $bsy1(x)$ is actually performed.

\begin{figure}
 \begin{center}
 \includegraphics[width=0.9\textwidth]{bessely-mono.pdf}
 \caption{Bessel $Y_n$ functions family}
 \label{besselynfamily}
 \end{center}
\end{figure}


Specific interfaces : \verb'besyn,dbesyn'

\subsubsection{bsi0}
\begin{verbatim}
 bsi0(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the Bessel function of the third kind of order 0 defined by equation \ref{bsi0}.
\begin{equation}
 \label{bsi0}
bsi0(x)=I_0(x)={1 \over \pi} \int _0 ^\pi cosh(x cos(\theta))d\theta
\end{equation}

Specific interfaces : \verb'besi0,dbesi0'

\subsubsection{bsi1}
\begin{verbatim}
 bsi1(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the Bessel function of the third kind of order 1 defined by equation \ref{bsi1}.
\begin{equation}
 \label{bsi1}
bsi1(x)=I_1(x)={1 \over \pi} \int _0 ^\pi e^{x cos(\theta)} cos(\theta)d\theta
\end{equation}

Specific interfaces : \verb'besi1,dbesi1'

\subsubsection{bsin}
\begin{verbatim}
 bsin(n,x,factor,big)
\end{verbatim}
\begin{itemize}
 \item n is an integer
 \item x is a real
 \item factor is an optional integer
 \item big is an optional real
\end{itemize}
This function evaluates Bessel function of the third kind of order n (plotted in figure \ref{besselinfamily}). These functions satisfy the recuurent relation \ref{bsin}
\begin{equation}
 \label{bsin}
I_{n+1}(x)=-{{2n}\over x}I_n(x) + I_{n-1}(x)
\end{equation}
This relation is unstable in the upward direction, therefore a Miller's Algorithm is used to evaluate the function. Even if there's a usable normalization relation \ref{bsinnorm}, it is not used in the routine, instead normalization  is done by a simple call to \verb'bsi0(x)'.
\begin{equation}
 \label{bsinnorm}
 1=I_0-2I_2+2I_4-2I_6+....
\end{equation}
The optional parameters \verb'factor' and \verb'big' can be used to modify the behaviour of the algorithm. \verb'factor' is used in determining the arbitrary starting order ( an even integer near $n+\sqrt{factor~n}$), the default $factor$ value is 40 for single precision and 150 for double precision. \verb'big' is a real determining the threshold for which anti-overflow counter measures has to be taken, default value is $1.10^{10}$

By convenience, the routine accept $n=0$ and $n=1$, in that cases a call to $bsi0(x)$ or $bsi1(x)$ is actually performed.

\begin{figure}
 \begin{center}
 \includegraphics[width=0.9\textwidth]{besseli-mono.pdf}
 \caption{Bessel $I_n$ functions family}
 \label{besselinfamily}
 \end{center}
\end{figure}


Specific interfaces : \verb'besin,dbesin'


\subsubsection{besri}
\begin{verbatim}
 call besri(x,n,b)
\end{verbatim}
\begin{itemize}
 \item x (in) is a real
 \item n (in) is an integer
 \item b (out) is a real array of dimension n
\end{itemize}
This subroutine evaluates Bessel function of the third kind of order \verb'0' to \verb'n-1' for argument x and return the result in array b, which then contain \verb'b(1)='$I_0(x)$,\verb'b(2)='$I_1(x)$,...,\verb'b(n)='$I_{n-1}(x)$.

The algorithm is different from the one used in \verb'bsin' and tends to be more accurate, the choice between the two depends on accuracy and timing considerations, there are test programs (\verb'test_besri' and \verb'test_bestime') in the \verb'fvn_test' directory which can help choosing the good one. 

Specific interfaces : \verb'besri,dbesri'


\subsubsection{bsk0}
\begin{verbatim}
 bsk0(x)
\end{verbatim}
\begin{itemize}
 \item x is a strictly positive real
\end{itemize}
This function evaluates the modified Bessel function of the second kind of order 0 defined by equation \ref{bsk0}
\begin{equation}
 \label{bsk0}
bsk0(x)=K_0(x)=\int _0 ^\infty cos(x sinh(t))dt
\end{equation}

Specific interfaces : \verb'besk0,dbesk0'

\subsubsection{bsk1}
\begin{verbatim}
 bsk1(x)
\end{verbatim}
\begin{itemize}
 \item x is a strictly positive real
\end{itemize}
This function evaluates the modified Bessel function of the second kind of order 1 defined by equation \ref{bsk1}
\begin{equation}
 \label{bsk1}
bsk1(x)=K_1(x)=\int _0 ^\infty sin(x sinh(t))sinh(t)dt
\end{equation}

Specific interfaces : \verb'besk1,dbesk1'

\subsubsection{bskn}
\begin{verbatim}
 bskn(n,x)
\end{verbatim}
\begin{itemize}
 \item n is an integer
 \item x is strictly positive real
\end{itemize}
This function evaluates the modified Bessel function of the second kind of order n (plotted in figure \ref{besselknfamily}). These functions satisfy the recurrent relation \ref{bskn}

\begin{equation}
 \label{bskn}
K_{n+1}(x)={{2n}\over x}K_n(x)+K_{n-1}(x)
\end{equation}
This recurrent relation is directly used in the upward direction to compute $K_n(x)$.

By convenience, the routine accept $n=0$ and $n=1$, in that cases a call to $bsk0(x)$ or $bsk1(x)$ is actually performed.

\begin{figure}
 \begin{center}
 \includegraphics[width=0.9\textwidth]{besselk-mono.pdf}
 \caption{Bessel $K_n$ functions family}
 \label{besselknfamily}
 \end{center}
\end{figure}


Specific interface : \verb'beskn,dbeskn'


\subsubsection{bsi0e}
\begin{verbatim}
 bsi0e(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates $e^{-|x|}I_0(x)$

Specific interfaces : \verb'besi0e,dbsi0e'

\subsubsection{bsi1e}
\begin{verbatim}
 bsi1e(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates $e^{-|x|}I_1(x)$

Specific interfaces : \verb'besi1e,dbsi1e'

\subsubsection{bsk0e}
\begin{verbatim}
 bsk0e(x)
\end{verbatim}
\begin{itemize}
 \item x is a strictly positive real
\end{itemize}
This function evaluates $e^x K_0(x)$

Specific interfaces : \verb'besk0e,dbsk0e'

\subsubsection{bsk1e}
\begin{verbatim}
 bsk1e(x)
\end{verbatim}
\begin{itemize}
 \item x is a strictly positive real
\end{itemize}
This function evaluates $e^x K_1(x)$

Specific interfaces : \verb'besk1e,dbsk1e'

\subsubsection{bsks}
\begin{verbatim}
 call bsks(xnu,x,nin,bk)
\end{verbatim}
\begin{itemize}
 \item xnu (in) is a real with $|xnu|<1$. It's the fractional order
 \item x (in) is a real. The value for which the sequence of Bessel functions is to be evaluated.
 \item nin (in) is an integer.
 \item bk (out) is a real vector of length abs(nin), containing the values of the function.
\end{itemize}
This subroutine evaluates a sequence of modified Bessel function of the second kind of fractional order.

If nin is positive, on completion $bk(1)=K_\nu(x)$,$bk(2)=K_{\nu+1}(x)$,...,$bk(nin)=K_{\nu+nin-1}(x)$. If nin is negative, on completion $bk(1)=K_\nu(x)$,$bk(2)=K_{\nu-1}(x)$,...,$bk(|nin|)=K_{\nu+nin+1}(x)$.

Specific interfaces : \verb'besks,dbesks'

\subsubsection{bskes}
\begin{verbatim}
 call bskes(xnu,x,nin,bke)
\end{verbatim}
\begin{itemize}
 \item xnu (in) is a real with $|xnu|<1$. It's the fractional order
 \item x (in) is a real. The value for which the sequence of exponentialy scaled Bessel functions is to be evaluated.
 \item nin (in) is an integer. Number of elements in the sequence.
 \item bke (out) is a real vector of length abs(nin), containing the values of the function.
\end{itemize}
This subroutine evaluates a sequence of exponentially scaled modified Bessel function of the second kind of fractional order.

If nin is positive, on completion $bk(1)=e^x K_\nu(x)$,$bk(2)=e^x K_{\nu+1}(x)$,...,$bk(nin)=e^x K_{\nu+nin-1}(x)$. If nin is negative, on completion $bk(1)=e^x K_\nu(x)$,$bk(2)=e^x K_{\nu-1}(x)$,...,$bk(|nin|)=e^x K_{\nu+nin+1}(x)$.

Specific interfaces : \verb'beskes,dbskes'

\subsection{Airy function and related}
\subsubsection{ai}
\begin{verbatim}
 ai(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the airy function defined by equation \ref{ai}
\begin{equation}
 \label{ai}
Ai(x)={1 \over \pi} \int _0 ^\infty cos(xt+ {1 \over 3}t^3)dt
\end{equation}

Specific interfaces : \verb'ai,dai'

\subsubsection{bi}
\begin{verbatim}
 bi(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the Airy function of the second kind defined by equation \ref{bi}
\begin{equation}
 \label{bi}
Bi(x)={1 \over \pi} \int _0 ^\infty e^{xt- {1 \over 3}t^3}dt +
{1 \over \pi} \int _0 ^\infty sin(xt+ {1 \over 3}t^3)dt
\end{equation}

Specific interfaces : \verb'bi,dbi'


\subsubsection{aid}
\begin{verbatim}
 aid(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the derivative of the Airy function, $aid(x)={d \over {dx}}Ai(x)$.

Specific interface : \verb'aid,daid'

\subsubsection{bid}
\begin{verbatim}
 bid(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the derivative of the Airy function of the second kind, $bid(x)={d \over {dx}}Bi(x)$.

Specific interfaces : \verb'bid,dbid'


\subsubsection{aie}
\begin{verbatim}
 aie(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}

This function evaluates the exponentially scaled Airy function defined in equation \ref{aie}.
%\begin{equation}
% \label{aie}
%aie(x)=Ai(x) \textrm{~if~}x\leq0 \qquad \qquad aie(x)=e^{{2\over3}x^{3\over2}}Ai(x)\rm{~if~}x>0
%\end{equation}
\begin{equation}
\label{aie}
  aie(x)=\left\{\begin{array}{ll}
  Ai(x) & \textrm{~if~}x\leq0 \\
  e^{{2\over3}x^{3\over2}}Ai(x) & \textrm{~if~}x>0
\end{array}\right.
\end{equation}

Specific interfaces : \verb'aie,daie'

\subsubsection{bie}

\begin{verbatim}
 bie(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the exponentially scaled Airy function of the second kind defined in equation \ref{bie}.
%\begin{equation}
% \label{bie}
%bie(x)=Bi(x)\textrm{~if~}x\leq0 \qquad \qquad bie(x)=e^{-{2\over3}x^{3\over2}}Bi(x)\rm{~if~}x>0
%\end{equation}
\begin{equation}
 \label{bie}
bie(x)=\left\{\begin{array}{ll}
Bi(x) & \textrm{~if~}x\leq0 \\
e^{-{2\over3}x^{3\over2}}Bi(x) & \textrm{~if~}x>0
\end{array}\right.
\end{equation}

Specific interfaces : \verb'bie,dbie'

\subsubsection{aide}
\begin{verbatim}
 aide(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the exponentially scaled derivative of the Airy function as defined in equation \ref{aide}.
%\begin{equation}
% \label{aide}
%aie(x)=Ai^\prime(x) \textrm{~if~}x\leq0 \qquad \qquad aie(x)=e^{{2\over3}x^{3\over2}}Ai^\prime(x)\rm{~if~}x>0
%\end{equation}
\begin{equation}
 \label{aide}
aide(x)=\left\{\begin{array}{ll}
Ai^\prime(x) & \textrm{~if~}x\leq0 \\
e^{{2\over3}x^{3\over2}}Ai^\prime(x) & \textrm{~if~}x>0
\end{array}\right.
\end{equation}


Specific interfaces : \verb'aide,daide'


\subsubsection{bide}
\begin{verbatim}
 bide(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates the exponentially scaled derivative of the Airy function of the second kind as defined in equation \ref{bide}.
\begin{equation}
 \label{bide}
bie(x)=\left\{\begin{array}{ll}
Bi^\prime(x) & \textrm{~if~}x\leq0 \\
e^{-{2\over3}x^{3\over2}}Bi^\prime(x) & \textrm{~if~}x>0
\end{array}\right.
\end{equation}

Specific interfaces : \verb'bide,dbide'


\subsection{Miscellanous functions}
\subsubsection{spenc}
\begin{verbatim}
 spenc(x)
\end{verbatim}
\begin{itemize}
 \item x is a real
\end{itemize}
This function evaluates Spence function defined in equation \ref{spenc}.
\begin{equation}
 \label{spenc}
spenc(x)=s(x)=- \int_0^x {{ln(|1-t|)}\over t}dt
\end{equation}

Specific interfaces : \verb'spenc,dspenc'


\subsubsection{inits}
\begin{verbatim}
 inits(os,nos,eta)
\end{verbatim}
\begin{itemize}
 \item os is a real vector of length nos, containing the coefficients in an orthogonal series.
 \item nos is an integer
 \item eta is a real (Warning eta is a real(4) even with the double precision version) representing the requested accuracy.
\end{itemize}
This function initialize the orthogonal series so that inits is the number of terms needed to insure the error is no larger than eta.

Specific interfaces : \verb'inits,initds'


\subsubsection{csevl}
\begin{verbatim}
 csevl(x,cs,n)
\end{verbatim}
\begin{itemize}
 \item x is a real in [-1,1]
 \item cs is a real vector of length n containing the coefficients of the Chebyshev serie.
 \item n is an integer
\end{itemize}
This function evaluates the Chebyshev series whose coefficients are stored in cs.

Specific interfaces : \verb'csevl,dcsevl'




\end{document}