[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / reli264.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:11  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
#if defined(CERNLIB_DOUBLE)
      FUNCTION DELI2(X,AKP,A,B,MODE)
C
#include "gen/imp64.inc"
C
      CHARACTER*(*) NAME
      PARAMETER(NAME='RELI2/DELI2')
#endif
#if !defined(CERNLIB_DOUBLE)
      FUNCTION RELI2(X,AKP,A,B,MODE)
C
      CHARACTER*(*) NAME
      PARAMETER(NAME='RELI2')
#endif
C
C     Translation of Algol procedure el2(x,kc,a,b) in
C      R. BULIRSCH Numerical Calculation of Elliptic Integrals and
C      Elliptic Functions,  Numer. Math. 7 (1965) 78-90,
C      extended for (k sub c)**2 < 0  (MODE = -1)
C
      PARAMETER (ID = 16)
      PARAMETER (PI = 3.14159 26535 89793 24D0)
      PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2)
      PARAMETER (CA = Z10**(-ID/2), CB = Z10**(-(ID+2)))

      CHARACTER*80 ERRTXT

      IF(X .EQ. 0 .OR. AKP .EQ. 0) THEN
#if defined(CERNLIB_DOUBLE)
       H=B*DASINH(X)+(A-B)*X/SQRT(1+X**2)
#endif
#if !defined(CERNLIB_DOUBLE)
       H=B*ASINH(X)+(A-B)*X/SQRT(1+X**2)
#endif
      ELSE
       IF(MODE .EQ. 1) THEN
        XX=X
        BKP=AKP
        AA=A
        BB=B
       ELSEIF(MODE .EQ. -1) THEN
        W=1-(AKP*X)**2
        IF(W .GT. 0) THEN
         W1=1/SQRT(1+AKP**2)
         XX=X/(W1*SQRT(W))
         BKP=AKP*W1
         AA=A*W1
         BB=(B+A*AKP**2)*W1**3
        ELSE
         H=0
         WRITE(ERRTXT,101) X,AKP
         CALL MTLPRT(NAME,'C346.1',ERRTXT)
         GO TO 9
        ENDIF
       ELSE
        H=0
        WRITE(ERRTXT,102) MODE
        CALL MTLPRT(NAME,'C346.2',ERRTXT)
        GO TO 9
       ENDIF
       C=XX**2
       D=1+C
       PP=SQRT((1+BKP**2*C)/D)
       D=XX/D
       C=D/(2*PP)
       Z=AA-BB
       XI=AA
       AA=HF*(AA+BB)
       Y=ABS(1/XX)
       F=0
       L=0
       XM=1
       YKP=ABS(BKP)
    1  BB=XI*YKP+BB
       E=XM*YKP
       G=E/PP
       D=F*G+D
       F=C
       XI=AA
       PP=G+PP
       C=HF*(D/PP+C)
       G=XM
       XM=YKP+XM
       AA=HF*(BB/XM+AA)
       Y=Y-E/Y
       IF(Y .EQ. 0) Y=SQRT(E)*CB
       IF(ABS(G-YKP) .GT. CA*G) THEN
        YKP=2*SQRT(E)
        L=2*L
        IF(Y .LT. 0) L=L+1
        GO TO 1
       ENDIF
       IF(Y .LT. 0) L=L+1
       E=(ATAN(XM/Y)+PI*L)*AA/XM
       IF(XX .LT. 0) E=-E
       H=E+C*Z
      ENDIF
#if defined(CERNLIB_DOUBLE)
    9 DELI2=H
#endif
#if !defined(CERNLIB_DOUBLE)
    9 RELI2=H
