solve_cubic.f90

Go to the documentation of this file.

subroutine solve_cubic(c0, c1, c2, c3, x1, x2, x3, ifail)
!-----------------------------------------------------------------------
! solves a cubic equation with a solution with -1.< x < 1
! cn : the coefficient of x**n, x : the real solution with -1.< x < 1.
!-----------------------------------------------------------------------

  use itm_types
  use math_functions, only : root
  
  implicit none
  
  real(r8), intent(in) :: c0, c1, c2, c3
  
  real(r8), intent(out) :: x1, x2, x3
  integer(itm_i4), intent(out) :: ifail
  
  real(r8) :: tol
  real(r8) :: aa, bb, cc, det
  real(r8) :: pi
  real(r8) :: a0, a1, a2
  real(r8) :: p, q
  real(r8) :: u, v
  real(r8) :: angle, dum
  
  x1 = 99._r8
  x2 = 999._r8
  x3 = 9999._r8
  tol = 1.e-8_r8
  ifail = 0
!------------------------------------- 2nd order poly for small c3
  if (abs(c3) / (abs(c1) + abs(c2) + abs(c3)) < 1.e-9_r8) then
    aa = c2
    bb = c1
    cc = c0
    det = bb**2 - 4._r8 * aa * cc
    if (det > 0._r8) then
      x1 = root(aa, bb, cc, det, 1._r8)
      if (abs(x1) > 1._r8 + tol) then
        x1 = root(aa, bb, cc, det, -1._r8)
      end if
    else
      ifail = 1
    end if
  else
!------------------------------------- 3rd order poly solution
    pi = 2._r8 * asin(1._r8)
    a0 = c0 / c3
    a1 = c1 / c3
    a2 = c2 / c3
    p = -(a2**2) / 3._r8 + a1
    q = 2._r8 / 27._r8 * (a2**3) - a2 * a1 / 3._r8 + a0
    det = (p / 3._r8)**3 + (q / 2._r8)**2
    if (det > 0._r8) then
      u = sign(1._r8, -q / 2._r8 + sqrt(det)) * abs(-q / 2._r8 &
       + sqrt(det))**(1._r8 / 3._r8)
      v = sign(1._r8, -q / 2._r8 - sqrt(det)) * abs(-q / 2._r8 &
       - sqrt(det))**(1._r8 / 3._r8)
      x1 =  u + v - a2 / 3._r8
      if (abs(x1) >= 1._r8 + tol) ifail = 1
    else
      p = -p
      angle = sign(1._r8, p) * acos((q / 2._r8) / sqrt(abs(p) &
       / 3._r8)**3)
      x1 = -2._r8 * sqrt(abs(p) / 3._r8) * cos(angle / 3._r8) &
       - a2 / 3._r8
      x2 = -2._r8 * sqrt(abs(p) / 3._r8) * cos(2._r8 * pi &
       / 3._r8 - angle / 3._r8) - a2 / 3._r8
      x3 = -2._r8 * sqrt(abs(p) / 3._r8) * cos(2._r8 * pi &
       / 3._r8 + angle / 3._r8) - a2 / 3._r8
    end if
    if (abs(x1) > abs(x2)) then
      dum = x1
      x1 = x2
      x2 = dum
    end if
    if (abs(x2) > abs(x3)) then
      dum = x2
      x2 = x3
      x3 = dum
    end if
    if (abs(x1) > abs(x2)) then
      dum = x1
      x1 = x2
      x2 = dum
    end if
  end if
  if (abs(x1) > 1._r8 + tol) ifail = 1
  
  return
end subroutine solve_cubic