Télécharger quarti.eso

Retour à la liste

Numérotation des lignes :

  1. C QUARTI SOURCE CHAT 05/01/13 02:41:18 5004
  2. SUBROUTINE QUARTI (CA4, CA3, CA2, CA1, CA0, Q1, Q2, Q3,Q4,NRoot)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8(A-H,O-Z)
  5. DIMENSION QQ(4)
  6. D1 = CA3 / CA4
  7. D2 = CA2 / CA4
  8. D3 = CA1 / CA4
  9. D4 = CA0 / CA4
  10. D5 = -.25d0* D1
  11. A = D5 * D5
  12. B = A * D5
  13. C = B * D5
  14. P = 6.d0 * A + 3.d0 * D1 * D5 + D2
  15. Q = 4.d0 * B + 3.d0 * D1 * A + 2.d0 * D2 * D5 + D3
  16. R = C + D1 * B + D2 * A + D3 * D5 + D4
  17. IF( ABS(Q) . LE. 1.e-40) THEN
  18. NROOT=0
  19. A = P * P - 4.d0 * R
  20. IF( A . GE . 0.D0) THEN
  21. Y1= (-P + SQRT( A))/2.d0
  22. Y2= (-P - SQRT( A))/2.D0
  23. IF( Y1. GE. 0.D0) THEN
  24. NROOT = 2
  25. QQ(1)= SQRT(Y1)
  26. QQ(2) = -QQ(1)
  27. ENDIF
  28. IF( Y2.GE.0.D0) THEN
  29. QQ(NROOT+1)= SQRT(Y2)
  30. QQ(NROOT+2)=-QQ(NROOT+1)
  31. NROOT=NROOT+2
  32. ENDIF
  33. ENDIF
  34. ELSEIF( ABS( R) . LE. 1.e-40) THEN
  35. CALL CUBIC( 1.d0,0.d0, P,Q,QQ(1),QQ(2),QQ(3),NROOT)
  36. NROOT=NROOT+1
  37. QQ(NROOT) = 0.D0
  38. ELSE
  39. GO TO 12
  40. ENDIF
  41. Q1 = QQ(1) + D5
  42. Q2 = QQ(2) + D5
  43. Q3 = QQ(3) + D5
  44. Q4 = QQ(4) + D5
  45. RETURN
  46. *
  47. * cas general
  48. *
  49. 12 CONTINUE
  50. CB3 = 1.d0
  51. CB2 = 2.d0 * P
  52. CB1 = P * P - 4.d0 * R
  53. CB0 = -1.d0 * Q * Q
  54. CALL CUBIC(CB3, CB2, CB1, CB0, X1, X2, X3, NCUBIC)
  55. IF(NCUBIC.EQ.0) THEN
  56. RETURN
  57. ELSEIF(NCUBIC.EQ.1) THEN
  58. X4 = X1
  59. ELSEIF(NCUBIC.EQ.2) THEN
  60. X4 = X1
  61. IF( X2.GT. X4 ) X4 = X2
  62. ELSEIF(NCUBIC.EQ.3) THEN
  63. X4 = X1
  64. IF( X2.GT. X4) X4 = X2
  65. IF( X3.GT. X4) X4 = X3
  66. ENDIF
  67. X5 = SQRT(ABS(X4))
  68. CC2 = 1.d0
  69. CC1 = X5
  70. IF( X5 . NE . 0.d0) THEN
  71. CC0 = .5d0 * (P + X4 - Q / X5)
  72. ELSE
  73. CC0 = -1.d0
  74. ENDIF
  75. CALL QUADRA(CC2, CC1, CC0, R1, R2, NQUAD1)
  76. CC2 = 1.
  77. CC1 = -X5
  78. IF( X5 . NE . 0.d0) THEN
  79. CC0 = .5d0 * (P + X4 + Q / X5)
  80. ELSE
  81. CC0 = 1.d0
  82. ENDIF
  83. CALL QUADRA(CC2, CC1, CC0, R3, R4, NQUAD2)
  84. IF (NQUAD1 + NQUAD2 .EQ. 0) RETURN
  85. IF (NQUAD1 + NQUAD2.EQ. 4) THEN
  86. Q1 = R1 + D5
  87. Q2 = R2 + D5
  88. Q3 = R3 + D5
  89. Q4 = R4 + D5
  90. NRoot = 4
  91. RETURN
  92. ELSEIF( NQUAD1 .EQ. 0 )THEN
  93. Q1 = R3 + D5
  94. Q2 = R4 + D5
  95. NRoot = 2
  96. RETURN
  97. ELSE
  98. Q1 = R1 + D5
  99. Q2 = R2 + D5
  100. NRoot = 2
  101. ENDIF
  102. RETURN
  103. END
  104.  
  105.  
  106.  

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