Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2020 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 math_lib
21 :
22 : use const_def, only: dp
23 :
24 : use math_pown
25 : use math_io
26 : use utils_lib
27 : use math_def
28 :
29 : use crmath
30 :
31 : implicit none
32 :
33 : character(LEN=16), parameter :: MATH_BACKEND = 'CRMATH'
34 :
35 : real(dp), save :: ln10_m
36 :
37 : interface safe_sqrt
38 : module procedure safe_sqrt_
39 : end interface safe_sqrt
40 :
41 : interface safe_log
42 : module procedure safe_log_
43 : end interface safe_log
44 :
45 : interface safe_log10
46 : module procedure safe_log10_
47 : end interface safe_log10
48 :
49 : interface exp10
50 : module procedure exp10_
51 : end interface exp10
52 :
53 : interface pow
54 : module procedure pow_i_
55 : module procedure pow_r_
56 : end interface pow
57 :
58 : private
59 :
60 : public :: MATH_BACKEND
61 :
62 : public :: math_init
63 : public :: safe_sqrt
64 : public :: log
65 : public :: safe_log
66 : public :: log10
67 : public :: safe_log10
68 : public :: log1p
69 : public :: log2
70 : public :: exp
71 : public :: exp10
72 : public :: expm1
73 : public :: powm1
74 : public :: pow
75 : public :: pow2
76 : public :: pow3
77 : public :: pow4
78 : public :: pow5
79 : public :: pow6
80 : public :: pow7
81 : public :: pow8
82 : public :: cos
83 : public :: sin
84 : public :: tan
85 : public :: cospi
86 : public :: sinpi
87 : public :: tanpi
88 : public :: acos
89 : public :: asin
90 : public :: atan
91 : public :: acospi
92 : public :: asinpi
93 : public :: atanpi
94 : public :: cosh
95 : public :: sinh
96 : public :: tanh
97 : public :: acosh
98 : public :: asinh
99 : public :: atanh
100 : public :: str_to_vector
101 : public :: str_to_double
102 : public :: double_to_str
103 :
104 : contains
105 :
106 16 : subroutine math_init ()
107 :
108 16 : call crmath_init()
109 :
110 16 : ln10_m = log(10._dp)
111 :
112 16 : call precompute_some_zs()
113 :
114 16 : end subroutine math_init
115 :
116 :
117 0 : elemental function safe_sqrt_ (x) result (sqrt_x)
118 :
119 : real(dp), intent(in) :: x
120 : real(dp) :: sqrt_x
121 :
122 0 : sqrt_x = SQRT(MAX(x, 0._dp))
123 :
124 0 : end function safe_sqrt_
125 :
126 :
127 5 : elemental function safe_log_ (x) result (log_x)
128 :
129 : real(dp), intent(in) :: x
130 : real(dp) :: log_x
131 :
132 5 : if (is_nan(x) .or. is_inf(x)) then
133 :
134 : log_x = -99._dp
135 :
136 : else
137 :
138 5 : log_x = log(MAX(1.E-99_dp, x))
139 :
140 : end if
141 :
142 5 : end function safe_log_
143 :
144 :
145 12047 : elemental function safe_log10_ (x) result (log10_x)
146 :
147 : real(dp), intent(in) :: x
148 : real(dp) :: log10_x
149 :
150 12047 : if (is_nan(x) .or. is_inf(x)) then
151 :
152 : log10_x = -99._dp
153 :
154 : else
155 :
156 12047 : log10_x = log10(MAX(1.E-99_dp, x))
157 :
158 : end if
159 :
160 12047 : end function safe_log10_
161 :
162 :
163 11574057 : elemental function exp10_ (x) result (exp10_x)
164 :
165 : real(dp), intent(in) :: x
166 : real(dp) :: exp10_x
167 :
168 : integer :: ix
169 : integer :: i
170 :
171 11574057 : ix = FLOOR(x)
172 :
173 11574057 : if (x == ix) then ! integer power of 10
174 :
175 113043 : exp10_x = 1._dp
176 :
177 920014 : do i = 1, ABS(ix)
178 920014 : exp10_x = exp10_x*10._dp
179 : end do
180 :
181 113043 : if (ix < 0) exp10_x = 1._dp/exp10_x
182 :
183 : else
184 :
185 11461014 : exp10_x = exp(x*ln10_m)
186 :
187 : end if
188 :
189 11574057 : end function exp10_
190 :
191 :
192 89170 : elemental function pow_i_ (x, iy) result (pow_x)
193 :
194 : real(dp), intent(in) :: x
195 : integer, intent(in) :: iy
196 : real(dp) :: pow_x
197 :
198 : integer :: i
199 :
200 89170 : if (x == 0._dp) then
201 :
202 : pow_x = 0._dp
203 :
204 : else
205 :
206 89170 : pow_x = 1._dp
207 :
208 526110 : do i = 1, ABS(iy)
209 526110 : pow_x = pow_x*x
210 : end do
211 :
212 89170 : if (iy < 0) pow_x = 1._dp/pow_x
213 :
214 : end if
215 :
216 89170 : end function pow_i_
217 :
218 :
219 22441851 : elemental function pow_r_ (x, y) result (pow_x)
220 :
221 : real(dp), intent(in) :: x
222 : real(dp), intent(in) :: y
223 : real(dp) :: pow_x
224 :
225 : integer :: iy
226 : integer :: i
227 :
228 22441851 : if (x == 0._dp) then
229 :
230 : pow_x = 0._dp
231 :
232 : else
233 :
234 22441801 : iy = floor(y)
235 :
236 22441801 : if (y == iy .AND. ABS(iy) < 100) then ! integer power of x
237 :
238 733812 : pow_x = 1._dp
239 :
240 3398033 : do i = 1, ABS(iy)
241 3398033 : pow_x = pow_x*x
242 : end do
243 :
244 733812 : if (iy < 0) pow_x = 1._dp/pow_x
245 :
246 : else
247 :
248 21707989 : pow_x = exp(log(x)*y)
249 :
250 : end if
251 :
252 : end if
253 :
254 22441851 : end function pow_r_
255 :
256 :
257 : include 'precompute_zs.inc'
258 :
259 :
260 : end module math_lib
|
