-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmathtool.f90
126 lines (122 loc) · 4.05 KB
/
mathtool.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
116
117
118
119
120
121
122
123
124
125
126
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
function rombint(f,a,b,tol)
! Rombint returns the integral from a to b of using Romberg integration.
! The method converges provided that f(x) is continuous in (a,b).
! f must be real(hp) and must be declared external in the calling
! routine. tol indicates the desired relative accuracy in the integral.
!
use Precision, only : hp
implicit none
integer, parameter :: MAXITER=20
integer, parameter :: MAXJ=5
dimension g(MAXJ+1)
real(hp) f
external f
real(hp) :: rombint
real(hp), intent(in) :: a,b,tol
integer :: nint, i, k, jmax, j
real(hp) :: h, gmax, error, g, g0, g1, fourj
!
h=0.5d0*(b-a)
gmax=h*(f(a)+f(b))
g(1)=gmax
nint=1
error=1.0d20
i=0
10 i=i+1
if (i.gt.MAXITER.or.(i.gt.5.and.abs(error).lt.tol)) &
go to 40
! Calculate next trapezoidal rule approximation to integral.
g0=0._hp
do 20 k=1,nint
g0=g0+f(a+(k+k-1)*h)
20 continue
g0=0.5d0*g(1)+h*g0
h=0.5d0*h
nint=nint+nint
jmax=min(i,MAXJ)
fourj=1._hp
do 30 j=1,jmax
! Use Richardson extrapolation.
fourj=4._hp*fourj
g1=g0+(g0-g(j))/(fourj-1._hp)
g(j)=g0
g0=g1
30 continue
if (abs(g0).gt.tol) then
error=1._hp-gmax/g0
else
error=gmax
end if
gmax=g0
g(jmax+1)=g0
go to 10
40 rombint=g0
if (i.gt.MAXITER.and.abs(error).gt.tol) then
write(*,*) 'Warning: Rombint failed to converge; '
write (*,*)'integral, error, tol:', rombint,error, tol
end if
end function rombint
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
function rombint2(f,a,b,tol)
!THE SAME AS ABOVE, ONLY BECAUSE THAT UNWANTED 'RECURSIVE USAGE' WOULD CAUSE PROBLEMS
! Rombint returns the integral from a to b of using Romberg integration.
! The method converges provided that f(x) is continuous in (a,b).
! f must be real(hp) and must be declared external in the calling
! routine. tol indicates the desired relative accuracy in the integral.
!
use Precision, only : hp
implicit none
integer, parameter :: MAXITER=20
integer, parameter :: MAXJ=5
dimension g(MAXJ+1)
real(hp) f
external f
real(hp) :: rombint2
real(hp), intent(in) :: a,b,tol
integer :: nint, i, k, jmax, j
real(hp) :: h, gmax, error, g, g0, g1, fourj
!
h=0.5d0*(b-a)
gmax=h*(f(a)+f(b))
g(1)=gmax
nint=1
error=1.0d20
i=0
10 i=i+1
if (i.gt.MAXITER.or.(i.gt.5.and.abs(error).lt.tol)) &
go to 40
! Calculate next trapezoidal rule approximation to integral.
g0=0._hp
do 20 k=1,nint
g0=g0+f(a+(k+k-1)*h)
20 continue
g0=0.5d0*g(1)+h*g0
h=0.5d0*h
nint=nint+nint
jmax=min(i,MAXJ)
fourj=1._hp
do 30 j=1,jmax
! Use Richardson extrapolation.
fourj=4._hp*fourj
g1=g0+(g0-g(j))/(fourj-1._hp)
g(j)=g0
g0=g1
30 continue
if (abs(g0).gt.tol) then
error=1._hp-gmax/g0
else
error=gmax
end if
gmax=g0
g(jmax+1)=g0
go to 10
40 rombint2=g0
if (i.gt.MAXITER.and.abs(error).gt.tol) then
write(*,*) 'Warning: Rombint2 failed to converge; '
write (*,*)'integral, error, tol:', rombint2,error, tol, a, b
call sleep(10)
end if
end function rombint2
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc