Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module mod_integrate
21 : use const_def, only: dp
22 : use math_lib
23 : use num_def
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : contains
29 :
30 461 : recursive function integrator(func, minx, maxx, args, atol, rtol, max_steps, ierr) result(res)
31 : procedure(integrand) :: func
32 : real(dp),intent(in) :: minx,maxx ! Min and max values to integrate over
33 : real(dp), intent(in) :: args(:) ! Extra args passed to func
34 : real(dp), intent(in) :: atol, rtol ! Absolute and relative tolerances
35 : integer, intent(in) :: max_steps ! Max number of sub-steps
36 : integer, intent(inout) :: ierr ! Error code
37 : real(dp) :: res
38 :
39 461 : real(dp) :: val1, val2
40 461 : real(dp) :: xlow,xhigh,xmid
41 :
42 461 : if(max_steps < 1) then
43 0 : ierr = -1
44 0 : return
45 : end if
46 :
47 461 : ierr = 0
48 :
49 461 : xlow = minx
50 461 : xhigh = maxx
51 461 : xmid = (xhigh+xlow)/2.d0
52 :
53 461 : if(xhigh < xlow) then
54 0 : ierr = -1
55 0 : return
56 : end if
57 :
58 461 : val1 = simp38(func, xlow, xhigh, args, ierr)
59 461 : if(ierr/=0) return
60 :
61 461 : val2 = simp38(func, xlow, xmid, args, ierr) + simp38(func, xmid, xhigh, args, ierr)
62 461 : if(ierr/=0) return
63 :
64 461 : if(val1==0d0 .or. val2 == 0d0) then
65 461 : res = val2
66 : return
67 : end if
68 :
69 433 : if(abs(val1-val2) < atol .or. abs(val1-val2)/val1 < rtol ) then
70 : res = val2
71 : else
72 227 : val1 = integrator(func, xlow, xmid, args, atol, rtol, max_steps-1, ierr)
73 227 : val2 = integrator(func, xmid, xhigh, args, atol, rtol, max_steps-1, ierr)
74 :
75 227 : res = val1+val2
76 : if(ierr/=0) return
77 : end if
78 :
79 : end function integrator
80 :
81 :
82 1383 : real(dp) function simp38(func, minx, maxx, args, ierr)
83 : ! Simpsons 3/8 rule
84 : procedure(integrand) :: func
85 : real(dp),intent(in) :: minx,maxx ! Min and max values to integrate over
86 : real(dp), intent(in) :: args(:) ! Extra args passed to func
87 : integer, intent(inout) :: ierr
88 :
89 : real(dp) :: x
90 :
91 1383 : ierr = 0
92 :
93 1383 : x = minx
94 1383 : simp38 = func(x, args, ierr)
95 1383 : if(ierr/=0) return
96 :
97 1383 : x = (2*minx + maxx)/3.d0
98 1383 : simp38 = simp38 + 3*func(x, args, ierr)
99 1383 : if(ierr/=0) return
100 :
101 1383 : x = (minx + 2*maxx)/3.d0
102 1383 : simp38 = simp38 + 3*func(x, args, ierr)
103 1383 : if(ierr/=0) return
104 :
105 1383 : x = maxx
106 1383 : simp38 = simp38 + func(x, args, ierr)
107 1383 : if(ierr/=0) return
108 :
109 1383 : simp38 = (maxx-minx)/8.d0 * simp38
110 :
111 1383 : end function simp38
112 :
113 :
114 : end module mod_integrate
|
