Esope source files

   1 : C CUBIC     SOURCE    CHAT      05/01/12    22:32:44     5004           
   2 :       SUBROUTINE CUBIC (CB3, CB2, CB1, CB0, X1, X2, X3, NRoot) 
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8(A-H,O-Z)
   5 : -INC CCREEL
   6 : 
   7 : C ... Attention à la valeur arbitraire de 1.d-12 ! ...
   8 :       DATA EPSILO /1.D-12/
   9 : 
  10 :       PI = XPI
  11 :       PIDIV2 = PI / 2.d0
  12 :       PI2 = 2.d0 * PI
  13 : 
  14 :       NRoot = 0
  15 :       X1 = 0.d0
  16 :       X2 = 0.d0
  17 :       X3 = 0.d0
  18 : 
  19 :       A = CB2 / CB3
  20 :       B = CB1 / CB3
  21 :       C = CB0 / CB3
  22 : 
  23 :       D  = 1.d0 / 3.d0
  24 :       XJ =   A  / 3.d0
  25 : 
  26 :       IF ( DABS(C) .LE. EPSILO) THEN
  27 : C    ... On a alors affaire à une une équation quadratique ...
  28 :          CC2 = 1.d0
  29 :          CC1 = A
  30 :          CC0 = B
  30 :          CALL QUADRA(CC2, CC1, CC0, X2, X3, NRoot)
  31 :          NROOT=NROOT+1
  32 :          RETURN
  33 :       ELSE
  34 : C    ... Sinon, on élimine le terme x^2 en posant x' = x - A/3 ...
  35 : C    ... Nouvelle équation : x^3 + E x + F = 0 ...
  36 :          E = B - A*A / 3.d0
  37 :          F = A * (2.d0 * A * A - 9.d0 * B) / 27.d0 + C
  38 : C    ... G est son déterminant ...
  39 :          G = E*E*E  / 27.d0 + F * F / 4.d0
  40 :       ENDIF
  41 : 
  42 : C ... Si G = 0 il y a 3 racines réelles, dont une double ...
  43 :       IF ( DABS(G) .LT. EPSILO) THEN
  44 :          NRoot = 3
  45 :          XI = F / 2.d0
  46 :          ALPHA = DSIGN(1.D0,-XI) * (DABS(XI)**D)
  47 :          X1 =  2.d0 * ALPHA - XJ
  48 :          X2 = -1.d0 * ALPHA - XJ
  49 :          X3 = X2
  50 :          RETURN
  51 :       ENDIF
  52 : 
  53 : C ... Sinon (G != 0) il y a 2 possibilités ...
  54 :       IF(DSIGN(1.D0,G). LT. 0.) THEN
  55 : C    ... G < 0 => 3 racines réelles ...
  56 :          NRoot = 3
  57 :          H = 2.d0 * SQRT(-E / 3.d0)
  58 :          B = -.5d0 * F / SQRT(-E * E * E / 27.d0)
  59 :          A = DSQRT(1.d0 - B * B)
  60 :          IF( DABS(A).LT. EPSILO) THEN
  61 :             C = (1.d0 - DSIGN(1.D0,B)) * PIDIV2
  62 :          ELSE
  63 :             C = (2.d0 - DSIGN(1.D0,A)) * PIDIV2
  64 :             IF( DABS(B).GT. EPSILO) THEN
  65 :                C = C + DSIGN(1.D0,A)*DSIGN(1.D0,B)*
  66 :      &             (ABS(ATAN(A / B)) - PIDIV2)
  67 :             ENDIF
  68 :          ENDIF
  69 :          X1 = H * COS(C / 3.d0) - XJ
  70 :          X2 = H * COS(C / 3.d0 + PI2 / 3.d0) - XJ
  71 :          X3 = H * COS(C / 3.d0 + 2.d0 * PI2 / 3.d0) - XJ
  72 :          RETURN
  73 :       ELSE IF ( DSIGN(1.D0,G) .GT. 0.) THEN
  74 : C    ... G > 0 => 1 racine réelle + 2 racines complexes conjuguées ...
  75 :          NRoot = 1
  76 :          XI = F / 2.d0
  77 :          XL = -XI + DSQRT(G)
  78 :          XM = -XI - DSQRT(G)
  79 :          X1 = DSIGN(1.D0,XL)*(DABS(XL)**D) +
  80 :      &        DSIGN(1.D0,XM)*(DABS(XM)**D) - XJ
  81 :       ENDIF
  82 : 
  83 :       RETURN
  84 :       END 
  85 :  
  86 :  
  87 :