Sources Esope | Cast3M

cubic
C CUBIC     SOURCE    CB215821  17/11/30    21:15:48     9639           
      SUBROUTINE CUBIC (CB3, CB2, CB1, CB0, X1, X2, X3, NRoot)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC CCREEL
 
C ... Attention à la valeur arbitraire de 1.d-12 ! ...
      DATA EPSILO /1.D-12/
 
      PI = XPI
      PIDIV2 = PI / 2.d0
      PI2 = 2.d0 * PI
 
      NRoot = 0
      X1 = 0.d0
      X2 = 0.d0
      X3 = 0.d0
 
      A = CB2 / CB3
      B = CB1 / CB3
      C = CB0 / CB3
 
      D  = 1.d0 / 3.d0
      XJ =   A  / 3.d0
 
      IF ( ABS(C) .LE. EPSILO) THEN
C    ... On a alors affaire à une une équation quadratique ...
         CC2 = 1.d0
         CC1 = A
         CC0 = B
         CALL QUADRA(CC2, CC1, CC0, X2, X3, NRoot)
         NROOT=NROOT+1
         RETURN
      ELSE
C    ... Sinon, on élimine le terme x^2 en posant x' = x - A/3 ...
C    ... Nouvelle équation : x^3 + E x + F = 0 ...
         E = B - A*A / 3.d0
         F = A * (2.d0 * A * A - 9.d0 * B) / 27.d0 + C
C    ... G est son déterminant ...
         G = E*E*E  / 27.d0 + F * F / 4.d0
      ENDIF
 
C ... Si G = 0 il y a 3 racines réelles, dont une double ...
      IF ( ABS(G) .LT. EPSILO) THEN
         NRoot = 3
         XI = F / 2.d0
         ALPHA =SIGN(1.D0,-XI) * (ABS(XI)**D)
         X1 =  2.d0 * ALPHA - XJ
         X2 = -1.d0 * ALPHA - XJ
         X3 = X2
         RETURN
      ENDIF
 
C ... Sinon (G != 0) il y a 2 possibilités ...
      IF(SIGN(1.D0,G). LT. 0.) THEN
C    ... G &lt; 0 => 3 racines réelles ...
         NRoot = 3
         H = 2.d0 * SQRT(-E / 3.d0)
         B = -.5d0 * F / SQRT(-E * E * E / 27.d0)
         A = SQRT(1.d0 - B * B)
         IF( ABS(A).LT. EPSILO) THEN
            C = (1.d0 -SIGN(1.D0,B)) * PIDIV2
         ELSE
            C = (2.d0 -SIGN(1.D0,A)) * PIDIV2
            IF( ABS(B).GT. EPSILO) THEN
               C = C +SIGN(1.D0,A)*SIGN(1.D0,B)*
     &             (ABS(ATAN(A / B)) - PIDIV2)
            ENDIF
         ENDIF
         X1 = H * COS(C / 3.d0) - XJ
         X2 = H * COS(C / 3.d0 + PI2 / 3.d0) - XJ
         X3 = H * COS(C / 3.d0 + 2.d0 * PI2 / 3.d0) - XJ
         RETURN
      ELSE IF (SIGN(1.D0,G) .GT. 0.) THEN
C    ... G > 0 => 1 racine réelle + 2 racines complexes conjuguées ...
         NRoot = 1
         XI = F / 2.d0
         XL = -XI + SQRT(G)
         XM = -XI - SQRT(G)
         X1 =SIGN(1.D0,XL)*(ABS(XL)**D) +
     &       SIGN(1.D0,XM)*(ABS(XM)**D) - XJ
      ENDIF
 
      RETURN
      END