*
* $Id: expint.F,v 1.1.1.1 1996/02/15 17:49:10 mclareni Exp $
*
* $Log: expint.F,v $
* Revision 1.1.1.1  1996/02/15 17:49:10  mclareni
* Kernlib
*
*
#include "kernnum/pilot.h"
      REAL FUNCTION EXPINT(RX)
      REAL RX
      CHARACTER*6 ENAME
      LOGICAL MFLAG,RFLAG
#if defined(CERNLIB_NUMHIPRE)
      REAL P1(5),Q1(5),P2(7),Q2(7),P3(6),Q3(6),P4(8),Q4(8)
      REAL A1(8),B1(8),A2(8),B2(8),A3(6),B3(6),XL(6)
      REAL DEXPIN,X,Y,ZERO,ONE,TWO,THREE,AP,BP,DP,AQ,BQ,DQ,X0,V,D
#endif
#if defined(CERNLIB_NUMLOPRE)
      DOUBLE PRECISION P1(5),Q1(5),P2(7),Q2(7),P3(6),Q3(6),P4(8),Q4(8)
      DOUBLE PRECISION A1(8),B1(8),A2(8),B2(8),A3(6),B3(6),XL(6)
      DOUBLE PRECISION X,Y,ZERO,ONE,TWO,THREE,AP,BP,DP,AQ,BQ,DQ,X0,V,D
      DOUBLE PRECISION DEXPIN,DX
#endif
 
      DATA ZERO /0.0D0/, ONE /1.0D0/, TWO /2.0D0/, THREE /3.0D0/
      DATA  X0 /0.37250 74107 8137D0/
      DATA XL /-24.0D0,-12.0D0,-6.0D0,0.0D0,1.0D0,4.0D0/
      DATA P1
     1/+4.29312 52343 210D+0, +3.98941 53870 321D+1,
     2 +2.92525 18866 921D+2, +4.25696 82638 592D+2,
     3 -4.34981 43832 952D+2/
      DATA Q1
     1/+1.00000 00000 000D+0, +1.88992 88395 003D+1,
     2 +1.50950 38744 251D+2, +5.68052 52718 987D+2,
     3 +7.53585 64359 843D+2/
      DATA P2
     1/+4.30967 83946 939D-1, +6.90522 52278 444D+0,
     2 +2.30192 55939 133D+1, +2.43784 08879 132D+1,
     3 +9.04161 55694 633D+0, +9.99979 57705 159D-1,
     4 +4.65627 10797 510D-7/
      DATA Q2
     1/+1.03400 13040 487D-1, +3.31909 21359 330D+0,
     2 +2.04494 78501 379D+1, +4.12807 84189 142D+1,
     3 +3.24264 21069 514D+1, +1.00411 64382 905D+1,
     4 +1.00000 00000 000D+0/
      DATA P3
     1/-2.39099 64453 136D+0, -1.47982 19500 504D+2,
     2 -2.54376 33976 890D+2, -1.19557 61038 372D+2,
     3 -1.96304 08535 939D+1, -9.99999 99990 360D-1/
      DATA Q3
     1/+1.77600 70940 351D+2, +5.30685 09610 812D+2,
     2 +4.62230 27156 148D+2, +1.56818 43364 539D+2,
     3 +2.16304 08494 238D+1, +1.00000 00000 000D+0/
      DATA P4
     1/-8.66937 33995 107D+0, -5.49142 26552 109D+2,
     2 -4.21001 61535 707D+3, -2.49301 39345 865D+5,
     3 -1.19623 66934 925D+5, -2.21744 62775 885D+7,
     4 +3.89280 42131 120D+6, -3.91546 07380 910D+8/
      DATA Q4
     1/+3.41718 75000 000D+1, -1.60708 92658 722D+3,
     2 +3.57300 29805 851D+4, -4.83547 43616 216D+5,
     3 +4.28559 62461 175D+6, -2.49033 37574 054D+7,
     4 +8.91925 76757 561D+7, -1.65254 29972 521D+8/
      DATA A1
     1/-2.18086 38152 072D+0, -2.19010 23385 488D+1,
     2 +9.30816 38566 217D+0, +2.50762 81129 356D+1,
     3 -3.31842 53199 722D+1, +6.01217 99083 008D+1,
     4 -4.32531 13287 813D+1, +1.00443 10922 808D+0/
      DATA B1
     1/+0.00000 00000 000D+0, +3.93707 70185 272D+0,
     2 +3.00892 64837 292D+2, -6.25041 16167 188D+0,
     3 +1.00367 43951 673D+3, +1.43256 73812 194D+1,
     4 +2.73624 11988 933D+3, +5.27468 85196 291D-1/
      DATA A2
     1/-3.48334 65360 285D+0, -1.86545 45488 340D+1,
     2 -8.28561 99414 064D+0, -3.23467 33030 540D+1,
     3 +1.79601 68876 925D+1, +1.75656 31546 961D+0,
     4 -1.95022 32128 966D+0, +9.99994 29607 471D-1/
      DATA B2
     1/+0.00000 00000 000D+0, +6.95000 65588 743D+1,
     2 +5.72837 19383 732D+1, +2.57776 38423 844D+1,
     3 +7.60761 14800 773D+2, +2.89516 72792 514D+1,
     4 -3.43942 26689 987D+0, +1.00083 86740 264D+0/
      DATA A3
     1/-2.77809 28934 438D+1, -1.01047 90815 760D+1,
     2 -9.14830 08216 736D+0, -5.02233 17461 851D+0,
     3 -3.00000 77799 358D+0, +1.00000 00000 704D+0/
      DATA B3
     1/+0.00000 00000 000D+0, +1.22399 93926 823D+2,
     2 +2.72761 00778 779D+0, -7.18975 18395 045D+0,
     3 -2.99901 18065 262D+0, +1.99999 99428 260D+0/
