1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
|
*DECK TABEN
FUNCTION TABEN(L,X)
*
*-----------------------------------------------------------------------
*
* COMPUTE AN EXPONENTIAL INTEGRAL OF ORDER L.
*
*----------------------------------------------------------- R. ROY ----
*
IMPLICIT NONE
INTEGER IOUT,NA,NB,NC
CHARACTER NAMSBR*6
PARAMETER(IOUT=6,NA=6,NB=4,NC=4,NAMSBR='TABEN ')
C----
C ROUTINE PARAMETERS
C-----
INTEGER L
REAL X
C----
C FUNCTION TYPE
C----
REAL TABEN
C----
C LOCAL PARAMETERS
C----
INTEGER ITERM
REAL EX,P,R,S
C----
C DATA
C----
REAL A(6),B(4),C(4)
SAVE A,B,C
DATA A
> /-.5772156649,.99999193,-.24991055,.05519968,
> -.00976004,.00107857/
DATA B
> /8.5733287401,18.0590169730,8.634760825,.2677737343/
DATA C
> /9.5733223454,25.6329561486,21.0996530827,3.9584969228/
IF(X .LE. 0.)THEN
IF(L .LE. 1)THEN
TABEN= 1.0E20
ELSE
TABEN=1.0/REAL(L-1)
ENDIF
ELSE
IF (X .LT. 50.) THEN
EX=EXP(-X)
ELSE
EX=0.0
ENDIF
IF (L .EQ. 0) THEN
TABEN=EX/X
ELSE
IF (X .LE. 1.0) THEN
P=A(1)+X*(A(2)+X*(A(3)+X*(A(4)+X*(A(5)+X*A(6)))))
TABEN=P-LOG(X)
ELSE
R=B(4)+X*(B(3)+X*(B(2)+X*(B(1)+X)))
S=C(4)+X*(C(3)+X*(C(2)+X*(C(1)+X)))
TABEN=R/S*EX/X
ENDIF
DO 100 ITERM=1,L-1
TABEN=(EX-X*TABEN)/REAL(ITERM)
100 CONTINUE
ENDIF
ENDIF
RETURN
END
|
