gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ALGEBRA: Make the roots of function quadratic() consistent.

Browse Source 
If there are two roots, make sure that x(1) corresponds to the formula
with '-' sign, and x(2) corresponds to the formula with '+' sign, i.e.

 Δ  = a₁² - 4 a₂ a₀
 x₁ = - (a₁ - sqrt(Δ)) / (2 a₂)
 x₂ = - (a₁ + sqrt(Δ)) / (2 a₂)

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2017-05-02 13:30:30 -03:00
parent 0c7034c5e8
commit 3d4beabc36
1 changed files with 14 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
src/algebra.F90
View File

@ -65,7 +65,8 @@ module algebra
! Arguments:
!
! a₂, a₁, a₀ - the quadratic equation coefficients;
! x - the root array;
! x - the root array; if there are two roots, x(1) corresponds to
! the one with '-' and x(2) to the one with '+';
!
! Return value:
!
@ -127,27 +128,23 @@ module algebra
! Δ > 0, so the quadratic has two real roots
!
if (b(2) /= 0.0d+00) then
tm = - (bh + sign(dr, bh))
if (b(2) > 0.0d+00) then
tm = - bh - dr
x(1) = b(1) / tm
x(2) = tm / b(3)
else if (b(2) < 0.0d+00) then
tm = - bh + dr
x(1) = tm / b(3)
x(2) = b(1) / tm
else
x(2) = dr / b(3)
x(1) = - x(2)
x(1) = dr / b(3)
x(2) = - x(1)
end if
! update the number of roots
!
nr = 2
! sort roots
!
if (x(1) > x(2)) then
tm = x(1)
x(1) = x(2)
x(2) = tm
end if
else if (dl == 0.0d+00) then ! Δ = 0
! Δ = 0, so the quadratic has two identical real roots
@ -180,12 +177,12 @@ module algebra
! (a₂ x + a₁) x = 0
!
tm = - b(2) / b(3)
if (tm < 0.0d+00) then
x(1) = tm
x(2) = 0.0d+00
else
if (tm >= 0.0d+00) then
x(1) = 0.0d+00
x(2) = tm
else
x(1) = tm
x(2) = 0.0d+00
end if
! update the number of roots