#if defined(CERNLIB_NUMHIPRE)
      ROUND(D)  =  D
#endif
#if defined(CERNLIB_NUMLOPRE)
      ROUND(D)  =  SNGL(D+(D-DBLE(SNGL(D))))
#endif
 
      ENAME='EXPINT'
      X=RX
#if defined(CERNLIB_NUMLOPRE)
      GOTO 9
      ENTRY DEXPIN(DX)
      ENAME='DEXPIN'
      X=DX
#endif
 
    9 IF(X .LE. XL(1)) THEN
       AP=A3(1)-X
       DO 1 I = 2,5
    1  AP=A3(I)-X+B3(I)/AP
       Y=(EXP(-X)/X)*(ONE-(A3(6)+B3(6)/AP)/X)
      ELSEIF(X .LE. XL(2)) THEN
       AP=A2(1)-X
       DO 2 I = 2,7
    2     AP=A2(I)-X+B2(I)/AP
       Y=(EXP(-X)/X)*(A2(8)+B2(8)/AP)
      ELSEIF(X .LE. XL(3)) THEN
       AP=A1(1)-X
       DO 3 I = 2,7
    3     AP=A1(I)-X+B1(I)/AP
       Y=(EXP(-X)/X)*(A1(8)+B1(8)/AP)
      ELSEIF(X .LT. XL(4)) THEN
       V=-TWO*(X/THREE+ONE)
       BP=ZERO
       DP=P4(1)
       DO 4 I = 2,8
          AP=BP
          BP=DP
    4     DP=P4(I)-AP+V*BP
       BQ=ZERO
       DQ=Q4(1)
       DO 14 I = 2,8
          AQ=BQ
          BQ=DQ
   14     DQ=Q4(I)-AQ+V*BQ
       Y=-LOG(-X/X0)+(X+X0)*(DP-AP)/(DQ-AQ)
      ELSEIF(X .EQ. XL(4)) THEN
       CALL KERMTR('C337.1',LGFILE,MFLAG,RFLAG)
       IF(MFLAG) THEN
        IF(LGFILE .EQ. 0) THEN
         WRITE(*,100) ENAME
        ELSE
         WRITE(LGFILE,100) ENAME
        ENDIF
       ENDIF
       IF(.NOT.RFLAG) CALL ABEND
       IF(ENAME .EQ. 'EXPINT')  THEN
          EXPINT=ZERO
       ELSE
          DEXPIN=ZERO
       ENDIF
       RETURN
      ELSEIF(X .LT. XL(5)) THEN
       AP=P1(1)
       AQ=Q1(1)
       DO 5 I = 2,5
          AP=P1(I)+X*AP
    5     AQ=Q1(I)+X*AQ
       Y=-LOG(X)+AP/AQ
      ELSEIF(X .LE. XL(6)) THEN
       Y=ONE/X
       AP=P2(1)
       AQ=Q2(1)
       DO 6 I = 2,7
          AP=P2(I)+Y*AP
    6     AQ=Q2(I)+Y*AQ
       Y=EXP(-X)*AP/AQ
      ELSE
       Y=ONE/X
       AP=P3(1)
       AQ=Q3(1)
       DO 7 I = 2,6
          AP=P3(I)+Y*AP
    7     AQ=Q3(I)+Y*AQ
       Y=EXP(-X)*Y*(ONE+Y*AP/AQ)
      ENDIF
      IF(ENAME .EQ. 'EXPINT')  THEN
         EXPINT=ROUND(Y)
      ELSE
         DEXPIN=Y
      ENDIF
      RETURN
 
  100 FORMAT(7X,A6,' ... ARGUMENT ZERO')
      END