Line data Source code
1 : !---------------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module calc_roots
5 :
6 : implicit none
7 :
8 : public :: cubic_solve, &
9 : quadratic_solve, &
10 : cube_root
11 :
12 : private ! Set Default Scope
13 :
14 : contains
15 :
16 : !=============================================================================
17 0 : function cubic_solve( nz, a_coef, b_coef, c_coef, d_coef ) &
18 0 : result( roots )
19 :
20 : ! Description:
21 : ! Solve for the roots of x in a cubic equation.
22 : !
23 : ! The cubic equation has the form:
24 : !
25 : ! f(x) = a*x^3 + b*x^2 + c*x + d;
26 : !
27 : ! where a /= 0. When f(x) = 0, the cubic formula is used to solve:
28 : !
29 : ! a*x^3 + b*x^2 + c*x + d = 0.
30 : !
31 : ! The cubic formula is also called Cardano's Formula.
32 : !
33 : ! The three solutions for x are:
34 : !
35 : ! x(1) = -(1/3)*(b/a) + ( S + T );
36 : ! x(2) = -(1/3)*(b/a) - (1/2) * ( S + T ) + (1/2)i * sqrt(3) * ( S - T );
37 : ! x(3) = -(1/3)*(b/a) - (1/2) * ( S + T ) - (1/2)i * sqrt(3) * ( S - T );
38 : !
39 : ! where:
40 : !
41 : ! S = ( R + sqrt( D ) )^(1/3); and
42 : ! T = ( R - sqrt( D ) )^(1/3).
43 : !
44 : ! The determinant, D, is given by:
45 : !
46 : ! D = R^2 + Q^3.
47 : !
48 : ! The values of R and Q relate back to the a, b, c, and d coefficients:
49 : !
50 : ! Q = ( 3*(c/a) - (b/a)^2 ) / 9; and
51 : ! R = ( 9*(b/a)*(c/a) - 27*(d/a) - 2*(b/a)^3 ) / 54.
52 : !
53 : ! When D < 0, there are three unique, real-valued roots. When D = 0, there
54 : ! are three real-valued roots, but one root is a double root or a triple
55 : ! root. When D > 0, there is one real-valued root and there are two roots
56 : ! that are complex conjugates.
57 :
58 : ! References:
59 : ! http://mathworld.wolfram.com/CubicFormula.html
60 : !-----------------------------------------------------------------------
61 :
62 : use constants_clubb, only: &
63 : three, & ! Constant(s)
64 : two, &
65 : one_half, &
66 : one_third
67 :
68 : use clubb_precision, only: &
69 : core_rknd ! Variable(s)
70 :
71 : implicit none
72 :
73 : ! Input Variables
74 : integer, intent(in) :: &
75 : nz ! Number of vertical levels
76 :
77 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
78 : a_coef, & ! Coefficient a (of x^3) in a*x^3 + b*x^2 + c^x + d = 0 [-]
79 : b_coef, & ! Coefficient b (of x^2) in a*x^3 + b*x^2 + c^x + d = 0 [-]
80 : c_coef, & ! Coefficient c (of x) in a*x^3 + b*x^2 + c^x + d = 0 [-]
81 : d_coef ! Coefficient d in a*x^3 + b*x^2 + c^x + d = 0 [-]
82 :
83 : ! Return Variables
84 : complex( kind = core_rknd ), dimension(nz,3) :: &
85 : roots ! Roots of x that satisfy a*x^3 + b*x^2 + c*x + d = 0 [-]
86 :
87 : ! Local Variables
88 : real( kind = core_rknd ), dimension(nz) :: &
89 0 : cap_Q_coef, & ! Coefficient Q in cubic formula [-]
90 0 : cap_R_coef, & ! Coefficient R in cubic formula [-]
91 0 : determinant ! Determinant D in cubic formula [-]
92 :
93 : complex( kind = core_rknd ), dimension(nz) :: &
94 0 : sqrt_det, & ! Square root of determinant D in cubic formula [-]
95 0 : cap_S_coef, & ! Coefficient S in cubic formula [-]
96 0 : cap_T_coef ! Coefficient T in cubic formula [-]
97 :
98 : complex( kind = core_rknd ), parameter :: &
99 : i_cmplx = ( 0.0_core_rknd, 1.0_core_rknd ) ! i = sqrt(-1)
100 :
101 : complex( kind = core_rknd ) :: &
102 : sqrt_3, & ! Sqrt 3 (complex data type)
103 : one_half_cmplx, & ! 1/2 (complex data type)
104 : one_third_cmplx ! 1/3 (complex data type)
105 :
106 :
107 : ! Declare some constants as complex data types in order to prevent
108 : ! data-type conversion warning messages.
109 0 : sqrt_3 = cmplx( sqrt( three ), kind = core_rknd )
110 0 : one_half_cmplx = cmplx( one_half, kind = core_rknd )
111 0 : one_third_cmplx = cmplx( one_third, kind = core_rknd )
112 :
113 : ! Find the value of the coefficient Q; where
114 : ! Q = ( 3*(c/a) - (b/a)^2 ) / 9.
115 : cap_Q_coef = ( three * (c_coef/a_coef) - (b_coef/a_coef)**2 ) &
116 0 : / 9.0_core_rknd
117 :
118 : ! Find the value of the coefficient R; where
119 : ! R = ( 9*(b/a)*(c/a) - 27*(d/a) - 2*(b/a)^3 ) / 54.
120 : cap_R_coef = ( 9.0_core_rknd * (b_coef/a_coef) * (c_coef/a_coef) &
121 : - 27.0_core_rknd * (d_coef/a_coef) &
122 0 : - two * (b_coef/a_coef)**3 ) / 54.0_core_rknd
123 :
124 : ! Find the value of the determinant D; where
125 : ! D = R^2 + Q^3.
126 0 : determinant = cap_Q_coef**3 + cap_R_coef**2
127 :
128 : ! Calculate the square root of the determinant. This will be a complex
129 : ! number.
130 0 : sqrt_det = sqrt( cmplx( determinant, kind = core_rknd ) )
131 :
132 : ! Find the value of the coefficient S; where
133 : ! S = ( R + sqrt( D ) )^(1/3).
134 : cap_S_coef &
135 0 : = ( cmplx( cap_R_coef, kind = core_rknd ) + sqrt_det )**one_third_cmplx
136 :
137 : ! Find the value of the coefficient T; where
138 : ! T = ( R - sqrt( D ) )^(1/3).
139 : cap_T_coef &
140 0 : = ( cmplx( cap_R_coef, kind = core_rknd ) - sqrt_det )**one_third_cmplx
141 :
142 : ! Find the values of the roots.
143 : ! This root is always real-valued.
144 : ! x(1) = -(1/3)*(b/a) + ( S + T ).
145 : roots(:,1) &
146 : = - one_third_cmplx * cmplx( b_coef/a_coef, kind = core_rknd ) &
147 0 : + ( cap_S_coef + cap_T_coef )
148 :
149 : ! This root is real-valued when D < 0 (even though the square root of the
150 : ! determinant is a complex number), as well as when D = 0 (when it is part
151 : ! of a double or triple root). When D > 0, this root is a complex number.
152 : ! It is the complex conjugate of roots(3).
153 : ! x(2) = -(1/3)*(b/a) - (1/2) * ( S + T ) + (1/2)i * sqrt(3) * ( S - T ).
154 : roots(:,2) &
155 : = - one_third_cmplx * cmplx( b_coef/a_coef, kind = core_rknd ) &
156 : - one_half_cmplx * ( cap_S_coef + cap_T_coef ) &
157 0 : + one_half_cmplx * i_cmplx * sqrt_3 * ( cap_S_coef - cap_T_coef )
158 :
159 : ! This root is real-valued when D < 0 (even though the square root of the
160 : ! determinant is a complex number), as well as when D = 0 (when it is part
161 : ! of a double or triple root). When D > 0, this root is a complex number.
162 : ! It is the complex conjugate of roots(2).
163 : ! x(3) = -(1/3)*(b/a) - (1/2) * ( S + T ) - (1/2)i * sqrt(3) * ( S - T ).
164 : roots(:,3) &
165 : = - one_third_cmplx * cmplx( b_coef/a_coef, kind = core_rknd ) &
166 : - one_half_cmplx * ( cap_S_coef + cap_T_coef ) &
167 0 : - one_half_cmplx * i_cmplx * sqrt_3 * ( cap_S_coef - cap_T_coef )
168 :
169 :
170 0 : return
171 :
172 0 : end function cubic_solve
173 :
174 : !=============================================================================
175 0 : function quadratic_solve( nz, a_coef, b_coef, c_coef ) &
176 0 : result( roots )
177 :
178 : ! Description:
179 : ! Solve for the roots of x in a quadratic equation.
180 : !
181 : ! The equation has the form:
182 : !
183 : ! f(x) = a*x^2 + b*x + c;
184 : !
185 : ! where a /= 0. When f(x) = 0, the quadratic formula is used to solve:
186 : !
187 : ! a*x^2 + b*x + c = 0.
188 : !
189 : ! The two solutions for x are:
190 : !
191 : ! x(1) = ( -b + sqrt( b^2 - 4*a*c ) ) / (2*a); and
192 : ! x(2) = ( -b - sqrt( b^2 - 4*a*c ) ) / (2*a).
193 : !
194 : ! The determinant, D, is given by:
195 : !
196 : ! D = b^2 - 4*a*c.
197 : !
198 : ! When D > 0, there are two unique, real-valued roots. When D = 0, there
199 : ! are two real-valued roots, but they are a double root. When D < 0, there
200 : ! there are two roots that are complex conjugates.
201 :
202 : ! References:
203 : !-----------------------------------------------------------------------
204 :
205 : use constants_clubb, only: &
206 : four, & ! Constant(s)
207 : two
208 :
209 : use clubb_precision, only: &
210 : core_rknd ! Variable(s)
211 :
212 : implicit none
213 :
214 : ! Input Variables
215 : integer, intent(in) :: &
216 : nz ! Number of vertical levels
217 :
218 : real( kind = core_rknd ), dimension(nz), intent(in) :: &
219 : a_coef, & ! Coefficient a (of x^2) in a*x^2 + b*x + c = 0 [-]
220 : b_coef, & ! Coefficient b (of x) in a*x^2 + b*x + c = 0 [-]
221 : c_coef ! Coefficient c in a*x^2 + b*x + c = 0 [-]
222 :
223 : ! Return Variables
224 : complex( kind = core_rknd ), dimension(nz,2) :: &
225 : roots ! Roots of x that satisfy a*x^2 + b*x + c = 0 [-]
226 :
227 : ! Local Variables
228 : real( kind = core_rknd ), dimension(nz) :: &
229 0 : determinant ! Determinant D in quadratic formula [-]
230 :
231 : complex( kind = core_rknd ), dimension(nz) :: &
232 0 : sqrt_det ! Square root of determinant D in quadratic formula [-]
233 :
234 :
235 : ! Find the value of the determinant D; where
236 : ! D = b^2 - 4*a*c.
237 0 : determinant = b_coef**2 - four * a_coef * c_coef
238 :
239 : ! Calculate the square root of the determinant. This will be a complex
240 : ! number.
241 0 : sqrt_det = sqrt( cmplx( determinant, kind = core_rknd ) )
242 :
243 : ! Find the values of the roots.
244 : ! This root is real-valued when D > 0, as well as when D = 0 (when it is
245 : ! part of a double root). When D < 0, this root is a complex number. It is
246 : ! the complex conjugate of roots(2).
247 : ! x(1) = ( -b + sqrt( b^2 - 4*a*c ) ) / (2*a); and
248 : roots(:,1) = ( -cmplx( b_coef, kind = core_rknd ) + sqrt_det ) &
249 0 : / cmplx( two * a_coef, kind = core_rknd )
250 :
251 : ! This root is real-valued when D > 0, as well as when D = 0 (when it is
252 : ! part of a double root). When D < 0, this root is a complex number. It is
253 : ! the complex conjugate of roots(1).
254 : ! x(2) = ( -b - sqrt( b^2 - 4*a*c ) ) / (2*a).
255 : roots(:,2) = ( -cmplx( b_coef, kind = core_rknd ) - sqrt_det ) &
256 0 : / cmplx( two * a_coef, kind = core_rknd )
257 :
258 :
259 0 : return
260 :
261 0 : end function quadratic_solve
262 :
263 : !=============================================================================
264 0 : function cube_root( x )
265 :
266 : ! Description:
267 : ! Calculates the cube root of x.
268 : !
269 : ! When x >= 0, this code simply calculates x^(1/3). When x < 0, this code
270 : ! uses x^(1/3) = -|x|^(1/3). This eliminates numerical errors when the
271 : ! exponent of 1/3 is not treated as exactly 1/3, which would sometimes
272 : ! result in values of NaN.
273 : !
274 : ! References:
275 : !-----------------------------------------------------------------------
276 :
277 : use constants_clubb, only: &
278 : one_third, & ! Constant(s)
279 : zero
280 :
281 : use clubb_precision, only: &
282 : core_rknd ! Variable(s)
283 :
284 : implicit none
285 :
286 : ! Input Variables
287 : real( kind = core_rknd ), intent(in) :: &
288 : x ! Variable x
289 :
290 : ! Return Variables
291 : real( kind = core_rknd ) :: &
292 : cube_root ! Cube root of x
293 :
294 :
295 0 : if ( x >= zero ) then
296 0 : cube_root = x**one_third
297 : else ! x < 0
298 0 : cube_root = -abs(x)**one_third
299 : endif ! x >= 0
300 :
301 :
302 : return
303 :
304 : end function cube_root
305 :
306 : !===============================================================================
307 :
308 : end module calc_roots
|