!************************************************************** !* 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