Télécharger cubic.eso

Retour à la liste

Numérotation des lignes :

  1. C CUBIC SOURCE CB215821 16/04/21 21:16:13 8920
  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 ( ABS(C) .LE. EPSILO) THEN
  27. C ... On a alors affaire à une une équation quadratique ...
  28. CC2 = 1.d0
  29. CC1 = A
  30. CC0 = B
  31. CALL QUADRA(CC2, CC1, CC0, X2, X3, NRoot)
  32. NROOT=NROOT+1
  33. RETURN
  34. ELSE
  35. C ... Sinon, on élimine le terme x^2 en posant x' = x - A/3 ...
  36. C ... Nouvelle équation : x^3 + E x + F = 0 ...
  37. E = B - A*A / 3.d0
  38. F = A * (2.d0 * A * A - 9.d0 * B) / 27.d0 + C
  39. C ... G est son déterminant ...
  40. G = E*E*E / 27.d0 + F * F / 4.d0
  41. ENDIF
  42.  
  43. C ... Si G = 0 il y a 3 racines réelles, dont une double ...
  44. IF ( ABS(G) .LT. EPSILO) THEN
  45. NRoot = 3
  46. XI = F / 2.d0
  47. ALPHA = DSIGN(1.D0,-XI) * (ABS(XI)**D)
  48. X1 = 2.d0 * ALPHA - XJ
  49. X2 = -1.d0 * ALPHA - XJ
  50. X3 = X2
  51. RETURN
  52. ENDIF
  53.  
  54. C ... Sinon (G != 0) il y a 2 possibilités ...
  55. IF(DSIGN(1.D0,G). LT. 0.) THEN
  56. C ... G < 0 => 3 racines réelles ...
  57. NRoot = 3
  58. H = 2.d0 * SQRT(-E / 3.d0)
  59. B = -.5d0 * F / SQRT(-E * E * E / 27.d0)
  60. A = SQRT(1.d0 - B * B)
  61. IF( ABS(A).LT. EPSILO) THEN
  62. C = (1.d0 - DSIGN(1.D0,B)) * PIDIV2
  63. ELSE
  64. C = (2.d0 - DSIGN(1.D0,A)) * PIDIV2
  65. IF( ABS(B).GT. EPSILO) THEN
  66. C = C + DSIGN(1.D0,A)*DSIGN(1.D0,B)*
  67. & (ABS(ATAN(A / B)) - PIDIV2)
  68. ENDIF
  69. ENDIF
  70. X1 = H * COS(C / 3.d0) - XJ
  71. X2 = H * COS(C / 3.d0 + PI2 / 3.d0) - XJ
  72. X3 = H * COS(C / 3.d0 + 2.d0 * PI2 / 3.d0) - XJ
  73. RETURN
  74. ELSE IF ( DSIGN(1.D0,G) .GT. 0.) THEN
  75. C ... G > 0 => 1 racine réelle + 2 racines complexes conjuguées ...
  76. NRoot = 1
  77. XI = F / 2.d0
  78. XL = -XI + SQRT(G)
  79. XM = -XI - SQRT(G)
  80. X1 = DSIGN(1.D0,XL)*(ABS(XL)**D) +
  81. & DSIGN(1.D0,XM)*(ABS(XM)**D) - XJ
  82. ENDIF
  83.  
  84. RETURN
  85. END
  86.  
  87.  
  88.  
  89.  
  90.  
  91.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales