MINICERN/mathlib/gen/c/rteq364.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:52  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       SUBROUTINE RRTEQ3(R,S,T,X,D)
  12 #endif
  13 #if defined(CERNLIB_DOUBLE)
  14       SUBROUTINE DRTEQ3(R,S,T,X,D)
  15 #endif
  16 #include "gen/imp64.inc"
  17 #include "gen/defc64.inc"
  18      +  I,Z(0:2)
  19 #include "gen/def128.inc"
  20      + ZQ1,QQ,PP,Q1,Q2,Q3
  21       DIMENSION X(*),Y(0:2)
  22
  23       PARAMETER(EPS = 1D-6, DELTA = 1D-15)
  24       PARAMETER(I = (0,1), ZD1 = 1, ZQ1 = 1)
  25       PARAMETER(R1 = 2*ZD1/27, R2 = ZD1/2, R3 = ZD1/3)
  26       PARAMETER(W3 = 1.73205 08075 68877 294D0, R4 = W3/2)
  27       PARAMETER(Q1 = 2*ZQ1/27, Q2 = ZQ1/2, Q3 = ZQ1/3)
  28
  29       IF(S .EQ. 0 .AND. T .EQ. 0) THEN
  30        X(1)=-R
  31        X(2)=0
  32        X(3)=0
  33        D=0
  34        RETURN
  35       ENDIF
  36       P=S-R3*R**2
  37       Q=(R1*R**2-R3*S)*R+T
  38       D=(R2*Q)**2+(R3*P)**3
  39       IF(ABS(D) .LE. EPS) THEN
  40        PP=S-Q3*R**2
  41        QQ=(Q1*R**2-Q3*S)*R+T
  42        D=(Q2*QQ)**2+(Q3*PP)**3
  43        P=PP
  44        Q=QQ
  45       ENDIF
  46       H=R3*R
  47       H1=R2*Q
  48       IF(D .GT. DELTA) THEN
  49        H2=SQRT(D)
  50        U0=-H1+H2
  51        V0=-H1-H2
  52        U=SIGN(ABS(U0)**R3,U0)
  53        V=SIGN(ABS(V0)**R3,V0)
  54        X(1)=U+V-H
  55        X(2)=-R2*(U+V)-H
  56        X(3)=R4*ABS(U-V)
  57        IF(ABS(U0) .LE. EPS .OR. ABS(V0) .LE. EPS) THEN
  58         Y(0)=X(1)
  59         DO 1 K = 0,1
  60     1   Y(K+1)=Y(K)-(((Y(K)+R)*Y(K)+S)*Y(K)+T)/((3*Y(K)+2*R)*Y(K)+S)
  61         X(1)=Y(2)
  62         Z(0)=X(2)+I*X(3)
  63         DO 2 K = 0,1
  64     2   Z(K+1)=Z(K)-(((Z(K)+R)*Z(K)+S)*Z(K)+T)/((3*Z(K)+2*R)*Z(K)+S)
  65         X(2)=Z(2)
  66         X(3)=-I*Z(2)
  67        ENDIF
  68       ELSEIF(ABS(D) .LE. DELTA) THEN
  69        D=0
  70        U=SIGN(ABS(H1)**R3,-H1)
  71        X(1)=U+U-H
  72        X(2)=-U-H
  73        X(3)=X(2)
  74        IF(ABS(H1) .LE. EPS) THEN
  75         Y(0)=X(1)
  76         DO 5 K = 0,1
  77         H1=(3*Y(K)+2*R)*Y(K)+S
  78         IF(ABS(H1) .GT. DELTA) THEN
  79          Y(K+1)=Y(K)-(((Y(K)+R)*Y(K)+S)*Y(K)+T)/H1
  80         ELSE
  81          X(1)=-R3*R
  82          X(2)=X(1)
  83          X(3)=X(1)
  84          RETURN
  85         ENDIF
  86     5   CONTINUE
  87         X(1)=Y(2)
  88         X(2)=-R2*(R+X(1))
  89         X(3)=X(2)
  90        ENDIF
  91       ELSE
  92        H3=SQRT(ABS(R3*P)**3)
  93        H2=R3*ACOS(-H1/H3)
  94        H1=H3**R3
  95        U=H1*COS(H2)
  96        V=W3*H1*SIN(H2)
  97        X(1)=U+U-H
  98        X(2)=-U-V-H
  99        X(3)=-U+V-H
 100        IF(H3 .LE. EPS .OR. X(1) .LE. EPS .OR. X(2) .LE .EPS .OR.
 101      1                     X(3) .LE. EPS) THEN
 102         DO 3 J = 1,3
 103         Y(0)=X(J)
 104         DO 4 K = 0,1
 105     4   Y(K+1)=Y(K)-(((Y(K)+R)*Y(K)+S)*Y(K)+T)/((3*Y(K)+2*R)*Y(K)+S)
 106     3   X(J)=Y(2)
 107        ENDIF
 108       ENDIF
 109       RETURN
 110       END