#endif
      RETURN
  101 FORMAT('X = ',1P,D15.8,' AKP = ',D15.8,' ILLEGAL',
     1       '[(AKP * X)**2 >= 1]')
  102 FORMAT('MODE = ',I5,' ILLEGAL')
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1996/04/01 15:02:11 mclareni
	6	* Mathlib gen
	7	*
	8	*
	9	#include "gen/pilot.h"
	10	#if defined(CERNLIB_DOUBLE)
	11	FUNCTION DELI2(X,AKP,A,B,MODE)
	12	C
	13	#include "gen/imp64.inc"
	14	C
	15	CHARACTER() NAME
	16	PARAMETER(NAME='RELI2/DELI2')
	17	#endif
	18	#if !defined(CERNLIB_DOUBLE)
	19	FUNCTION RELI2(X,AKP,A,B,MODE)
	20	C
	21	CHARACTER() NAME
	22	PARAMETER(NAME='RELI2')
	23	#endif
	24	C
	25	C Translation of Algol procedure el2(x,kc,a,b) in
	26	C R. BULIRSCH Numerical Calculation of Elliptic Integrals and
	27	C Elliptic Functions, Numer. Math. 7 (1965) 78-90,
	28	C extended for (k sub c)**2 < 0 (MODE = -1)
	29	C
	30	PARAMETER (ID = 16)
	31	PARAMETER (PI = 3.14159 26535 89793 24D0)
	32	PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2)
	33	PARAMETER (CA = Z10(-ID/2), CB = Z10(-(ID+2)))
	34
	35	CHARACTER*80 ERRTXT
	36
	37	IF(X .EQ. 0 .OR. AKP .EQ. 0) THEN
	38	#if defined(CERNLIB_DOUBLE)
	39	H=BDASINH(X)+(A-B)X/SQRT(1+X**2)
	40	#endif
	41	#if !defined(CERNLIB_DOUBLE)
	42	H=BASINH(X)+(A-B)X/SQRT(1+X**2)
	43	#endif
	44	ELSE
	45	IF(MODE .EQ. 1) THEN
	46	XX=X
	47	BKP=AKP
	48	AA=A
	49	BB=B
	50	ELSEIF(MODE .EQ. -1) THEN
	51	W=1-(AKPX)*2
	52	IF(W .GT. 0) THEN
	53	W1=1/SQRT(1+AKP**2)
	54	XX=X/(W1*SQRT(W))
	55	BKP=AKP*W1
	56	AA=A*W1
	57	BB=(B+AAKP2)W1**3
	58	ELSE
	59	H=0
	60	WRITE(ERRTXT,101) X,AKP
	61	CALL MTLPRT(NAME,'C346.1',ERRTXT)
	62	GO TO 9
	63	ENDIF
	64	ELSE
65	H=0
66	WRITE(ERRTXT,102) MODE
67	CALL MTLPRT(NAME,'C346.2',ERRTXT)
68	GO TO 9
69	ENDIF
70	C=XX**2
71	D=1+C
72	PP=SQRT((1+BKP*2C)/D)
73	D=XX/D
74	C=D/(2*PP)
75	Z=AA-BB
76	XI=AA
77	AA=HF*(AA+BB)
78	Y=ABS(1/XX)
79	F=0
80	L=0
81	XM=1
82	YKP=ABS(BKP)
83	1 BB=XI*YKP+BB
84	E=XM*YKP
85	G=E/PP
86	D=F*G+D
87	F=C
88	XI=AA
89	PP=G+PP
90	C=HF*(D/PP+C)
91	G=XM
92	XM=YKP+XM
93	AA=HF*(BB/XM+AA)
94	Y=Y-E/Y
95	IF(Y .EQ. 0) Y=SQRT(E)*CB
96	IF(ABS(G-YKP) .GT. CA*G) THEN
97	YKP=2*SQRT(E)
98	L=2*L
99	IF(Y .LT. 0) L=L+1
100	GO TO 1
101	ENDIF
102	IF(Y .LT. 0) L=L+1
103	E=(ATAN(XM/Y)+PIL)AA/XM
104	IF(XX .LT. 0) E=-E
105	H=E+C*Z
106	ENDIF
107	#if defined(CERNLIB_DOUBLE)
108	9 DELI2=H
109	#endif
110	#if !defined(CERNLIB_DOUBLE)
111	9 RELI2=H
112	#endif
113	RETURN
114	101 FORMAT('X = ',1P,D15.8,' AKP = ',D15.8,' ILLEGAL',
115	1 '[(AKP * X)**2 >= 1]')
116	102 FORMAT('MODE = ',I5,' ILLEGAL')
117	END