Search:

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   ```

```