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
|
*DECK XDREXP
SUBROUTINE XDREXP(DX,NBX,PARAM,E00,E01,E10,E11)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Construct exponential tables for linear interpolation.
*
*Copyright:
* Copyright (C) 2002 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): R. Roy, R. Le Tellier
*
*Parameters: input
* DX step for tables (here, DX=0.02d0).
* NBX order of tables (here, NBX=7936).
*
*Parameters: output
* PARAM exponential table characteristics.
* E00 exponential table.
* E01 exponential table.
* E10 exponential table.
* E11 exponential table.
*
*Comments:
* Modified in order to tabulate (1-exp(-x))/x instead of (1-exp(-x))).
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NBX
DOUBLE PRECISION DX
REAL PARAM(3), E00(0:NBX), E01(0:NBX), E10(0:NBX),
> E11(0:NBX)
*----
* LOCAL VARIABLES
*----
INTEGER I
DOUBLE PRECISION X0, X1, X2, X3, X4, PAS, XLIM, EX00, EX01,
> EX10, EX11
DOUBLE PRECISION DEPS, DREF, DZERO, DONE, TAUDMIN
PARAMETER ( DEPS= 1.D-10, DREF= 1.D0/512.D0,
> DZERO= 0.D0, DONE=1.D0, TAUDMIN=2.D-2 )
*----
* SPACE AND VALUES FOR EXPONENTIAL TABLES
*----
INTEGER MEX1
PARAMETER ( MEX1=7936 )
*
IF( NBX.NE. MEX1 ) GO TO 97
IF( DX .GT. DREF+DEPS .OR. DX .LT. DREF-DEPS ) GO TO 98
PAS= DONE/DX
XLIM= DBLE(NBX)*DX
*----
* WE CONSTRUCT THE LINEAR TABLES USING ACCURATE *EXP* VALUES
*----
X0= DZERO
EX00= DZERO
EX10= DONE
DO 20 I= 0, NBX-1
*
* STORE CONSTANT VALUE:
X1= X0 + DX
IF (X1.LE.TAUDMIN) THEN
X2=0.5D0*X1
X3=X1/3.D0
X4=0.5D0*X2
EX11=DONE-X2*(DONE-X3*(DONE-X4))
ELSE
EX11=(DONE - EXP(-X1))/X1
ENDIF
EX01=(DONE - EXP(-X1))
*
* STORE STEP AND CONSTANT VALUES:
X2=(EX01-EX00)/(X1-X0)
X3=(EX11-EX10)/(X1-X0)
E01(I)= REAL(X2)
E00(I)= REAL(EX00-X0*X2)
E11(I)= REAL(X3)
E10(I)= REAL(EX10-X0*X3)
X0= X1
EX00= EX01
EX10= EX11
20 CONTINUE
PARAM(1)=REAL(PAS)
PARAM(2)=REAL(DX)
PARAM(3)=REAL(XLIM)
E10(NBX)= REAL(1.D0/XLIM)
E11(NBX)= 0.0
E00(NBX)= 1.0
E01(NBX)= 0.0
RETURN
*----
* ERROR SECTION
*----
97 CALL XABORT('XDREXP: EXP LINEAR TABLES HAVE MORE THAN 7936 ELEM'
> //'ENTS')
98 CALL XABORT('XDREXP: EXP LINEAR TABLES HAVE A STEP OF 1/512')
END
|
