]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |