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