C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+ C File: ccubsolv.f (Fortran 77 source for the ccubsolv routine) C Author: Fredrik Jonsson 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 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