C                                                            <HTML><PRE>
C Subroutine file 'cubic.for' containing subroutine CUBIC designed to
C solve the complex-valued cubic equation.
C
C Version: 8.10
C Date: 2022, November 17
C Coded by Ludek Klimes
C
C.......................................................................
C
      SUBROUTINE CUBIC(AR,AI,BR,BI,CR,CI,X1R,X1I,X2R,X2I,X3R,X3I)
      REAL AR,AI,BR,BI,CR,CI,X1R,X1I,X2R,X2I,X3R,X3I
C
C Subroutine solving complex-valued cubic equation
C X**3+A*X**2+B*X+C=0
C using the algorithm described in Section 5.6 of
C Press, Teukolsky, Vetterling & Flannery (2007).
C
C Input:
C     AR,AI.. Real and imaginary part of coefficient A.
C     BR,BI.. Real and imaginary part of coefficient B.
C     CR,CI.. Real and imaginary part of coefficient C.
C
C Output:
C     X1R,X1I... Real and imaginary part of the root with the smallest
C             real part.
C     X2R,X2I... Real and imaginary part of the root with the middle
C             real part.
C     X3R,X3I... Real and imaginary part of the root with the greatest
C             real part.
C
C No subroutines and external functions required.
C
C Reference:
C     Press, W.H., Teukolsky, S.A., Vetterling, W.V. & Flannery, B.P.:
C       Numerical Recipes.  Third Edition.  Cambridge Univ. Press,
C       New York, 2007, ISBN 978-0-521-88068-8.
C
C Date: 2022, November 17
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Temporary storage locations:
      REAL XR,QR,SR,RR,CCR,AAR,BBR,APBR,AMBR
      REAL XI,QI,SI,RI,CCI,AAI,BBI,APBI,AMBI
      REAL AAM,ARG
C
C.......................................................................
C
C     X=A/3.
      XR=AR/3.
      XI=AI/3.
C     Q=X*X-B/3.  (NR 5.6.10)
      QR=(XR+XI)*(XR-XI)-BR/3.
      QI=2.*XR*XI       -BI/3.
C     R=X*(Q-B/6.)+C/2.  (NR 5.6.10)
      SR=QR-BR/6.
      SI=QI-BI/6.
      RR=XR*SR-XI*SI+0.5*CR
      RI=XR*SI+XI*SR+0.5*CI
C     AA=CSQRT(R*R-Q*Q*Q)  (NR 5.6.13)
      CCR=(RR+RI)*(RR-RI)-QR*(QR*QR-3.*QI*QI)
      CCI=2.*RR*RI       +QI*(QI*QI-3.*QR*QR)
      IF(CCR.EQ.0.AND.CCI.EQ.0.) THEN
C       Double root
        AAR=0.
        AAI=0.
      ELSE IF(CCR.GE.0.) THEN
        AAR=SQRT(0.5*(SQRT(CCR*CCR+CCI*CCI)+CCR))
        AAI=0.5*CCI/AAR
      ELSE
        AAI=SQRT(0.5*(SQRT(CCR*CCR+CCI*CCI)-CCR))
        AAR=0.5*CCI/AAI
      END IF
C     Re(RR*AA).GE.0  (NR 5.6.14)
      IF(RR*AAR+RI*AAI.LT.0.) THEN
        AAR=-AAR
        AAI=-AAI
      END IF
C     AA=-R-CSQRT(R*R-Q*Q*Q)  (NR 5.6.13)
      AAR=-RR-AAR
      AAI=-RI-AAI
      IF(AAR.EQ.0.AND.AAI.EQ.0.) THEN
C       AA=0, then AA**(1./3.)=0 and BB=0  (NR 5.6.16)
        BBR=0.
        BBI=0.
      ELSE
C       AA=AA**(1./3.)  (NR 5.6.13)
        ARG=ATAN2(AAI,AAR)/3.
        AAM=(AAR*AAR+AAI*AAI)**0.166666667
        AAR=AAM*COS(ARG)
        AAI=AAM*SIN(ARG)
C       BB=Q/AA  (NR 5.6.16)
        AAM=AAM*AAM
        BBR=(QR*AAR+QI*AAI)/AAM
        BBI=(QI*AAR-QR*AAI)/AAM
      END IF
C
C     Roots
C     APB=AA+BB  (NR 5.6.17)
      APBR=AAR+BBR
      APBI=AAI+BBI
      X3R=APBR-XR
      X3I=APBI-XI
C     APB=-(AA+BB)/2.-A/3.  (NR 5.6.18)
      APBR=-0.5*APBR-XR
      APBI=-0.5*APBI-XI
C     AMB=SQRT(3/4)*i*(AA-BB)
      AMBR=0.866025403*(BBI-AAI)
      AMBI=0.866025403*(AAR-BBR)
      IF(AMBR.LT.0.) THEN
C       Sorting roots X1 and X2 according to their real parts
        AMBR=-AMBR
        AMBI=-AMBI
      END IF
      X2R=APBR+AMBR
      X2I=APBI+AMBI
      X1R=APBR-AMBR
      X1I=APBI-AMBI
C
C     Sorting root X3 according to its real part
      IF(X3R.LT.X2R) THEN
        XR=X3R
        XI=X3I
        X3R=X2R
        X3I=X2I
        X2R=XR
        X2I=XI
        IF(X2R.LT.X1R) THEN
          XR=X2R
          XI=X2I
          X2R=X1R
          X2I=X1I
          X1R=XR
          X1I=XI
        END IF
      END IF
C
      RETURN
      END
C
C=======================================================================
C                                                                 </PRE>
