C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
C     File:        ccubsolv.f (Fortran 77 source for the ccubsolv routine)
C     Author:      Fredrik Jonsson <fj@phys.soton.ac.uk>
C     Date:        December 29, 2005
C     Last change: December 29, 2006
C     Description: The CCUBSOLV(A,Z) routine solves the cubic polynomial
C
C                        z^3+c[2]*z^2+c[1]*z+c[0]=0
C
C                  for general complex coefficients c[k]. The routine
C                  takes a vector A(1..3) of COMPLEX*8 floating point
C                  precision as input, containing the c[k] coefficients
C                  as conforming to the convention
C
C                        A(1)=c[0],  A(2)=c[1],  A(3)=c[2],
C
C                  and returns the three complex-valued roots for z in the
C                  vector Z(1..3), also being of COMPLEX*8 floating point
C                  precision.
C
C     Copyright (C) 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
      subroutine ccubsolv(a,z)
      complex*8 a(3),z(3)
      real*4 pi,phi,onethird,twothird,absrr,absqq,absrr23,absqq32
      complex*8 p,q,rr,qq,cc,w
      integer k

      pi=3.14159265358979323846d0
      onethird=1.0/3.0
      twothird=2.0/3.0
      p=a(2)-onethird*(a(3)**2)
      q=2.0*(onethird*a(3))**3-onethird*a(2)*a(3)+a(1)
      rr=-0.5*q
      qq=onethird*p

C     *********************************************************************
C     The following construction is made in order to avoid numeric overflow
C     in evaluation of the discriminant, in a a manner similar to that used
C     in standard routines for evaluation of sqrt(x^2+y^2), for example.
C     *********************************************************************
      absrr=cabs(rr)
      absqq=cabs(qq)
      absrr23=absrr**twothird
      if ((absrr23).lt.(absqq)) then
         absqq32=absqq**1.5
         cc=rr+absqq32*csqrt((rr/absqq32)**2+(qq/absqq)**3)
      else
         cc=rr+absrr*csqrt((rr/absrr)**2+(qq/absrr23)**3)
      endif

C     *********************************************************************
C     Evaluate the solutions for w from the binomial equation w^3=c, and
C     form the solutions z(k) from the three obtained roots.
C     *********************************************************************
      cc=cc**onethird
      do 10,k=0,2
         phi=k*twothird*pi
         w=cc*cmplx(cos(phi),sin(phi))
         z(k+1)=w-onethird*(a(3)+p/w)
 10   continue

      return
      end