Télécharger la plaquette Cast3M


fr en

Télécharger quarti.eso

Retour à la liste

Numérotation des lignes : 

quarti

  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