Esope source files

   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
  34 :             CALL CUBIC( 1.d0,0.d0, P,Q,QQ(1),QQ(2),QQ(3),NROOT)
  35 :             NROOT=NROOT+1
  36 :             QQ(NROOT) = 0.D0
  37 :       ELSE 
  38 :             GO TO 12
  39 :       ENDIF
  40 :       Q1 = QQ(1) + D5
  41 :       Q2 = QQ(2) + D5 
  42 :       Q3 = QQ(3) + D5
  43 :       Q4 = QQ(4) + D5 
  44 :       RETURN
  45 : *
  46 : *   cas general
  47 : *
  48 :   12  CONTINUE
  49 :       CB3 = 1.d0
  50 :       CB2 = 2.d0 * P
  51 :       CB1 = P * P - 4.d0 * R
  52 :       CB0 = -1.d0 * Q * Q
  52 :       CALL CUBIC(CB3, CB2, CB1, CB0, X1, X2, X3, NCUBIC)
  53 :       IF(NCUBIC.EQ.0) THEN  
  54 :         RETURN
  55 :       ELSEIF(NCUBIC.EQ.1) THEN
  56 :          X4 = X1
  57 :       ELSEIF(NCUBIC.EQ.2) THEN
  58 :          X4 = X1
  59 :          IF( X2.GT. X4 ) X4 = X2
  60 :       ELSEIF(NCUBIC.EQ.3) THEN
  61 :          X4 = X1
  62 :          IF( X2.GT. X4) X4 = X2
  63 :          IF( X3.GT. X4)  X4 = X3
  64 :       ENDIF
  65 :       X5 = SQRT(ABS(X4))
  66 :       CC2 = 1.d0
  67 :       CC1 = X5
  68 :       IF( X5 . NE . 0.d0) THEN
  69 :             CC0 = .5d0 * (P + X4 - Q / X5)
  70 :        ELSE
  71 :           CC0 = -1.d0
  72 :        ENDIF
  72 :       CALL QUADRA(CC2, CC1, CC0, R1, R2, NQUAD1)
  73 :       CC2 = 1.
  74 :       CC1 = -X5
  75 :       IF( X5 . NE . 0.d0) THEN
  76 :       CC0 = .5d0 * (P + X4 + Q / X5)
  77 :       ELSE
  78 :             CC0 = 1.d0
  79 :       ENDIF
  79 :       CALL QUADRA(CC2, CC1, CC0, R3, R4, NQUAD2)
  80 :       IF (NQUAD1 + NQUAD2 .EQ. 0) RETURN
  81 :       IF (NQUAD1 + NQUAD2.EQ. 4) THEN
  82 :        Q1 = R1 + D5
  83 :        Q2 = R2 + D5
  84 :        Q3 = R3 + D5
  85 :        Q4 = R4 + D5
  86 :        NRoot = 4
  87 :        RETURN
  88 :       ELSEIF( NQUAD1 .EQ. 0 )THEN
  89 :        Q1 = R3 + D5
  90 :        Q2 = R4 + D5
  91 :        NRoot = 2
  92 :        RETURN
  93 :       ELSE
  94 :        Q1 = R1 + D5
  95 :        Q2 = R2 + D5
  96 :        NRoot = 2
  97 :       ENDIF
  98 :       RETURN
  99 :       END
 100 :  
 101 :