-
Notifications
You must be signed in to change notification settings - Fork 7
/
acc_F.f90
executable file
·115 lines (113 loc) · 3.84 KB
/
acc_F.f90
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
!> Determines how many terms have to be included for accurate evaluation of the generalized hypergeometric functions using a continued fraction formula, and whether arbitrary precision arithmetic has to be employed or not
!! \param z argument of the generalized hypergeometric function
!! \param ndp measure for the number of significant digits to be included for cases where arbitrary precision arithmetic is required (i.e. cases with ndp >0)
!! \param nfrac required number of terms to be included in continued fraction formula for an accurate evaluation of the generalized hypergeometric functions
subroutine acc_F(z,ndp,nfrac)
implicit none
integer :: ndp,nfrac
real(kind=16) :: z
if(abs(z).lt.0.04_16) then
ndp=0
nfrac=1
else if((abs(z).ge.0.04_16).and.(abs(z).lt.0.1_16)) then
ndp=0
nfrac=2
else if((abs(z).ge.0.1_16).and.(abs(z).lt.1.0_16)) then
ndp=0
nfrac=5
else if((abs(z).ge.1.0_16).and.(abs(z).lt.5.0_16)) then
ndp=0
nfrac=9
else if((abs(z).ge.5.0_16).and.(abs(z).lt.10.0_16)) then
ndp=0
nfrac=12
else if((abs(z).ge.10.0_16).and.(abs(z).lt.30.0_16)) then
ndp=0
nfrac=19
else if((abs(z).ge.30.0_16).and.(abs(z).lt.50.0_16)) then
ndp=0
nfrac=23
else if((abs(z).ge.50.0_16).and.(abs(z).lt.100.0_16)) then
ndp=0
nfrac=31
else if((abs(z).ge.100.0_16).and.(abs(z).lt.200.0_16)) then
ndp=0
nfrac=43
else if((abs(z).ge.200.0_16).and.(abs(z).lt.400.0_16)) then
ndp=0
nfrac=57
else if((abs(z).ge.400.0_16).and.(abs(z).lt.600.0_16)) then
ndp=0
nfrac=70
else if((abs(z).ge.600.0_16).and.(abs(z).lt.800.0_16)) then
ndp=30
nfrac=82
else if((abs(z).ge.800.0_16).and.(abs(z).lt.1000.0_16)) then
ndp=30
nfrac=91
else if((abs(z).ge.1000.0_16).and.(abs(z).lt.1500.0_16)) then
ndp=30
nfrac=110
else if((abs(z).ge.1500.0_16).and.(abs(z).lt.2000.0_16)) then
ndp=44
nfrac=126
else if((abs(z).ge.2000.0_16).and.(abs(z).lt.3000.0_16)) then
ndp=44
nfrac=153
else if((abs(z).ge.3000.0_16).and.(abs(z).lt.4000.0_16)) then
ndp=58
nfrac=176
else if((abs(z).ge.4000.0_16).and.(abs(z).lt.6000.0_16)) then
ndp=73
nfrac=215
else if((abs(z).ge.6000.0_16).and.(abs(z).lt.8000.0_16)) then
ndp=73
nfrac=247
else if((abs(z).ge.8000.0_16).and.(abs(z).lt.10000.0_16)) then
ndp=87
nfrac=276
else if((abs(z).ge.10000.0_16).and.(abs(z).lt.14000.0_16)) then
ndp=102
nfrac=326
else if((abs(z).ge.14000.0_16).and.(abs(z).lt.20000.0_16)) then
ndp=131
nfrac=388
else if((abs(z).ge.20000.0_16).and.(abs(z).lt.30000.0_16)) then
ndp=159
nfrac=474
else if((abs(z).ge.30000.0_16).and.(abs(z).lt.40000.0_16)) then
ndp=174
nfrac=547
else if((abs(z).ge.40000.0_16).and.(abs(z).lt.50000.0_16)) then
ndp=203
nfrac=611
else if((abs(z).ge.50000.0_16).and.(abs(z).lt.70000.0_16)) then
ndp=232
nfrac=723
else if((abs(z).ge.70000.0_16).and.(abs(z).lt.100000.0_16)) then
ndp=275
nfrac=863
else if((abs(z).ge.100000.0_16).and.(abs(z).lt.150000.0_16)) then
ndp=333
nfrac=1056
else if((abs(z).ge.150000.0_16).and.(abs(z).lt.200000.0_16)) then
ndp=391
nfrac=1219
else if((abs(z).ge.200000.0_16).and.(abs(z).lt.300000.0_16)) then
ndp=477
nfrac=1492
else if((abs(z).ge.300000.0_16).and.(abs(z).lt.400000.0_16)) then
ndp=550
nfrac=1722
else if((abs(z).ge.400000.0_16).and.(abs(z).lt.500000.0_16)) then
ndp=622
nfrac=1925
else if((abs(z).ge.500000.0_16).and.(abs(z).lt.700000.0_16)) then
ndp=723
nfrac=2277
else
write(*,*) 'Bessel argument very large!'
write(*,*) 'Increase precision manually to ndp>477 and nfrac>1492'
stop
endif
end subroutine acc_F