Contents of file 'ccubsolv.f':
1 C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
2 C File: ccubsolv.f (Fortran 77 source for the ccubsolv routine)
3 C Author: Fredrik Jonsson <fj@phys.soton.ac.uk>
4 C Date: December 29, 2005
5 C Last change: December 29, 2006
6 C Description: The CCUBSOLV(A,Z) routine solves the cubic polynomial
7 C
8 C z^3+c[2]*z^2+c[1]*z+c[0]=0
9 C
10 C for general complex coefficients c[k]. The routine
11 C takes a vector A(1..3) of COMPLEX*8 floating point
12 C precision as input, containing the c[k] coefficients
13 C as conforming to the convention
14 C
15 C A(1)=c[0], A(2)=c[1], A(3)=c[2],
16 C
17 C and returns the three complex-valued roots for z in the
18 C vector Z(1..3), also being of COMPLEX*8 floating point
19 C precision.
20 C
21 C Copyright (C) 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
22 C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
23 subroutine ccubsolv(a,z)
24 complex*8 a(3),z(3)
25 real*4 pi,phi,onethird,twothird,absrr,absqq,absrr23,absqq32
26 complex*8 p,q,rr,qq,cc,w
27 integer k
28
29 pi=3.14159265358979323846d0
30 onethird=1.0/3.0
31 twothird=2.0/3.0
32 p=a(2)-onethird*(a(3)**2)
33 q=2.0*(onethird*a(3))**3-onethird*a(2)*a(3)+a(1)
34 rr=-0.5*q
35 qq=onethird*p
36
37 C *********************************************************************
38 C The following construction is made in order to avoid numeric overflow
39 C in evaluation of the discriminant, in a a manner similar to that used
40 C in standard routines for evaluation of sqrt(x^2+y^2), for example.
41 C *********************************************************************
42 absrr=cabs(rr)
43 absqq=cabs(qq)
44 absrr23=absrr**twothird
45 if ((absrr23).lt.(absqq)) then
46 absqq32=absqq**1.5
47 cc=rr+absqq32*csqrt((rr/absqq32)**2+(qq/absqq)**3)
48 else
49 cc=rr+absrr*csqrt((rr/absrr)**2+(qq/absrr23)**3)
50 endif
51
52 C *********************************************************************
53 C Evaluate the solutions for w from the binomial equation w^3=c, and
54 C form the solutions z(k) from the three obtained roots.
55 C *********************************************************************
56 cc=cc**onethird
57 do 10,k=0,2
58 phi=k*twothird*pi
59 w=cc*cmplx(cos(phi),sin(phi))
60 z(k+1)=w-onethird*(a(3)+p/w)
61 10 continue
62
63 return
64 end
65
Generated by ::viewsrc::