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
|
*DECK XDRKIN
SUBROUTINE XDRKIN(DX,NBX,MLOG,BIV,PASV,XLIMV)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Construct Bickley tables for KI1(X), KI2(X), KI3(X), KI4(X), KI5(X),
* taking into account logarithmic singularities for KI1(X) and KI2(X).
*
*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
*
*Parameters: input
* DX step for tables (here, DX=0.02d0).
* NBX nb element in tables (here, NBX=600).
* MLOG interval for logarithmic singularities (suggested values:
* MLOG(1)=30, MLOG(2)=15, MLOG(3)= 0, MLOG(4)= 0, MLOG(5)= 0).
*
*Parameters: output
* BIV elements of quadratic BICKLEY table.
* PASV step quadratic of BICKLEY table.
* XLMV upper limit of quadratic BICKLEY table.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
***** CALLS: *AKIN10* ROUTINE FOR ACCURATE KIN(X) BICKLEY VALUES
* *AK0BES* ROUTINE FOR ACCURATE K0(X) BESSEL VALUES
* *AK1BES* ROUTINE FOR ACCURATE K1(X) BESSEL VALUES
*----
* SUBROUTINE ARGUMENTS
*----
DOUBLE PRECISION DX
INTEGER NBX,MLOG(5)
REAL BIV(0:NBX,3,5),PASV(5),XLIMV(5)
*----
* LOCAL VARIABLES
*----
INTEGER I,KI,IORD
DOUBLE PRECISION X, AKIN(-1:10), AK0BES, AK1BES
DOUBLE PRECISION GAMMA, PIO2, PAS, XLIM
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: C0, C1, C2
PARAMETER ( GAMMA=0.57721566490153D0,
> PIO2= 1.57079632679490D0 )
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(C0(5,0:NBX), C1(5,0:NBX), C2(5,0:NBX))
*
IF( MLOG(3).GT.0.OR.MLOG(4).GT.0.OR.MLOG(5).GT.0 )GOTO 99
*
IF( NBX .NE. 600 )GOTO 97
IF( DX .NE.0.02D0 )GOTO 98
PAS= 1.0D0/DX
XLIM= DBLE(NBX)*DX
DO 10 I=1,5
XLIMV(I)=REAL(XLIM)
PASV(I)=REAL(PAS)
10 CONTINUE
*----
* FIRST, WE CONSTRUCT THE TABLES USING ACCURATE *AKIN10* VALUES
*----
X= 0.D0
CALL AKIN10(X,AKIN(1))
AKIN( 0)= 0.D0
AKIN(-1)= 0.D0
DO 30 I= 0, NBX-1
DO 20 KI= 1, 5
C2(KI,I)= 0.5D0 * AKIN(KI-2)
C1(KI,I)= -(AKIN(KI-1)+X*AKIN(KI-2))
C0(KI,I)= AKIN(KI)+X*(AKIN(KI-1)+X*C2(KI,I))
20 CONTINUE
X= X + DX
CALL AKIN10(X,AKIN(1))
AKIN( 0)= AK0BES(X)
AKIN(-1)= AK1BES(X)
30 CONTINUE
DO 40 KI= 1, 5
C0(KI,NBX)= 0.D0
C1(KI,NBX)= 0.D0
C2(KI,NBX)= 0.D0
40 CONTINUE
*----
* KI1(X) ADJUSTMENTS
*----
X= 0.D0
DO 50 I= 1, MLOG(1)-1
X= X + DX
C0(1,I)= C0(1,I) + 0.5D0*X
C1(1,I)= C1(1,I) - LOG(X)
C2(1,I)= C2(1,I) - 0.5D0/X
50 CONTINUE
*
C1(1,0)= GAMMA-(LOG(2.D0)+1.D0)
*----
* KI2(X) ADJUSTMENTS
*----
X= 0.D0
DO 60 I= 1, MLOG(2)-1
X= X + DX
C0(2,I)= C0(2,I) + 0.25D0*X*X
C1(2,I)= C1(2,I) - X
C2(2,I)= C2(2,I) + 0.5D0*LOG(X) +0.75D0
60 CONTINUE
*
C1(2,0)= -PIO2
C2(2,0)= 0.5D0*(LOG(2.D0)+1.5D0-GAMMA)
*----
* OUTPUT VALUES
*----
DO 80 I= 0, NBX
DO 70 IORD=1,5
BIV(I,1,IORD)= REAL(C0(IORD,I))
BIV(I,2,IORD)= REAL(C1(IORD,I))
BIV(I,3,IORD)= REAL(C2(IORD,I))
70 CONTINUE
80 CONTINUE
*----
* SCRATCH STORAGE DEALLOCATION
*----
DEALLOCATE(C2, C1, C0)
RETURN
*----
* ERROR SECTION
*----
97 CALL XABORT('XDRKIN: KIN TABLES HAVE MORE THAN 600 ELEMENTS')
98 CALL XABORT('XDRKIN: KIN TABLES HAVE A STEP OF 0.02')
99 CALL XABORT('XDRKIN: NO LOG SINGULARITY TAKEN FOR KI345')
END
|
