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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
DOUBLE PRECISION FUNCTION DEBYE3(XVALUE)
C
C
C DEFINITION:
C
C This program calculates the Debye function of order 3, defined as
C
C DEBYE3(x) = 3*[Integral {0 to x} t^3/(exp(t)-1) dt] / (x^3)
C
C The code uses Chebyshev series whose coefficients
C are given to 20 decimal places.
C
C
C ERROR RETURNS:
C
C If XVALUE < 0.0 an error message is printed and the
C function returns the value 0.0
C
C
C MACHINE-DEPENDENT PARAMETERS:
C
C NTERMS - INTEGER - The no. of elements of the array ADEB3.
C The recommended value is such that
C ABS(ADEB3(NTERMS)) < EPS/100,
C subject to 1 <= NTERMS <= 18
C
C XLOW - DOUBLE PRECISION - The value below which
C DEBYE3 = 1 - 3x/8 + x*x/20 to machine precision.
C The recommended value is
C SQRT(8*EPSNEG)
C
C XUPPER - DOUBLE PRECISION - The value above which
C DEBYE3 = (18*zeta(4)/x^3) - 3*exp(-x)(x^3+3x^2+6x+6)/x^3.
C The recommended value is
C -LOG(2*EPS)
C
C XLIM1 - DOUBLE PRECISION - The value above which DEBYE3 = 18*zeta(4)/x^3
C The recommended value is
C -LOG(XMIN)
C
C XLIM2 - DOUBLE PRECISION - The value above which DEBYE3 = 0.0 to machine
C precision. The recommended value is
C CUBE ROOT(19/XMIN)
C
C For values of EPS, EPSNEG, and XMIN see the file MACHCON.TXT
C
C The machine-dependent constants are computed internally by
C using the D1MACH subroutine.
C
C
C OTHER MISCFUN SUBROUTINES USED:
C
C CHEVAL , ERRPRN, D1MACH
C
C
C INTRINSIC FUNCTIONS USED:
C
C AINT , EXP , INT , LOG , SQRT
C
C
C AUTHOR:
C Dr. Allan J. MacLeod,
C Dept. of Mathematics and Statistics,
C University of Paisley
C High St.
C PAISLEY
C SCOTLAND
C PA1 2BE
C
C (e-mail: macl_ms0@paisley.ac.uk )
C
C
C LATEST UPDATE: 23 January, 1996
C
INTEGER I,NEXP,NTERMS
DOUBLE PRECISION ADEB3(0:18),CHEVAL,DEBINF,EIGHT,EXPMX,FOUR,
& HALF,ONE,ONEHUN,PT375,RK,SEVP5,SIX,SUM,T,THREE,TWENTY,X,
& XK,XKI,XLIM1,XLIM2,XLOW,XUPPER,XVALUE,ZERO,D1MACH
C
C OTHER MISCFUN SUBROUTINES USED:
C
external CHEVAL, D1MACH
C
C
C INTRINSIC FUNCTIONS USED:
C
intrinsic AINT , EXP , INT , LOG , SQRT, abs
C
c*****CHARACTER FNNAME*6,ERRMSG*17
c*****DATA FNNAME/'DEBYE3'/
c*****DATA ERRMSG/'ARGUMENT NEGATIVE'/
DATA ZERO,PT375/0.0 D 0 , 0.375 D 0/
DATA HALF,ONE/0.5 D 0 , 1.0 D 0/
DATA THREE,FOUR,SIX/3.0 D 0 , 4.0 D 0 , 6.0 D 0/
DATA SEVP5,EIGHT,TWENTY/7.5 D 0 , 8.0 D 0 , 20.0 D 0/
DATA ONEHUN/100.0 D 0/
DATA DEBINF/0.51329 91127 34216 75946 D -1/
DATA ADEB3/2.70773 70683 27440 94526 D 0,
1 0.34006 81352 11091 75100 D 0,
2 -0.12945 15018 44408 6863 D -1,
3 0.79637 55380 17381 64 D -3,
4 -0.54636 00095 90823 8 D -4,
5 0.39243 01959 88049 D -5,
6 -0.28940 32823 5386 D -6,
7 0.21731 76139 625 D -7,
8 -0.16542 09994 98 D -8,
9 0.12727 96189 2 D -9,
X -0.98796 3459 D -11,
1 0.77250 740 D -12,
2 -0.60779 72 D -13,
3 0.48075 9 D -14,
4 -0.38204 D -15,
5 0.3048 D -16,
6 -0.244 D -17,
7 0.20 D -18,
8 -0.2 D -19/
C
C Start computation
C
X = XVALUE
C
C Error test
C
IF ( X .LT. ZERO ) THEN
c********CALL ERRPRN(FNNAME,ERRMSG)
DEBYE3 = ZERO
RETURN
ENDIF
C
C Compute the machine-dependent constants.
C
T = D1MACH(1)
XLIM1 = - LOG( T )
XK = ONE / THREE
XKI = (ONE/DEBINF) ** XK
RK = T ** XK
XLIM2 = XKI / RK
T = D1MACH(3)
XLOW = SQRT ( T * EIGHT )
XUPPER = - LOG( T + T )
T = T / ONEHUN
DO 10 NTERMS = 18 , 0 , -1
IF ( ABS(ADEB3(NTERMS)) .GT. T ) GOTO 19
10 CONTINUE
C
C Code for x <= 4.0
C
19 IF ( X .LE. FOUR ) THEN
IF ( X .LT. XLOW ) THEN
DEBYE3 = ( ( X - SEVP5 ) * X + TWENTY ) / TWENTY
ELSE
T = ( ( X * X / EIGHT ) - HALF ) - HALF
DEBYE3 = CHEVAL ( NTERMS , ADEB3 , T ) - PT375 * X
ENDIF
ELSE
C
C Code for x > 4.0
C
IF ( X .GT. XLIM2 ) THEN
DEBYE3 = ZERO
ELSE
DEBYE3 = ONE / ( DEBINF * X * X * X )
IF ( X .LT. XLIM1 ) THEN
EXPMX = EXP ( -X )
IF ( X .GT. XUPPER ) THEN
SUM = (((X+THREE)*X+SIX)*X+SIX) / (X*X*X)
ELSE
SUM = ZERO
RK = AINT ( XLIM1 / X )
NEXP = INT ( RK )
XK = RK * X
DO 100 I = NEXP,1,-1
XKI = ONE / XK
T = (((SIX*XKI+SIX)*XKI+THREE)*XKI+ONE) / RK
SUM = SUM * EXPMX + T
RK = RK - ONE
XK = XK - X
100 CONTINUE
ENDIF
DEBYE3 = DEBYE3 - THREE * SUM * EXPMX
ENDIF
ENDIF
ENDIF
RETURN
END
|
