!******************************************************************
!*      Purpose: This program computes the cosine and sine        *
!*               integrals using subroutine CISIA                 *
!*      Input :  x  --- Argument of Ci(x) and Si(x)               *
!*      Output:  CI --- Ci(x)                                     *
!*               SI --- Si(x)                                     *
!*      Example:                                                  * 
!*                   x         Ci(x)          Si(x)               *
!*                ------------------------------------            *
!*                  0.0     - ì             .00000000             *
!*                  5.0     -.19002975     1.54993124             *
!*                 10.0     -.04545643     1.65834759             *
!*                 20.0      .04441982     1.54824170             * 
!*                 30.0     -.03303242     1.56675654             *
!*                 40.0      .01902001     1.58698512             *
!* -------------------------------------------------------------- *
!* REFERENCE: "Fortran Routines for Computation of Special        *
!*             Functions jin.ece.uiuc.edu/routines/routines.html" *
!*                                                                *
!*                              F90 Release By J-P Moreau, Paris. *
!*                                      (www.jpmoreau.fr)         *
!******************************************************************
        PROGRAM MCISIA
        DOUBLE PRECISION CI,SI,X
        WRITE(*,*)'Please enter x '
        READ(*,*) X
        WRITE(*,*)'   x         Ci(x)          Si(x)'
        WRITE(*,*)'------------------------------------'
        CALL CISIA(X,CI,SI)
        IF (X.NE.0.0D0) WRITE(*,10)X,CI,SI
        IF (X.EQ.0.0D0) WRITE(*,20)
10      FORMAT(1X,F5.1,2F15.8)
20      FORMAT(3X,' .0',4X,' - ì',13X,'.00000000')
        END



        SUBROUTINE CISIA(X,CI,SI)

!       =============================================
!       Purpose: Compute cosine and sine integrals
!                Si(x) and Ci(x)  ( x ò 0 )
!       Input :  x  --- Argument of Ci(x) and Si(x)
!       Output:  CI --- Ci(x)
!                SI --- Si(x)
!       =============================================

        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION BJ(101)
        P2=1.570796326794897D0
        EL=.5772156649015329D0
        EPS=1.0D-15
        X2=X*X
        IF (X.EQ.0.0D0) THEN
           CI=-1.0D+300
           SI=0.0D0
        ELSE IF (X.LE.16.0D0) THEN
           XR=-.25D0*X2
           CI=EL+DLOG(X)+XR
           DO 10 K=2,40
              XR=-.5D0*XR*(K-1)/(K*K*(2*K-1))*X2
              CI=CI+XR
              IF (DABS(XR).LT.DABS(CI)*EPS) GO TO 15
10         CONTINUE
15         XR=X
           SI=X
           DO 20 K=1,40
              XR=-.5D0*XR*(2*K-1)/K/(4*K*K+4*K+1)*X2
              SI=SI+XR
              IF (DABS(XR).LT.DABS(SI)*EPS) RETURN
20         CONTINUE
        ELSE IF (X.LE.32.0D0) THEN
           M=INT(47.2+.82*X)
           XA1=0.0D0
           XA0=1.0D-100
           DO 25 K=M,1,-1
              XA=4.0D0*K*XA0/X-XA1
              BJ(K)=XA
              XA1=XA0
25            XA0=XA
           XS=BJ(1)
           DO 30 K=3,M,2
30            XS=XS+2.0D0*BJ(K)
           BJ(1)=BJ(1)/XS
           DO 35 K=2,M
35            BJ(K)=BJ(K)/XS
           XR=1.0D0
           XG1=BJ(1)
           DO 40 K=2,M
              XR=.25D0*XR*(2.0*K-3.0)**2/((K-1.0)*(2.0*K-1.0)**2)*X
40            XG1=XG1+BJ(K)*XR
           XR=1.0D0
           XG2=BJ(1)
           DO 45 K=2,M
              XR=.25D0*XR*(2.0*K-5.0)**2/((K-1.0)*(2.0*K-3.0)**2)*X
45            XG2=XG2+BJ(K)*XR
           XCS=DCOS(X/2.0D0)
           XSS=DSIN(X/2.0D0)
           CI=EL+DLOG(X)-X*XSS*XG1+2*XCS*XG2-2*XCS*XCS
           SI=X*XCS*XG1+2*XSS*XG2-DSIN(X)
        ELSE
           XR=1.0D0
           XF=1.0D0
           DO 50 K=1,9
              XR=-2.0D0*XR*K*(2*K-1)/X2
50            XF=XF+XR
           XR=1.0D0/X
           XG=XR
           DO 55 K=1,8
              XR=-2.0D0*XR*(2*K+1)*K/X2
55            XG=XG+XR
           CI=XF*DSIN(X)/X-XG*DCOS(X)/X
           SI=P2-XF*DCOS(X)/X-XG*DSIN(X)/X
        ENDIF
        RETURN
        END

!end of file mcisia.f90