Search:

Return to previous page

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   

Return to previous page

Generated by ::viewsrc::

Last modified Wednesday 15 Feb 2023