!**************************************************************
!*   Function sinintegral(x) to integrate the real function   *
!*   FUNC(t) = sin(t)/t from t=0 to t=x by Simpson's method   *
!* ---------------------------------------------------------- *
!*   REFERENCE:  (for Simpson's method)                       * 
!*              "Mathematiques en Turbo-Pascal (Part 1) by    *
!*               Marc Ducamp and Alain Reverchon, Eyrolles,   *
!*               Paris, 1987" [BIBLI 03].                     *
!* ---------------------------------------------------------- *
!*                                                            *
!*                                F90 Version By J-P Moreau.  *
!*                                    (www.jpmoreau.fr)       *
!**************************************************************
Module Sinint

CONTAINS

!Given function to integrate (2 kinds)
real*8 Function FUNC(kind,t)
real*8 t
  if (kind==1) then
    if (dabs(t)<1.d-10) then
      FUNC = 1.d0
    else
      FUNC = dsin(t)/t
    end if
  else
    if (dabs(t)<1.d-10) then
      FUNC = 0.d0
    else
      FUNC = (dcos(t)-1.d0)/t
    end if
  end if
  return
End Function FUNC

!*******************************************************
!* Integral of a function FUNC(X) by Simpson's method  *
!* --------------------------------------------------- *
!* INPUTS:                                             *
!*          kind   =1 for sinintegral                  *
!*                 =2 for cosintegral                  * 
!*          a      begin value of x variable           *
!*          b      end value of x variable             *
!*          n      number of integration steps         *
!*                                                     *
!* OUTPUT:                                             *
!*          res    the integral of FUNC(X) from a to b *
!*                                                     *
!*******************************************************
Subroutine Integral_Simpson(kind, a, b, n, res)
  real*8 a,b,res, step,r
  step=(b-a)/2/n
  r=FUNC(kind,a)
  res=(r+FUNC(kind,b))/2.d0
  do i=1, 2*n-1
    r=FUNC(kind,a+i*step)
    if (Mod(i,2).ne.0) then  
	  res = res + r + r
    else 
	  res = res + r
    end if
  end do
  res = res * step*2/3
  return
End Subroutine Integral_Simpson
          
real*8 Function sinintegral(x)
  real*8 x0,  &   !begin x value
  x               !end x value
  integer nstep   !number of integration steps
  real*8 res      !result of integral

  x0=0.d0
  nstep = int(200*dabs(x))   !this should ensure about 14 exact digits
  
  call Integral_Simpson(1,x0,x,nstep,res)   !kind=1

  sinintegral = res
  return
End Function sinintegral

real*8 Function cosintegral(x)
  real*8 x0,  &   !begin x value
  x               !end x value
  integer nstep   !number of integration steps
  real*8 res      !result of integral

  real*8, parameter :: gamma = 0.57721566490153286d0  !Euler's constant

  x0=0.d0
  nstep = int(200*dabs(x))   !this should ensure about 14 exact digits
  
  call Integral_Simpson(2,x0,x,nstep,res)   !kind=2

  cosintegral = gamma + dlog(x) + res
  return
End Function cosintegral

END MODULE

!End of file sinint.f90