IMPLEMENTATION MODULE SpezFunkt3;
(*------------------------------------------------------------------------*)
(* Stellt verschiedene mathematische Funktionen zur Verf"ugung *)
(* This library provides some special mathematical functions. *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 14.05.16, MRi: Prozeduren Bessel, Y01, BessJ0, H01 aus LMathLib her- *)
(* ausgenommen und Modul SpezFunkt3 erstellt *)
(* 29.05.16, MRi: Bessel in BesselJN, BessJ0 in BesselJ0, Y01 in *)
(* BesselY0 umbenannt. *)
(* 30.05.16, MRi: H01 in Hankel01 umbenannt, Fehler korrigiert *)
(* 11.09.17, MRi: ZeroBessel{0|1} eingefuegt *)
(* 09.05.18, MRi: Erstellen der ersten Version von DBesselJ0 *)
(* 09.05.18, MRi: Korrektur in BesselJ0 *)
(* 10.05.18, MRi: Erstellen der ersten neuen Versionen von BesselJ0 und *)
(* BesselJ1 *)
(* 12.05.18, MRi: Korrektur in BesselJn (Rekursion von J0/J1 fuer grosse *)
(* x, Korrektur in BesselJN (start) *)
(* 13.05.18, MRi: Korrektur in BesselY0, neue Version von BesselY1, die *)
(* ersten beiden Nullstellen in ZeroBesselJ{0|1} hart *)
(* kodiert. *)
(* 14.05.18, MRi: Hankel01 erweitert *)
(* 02.02.19, MRi: Kosmetik, I bzw IM durch "imag" ersetzt *)
(*------------------------------------------------------------------------*)
(* - Procedures BesselJ0, BesselJ1, BesselJN, BesselY01 are Modula-2 *)
(* adoptions of routines found *)
(* NUMAL Numerical Algol library (Stichting CWI) *)
(* - Procedure DBesselJ0 and DBesselY1 are derived form the *)
(* CERN 2006 library *)
(* - Procedure ZeroBessel{0|1} are derived from Pascal mathlib library *)
(* PMATH is provided by Jean Debord, l'Universite de Limoges, France *)
(* - Procedure and initialisation code in local module CHank is derived *)
(* from Fortran subroutines CHAE01 and CINITAL from Derek Nielson. *)
(*------------------------------------------------------------------------*)
(* Testroutinen in TstSpFun3.mod gegen SpezFunkt.f (Fortran) *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: SpezFunkt3.mod,v 1.7 2019/02/02 12:33:14 mriedl Exp mriedl $ *)
FROM LongMath IMPORT pi,sqrt,sin,cos,ln;
IMPORT LongComplexMath;
FROM LMathLib IMPORT MachEps;
FROM Errors IMPORT ErrOut;
IMPORT LMathLib;
IMPORT TIO;
CONST CEXP = LongComplexMath.exp;
CABS = LongComplexMath.abs;
scalarMult = LongComplexMath.scalarMult;
CONST zero = LongComplexMath.zero;
one = LongComplexMath.one;
(*
* PROCEDURE scalarMult(a : LONGREAL;
* x : LONGCOMPLEX) : LONGCOMPLEX;
*
* BEGIN
* RETURN CMPLX(a*RE(x),a*IM(x));
* END scalarMult;
*)
CONST debug = (1 = 0);
PROCEDURE DBesselJ0(x : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Source : CERN library 2006 *)
(* Fortran bearbeitet : MRi, 07.05.18 *)
(*----------------------------------------------------------------*)
CONST R32 = 1.0 / 32.0;
PI1 = 2.0 / pi;
PI4 = pi / 4.0;
imag = CMPLX(0.0,1.0);
TYPE RVek17 = ARRAY [1..17] OF LONGREAL;
CONST A = RVek17{
0.157727971474890120, -0.008723442352852221,
0.265178613203336810, -0.370094993872649779,
0.158067102332097261, -0.034893769411408885,
0.004819180069467605, -0.000460626166206275,
0.000032460328821005, -0.000001761946907762,
0.000000076081635924, -0.000000002679253531,
0.000000000078486963, -0.000000000001943835,
0.000000000000041253, -0.000000000000000759,
0.000000000000000012
};
TYPE CVek20 = ARRAY [1..20] OF LONGCOMPLEX;
CONST C = CVek20{
CMPLX( 0.998988089858965153, -0.012331520578544144),
CMPLX(-0.001338428549971856, -0.012249496281259475),
CMPLX(-0.000318789878061893, 0.000096494184993423),
CMPLX( 0.000008511232210657, 0.000013655570490357),
CMPLX( 0.000000691542349139, -0.000000851806644426),
CMPLX(-0.000000090770101537, -0.000000027244053414),
CMPLX( 0.000000001454928079, 0.000000009646421338),
CMPLX( 0.000000000926762487, -0.000000000683347518),
CMPLX(-0.000000000139166198, -0.000000000060627380),
CMPLX( 0.000000000003237975, 0.000000000021695716),
CMPLX( 0.000000000002535357, -0.000000000002304899),
CMPLX(-0.000000000000559090, -0.000000000000122554),
CMPLX( 0.000000000000041919, 0.000000000000092314),
CMPLX( 0.000000000000008733, -0.000000000000016778),
CMPLX(-0.000000000000003619, 0.000000000000000754),
CMPLX( 0.000000000000000594, 0.000000000000000462),
CMPLX(-0.000000000000000010, -0.000000000000000159),
CMPLX(-0.000000000000000024, 0.000000000000000025),
CMPLX( 0.000000000000000008, 0.000000000000000000),
CMPLX(-0.000000000000000001, -0.000000000000000001)
};
VAR alfa,b0,b1,b2,h,p,r,v : LONGREAL;
cb0,cb1,cb2,ctmp : LONGCOMPLEX;
i : CARDINAL;
BEGIN
v := ABS(x);
IF (v < 8.0) THEN
h := R32*v*v - 1.0;
alfa := h + h;
b1 := 0.0; b2 := 0.0;
FOR i:=17 TO 1 BY -1 DO
b0 := A[i] + alfa*b1 - b2;
b2 := b1;
b1 := b0;
END;
p := b0 - h*b2;
ELSE
r := 1.0 / v;
h := 10.0*r - 1.0;
alfa := h + h;
cb1 := zero;
cb2 := zero;
FOR i:=20 TO 1 BY -1 DO
cb0 := C[i] + scalarMult(alfa,cb1) - cb2;
cb2 := cb1;
cb1 := cb0;
END;
ctmp := CEXP(imag*CMPLX((v - PI4),0.0));
cb0 := scalarMult(sqrt(PI1*r),ctmp)*(cb0 - scalarMult(h,cb2)); (*red?*)
p := RE(cb0);
END;
RETURN p;
END DBesselJ0;
PROCEDURE besspq0( x : LONGREAL;
VAR p0,q0 : LONGREAL);
(*---------------------------------------------------------------*)
(* This procedure is an auxiliary procedure for the computation *)
(* of the ordinary bessel functions of order zero for large *)
(* values of their argument *)
(* *)
(* x : the argument, this argument should satisfy x > 0 *)
(* p0 : on exit p0 corresponds with the function p(x,0) defined *)
(* in [1], formulas 9.2.5 and 9.2.6 *)
(* q0 : on exit q0 corresponds with the function q(x,0) defined *)
(* in [1], formulas 9.2.5 and 9.2.6 *)
(* *)
(* [1] Abramowitz, Milton, Stegun, Irene, *)
(* "Handbook of mathematical functions with formulas, graphs *)
(* and mathematical tables", Appl. Math. 55, (1964). *)
(*---------------------------------------------------------------*)
TYPE Vek10 = ARRAY [1..10] OF LONGREAL;
Vek11 = ARRAY [1..11] OF LONGREAL;
CONST ar1 = Vek11{
-0.10012E-15, +0.67481E-15,
-0.506903E-14, +0.4326596E-13,
-0.43045789E-12, +0.516826239E-11,
-0.7864091377E-10, +0.163064646352E-08,
-0.5170594537606E-07, +0.30751847875195E-05,
-0.536522046813212E-03
};
CONST ar2 = Vek10{
-0.60999E-15, +0.425523E-14,
-0.3336328E-13, +0.30061451E-12,
-0.320674742E-11, +0.4220121905E-10,
-0.72719159369E-09, +0.1797245724797E-07,
-0.74144984110606E-06, +0.683851994261165E-04
};
VAR b,cosx,sinx,j0x,y0 : LONGREAL;
x2,b0,b1,b2,y : LONGREAL;
i : CARDINAL;
BEGIN
IF (x < 8.0) THEN
IF debug THEN TIO.WrStr(" (x < 8)"); END;
b:= sqrt(x)*1.25331413731550;
BesselY01(x,y0,j0x); j0x := BesselJ0(x);
x:= x - 0.785398163397448; cosx:= cos(x); sinx:= sin(x);
p0:= b*(y0*sinx + j0x*cosx);
q0:= b*(y0*cosx - j0x*sinx);
ELSE
IF debug THEN TIO.WrStr(" (x >= 8)"); END;
y:= 8.0/x; x:= 2.0*y*y - 1.0; x2 := x+x;
b1:=0.0; b2:=0.0;
FOR i:=1 TO 11 DO
b0 := x2*b1 - b2 + ar1[i];
b2 := b1; b1 := b0;
END;
p0:= x*b1 - b2 + 0.99946034934752;
(* computation of q0 *)
b1:=0.0; b2:=0.0;
FOR i:=1 TO 10 DO
b0 := x2*b1 - b2 + ar2[i];
b2 := b1; b1 := b0;
END;
q0:=(x*b1 - b2 -0.015555854605337)*y
END;
END besspq0;
PROCEDURE BesselJ0(x : LONGREAL) : LONGREAL;
TYPE Vek15 = ARRAY [1..15] OF LONGREAL;
CONST ar = Vek15{
-0.75885E-15, +0.4125321E-13,
-0.194383469E-11, +0.7848696314E-10,
-0.267925353056E-08, +0.7608163592419E-07,
-0.176194690776215E-05, +0.324603288210051E-04,
-0.46062616620628E-03, +0.48191800694676E-02,
-0.34893769411409E-01, +0.158067102332097,
-0.37009499387265E-00, +0.265178613203337,
-0.872344235285222E-02
};
CONST sqrt2oPi = 0.797884560802865355879892119868764;
VAR z,z2,b0,b1,b2,p0,q0 : LONGREAL;
cosx,sinx,ss,cc : LONGREAL;
i : CARDINAL;
BEGIN
IF debug THEN TIO.WrStr("Aufruf von bessj0"); TIO.WrLn; END;
IF (x = 0.0) THEN
RETURN 1.0;
ELSIF (ABS(x) < 8.0) THEN
x:= x/8.0; z:= 2.0*x*x - 1.0; z2 := z+z;
b1:=0.0; b2:=0.0;
FOR i:=1 TO 15 DO
b0 := z2*b1 - b2 + ar[i];
b2 := b1; b1 := b0;
END;
RETURN z*b1 - b2 + 0.15772797147489;
ELSE
x := ABS(x); (* Codefragment korrigiert, nicht original NumaL Algol *)
sinx := sin(x); cosx:=cos(x);
ss := sinx - cosx; cc := sinx + cosx;
besspq0(x,p0,q0);
RETURN (1.0/sqrt(pi*x))*(p0*cc - q0*ss);
END;
END BesselJ0;
PROCEDURE besspq1( x : LONGREAL;
VAR p1,q1 : LONGREAL);
(*---------------------------------------------------------------*)
(* This procedure is an auxiliary procedure for the computation *)
(* of the ordinary bessel functions of order one for large *)
(* values of their argument. *)
(* *)
(* x : the argument, this argument should satisfy x > 0 *)
(* p1 : on exit p1 corresponds with the function p(x,1) defined *)
(* in [1,formulas 9.2.5 and 9.2.6] *)
(* q1 : on exit q1 corresponds with the function q(x,1) defined *)
(* in [1,formulas 9.2.5 and 9.2.6] *)
(*---------------------------------------------------------------*)
TYPE Vek11 = ARRAY [1..11] OF LONGREAL;
CONST ar1 = Vek11{
+0.10668E-15,
-0.72212E-15, +0.545267E-14,
-0.4684224E-13, +0.46991955E-12,
-0.570486364E-11, +0.881689866E-10,
-0.187189074911E-08, +0.6177633960644E-07,
-0.39872843004889E-05, +0.89898983308594E-03
};
CONST ar2 = Vek11{
-0.10269E-15, +0.65083E-15,
-0.456125E-14, +0.3596777E-13,
-0.32643157E-12, +0.351521879E-11,
-0.4686363688E-10, +0.82291933277E-09,
-0.2095978138408E-07, +0.91386152579555E-06,
-0.96277235491571E-04
};
VAR b,cosx,sinx,j1x,y1 : LONGREAL;
x2,b0,b1,b2,y : LONGREAL;
i : CARDINAL;
BEGIN
IF debug THEN TIO.WrStr("Aufruf von besspq1"); TIO.WrLn; END;
IF (x < 8.0) THEN
b := sqrt(x)*1.25331413731550; (* Wie soll hier b gesetzt werden ??? *)
BesselY01(x,j1x,y1);
j1x := BesselJ1(x);
x:= x - 0.785398163397448; cosx:= cos(x); sinx:= sin(x);
p1:= b*(j1x*sinx - y1*cosx);
q1:= b*(j1x*cosx + y1*sinx)
ELSE
y:= 8.0 / x; x:= 2.0*y*y - 1.0; x2 := x+x;
(* computation of p1; *)
b1:=0.0; b2:=0.0;
FOR i:=1 TO 11 DO
b0:= b1*x2 - b2 + ar1[i];
b2:= b1; b1:= b0
END;
p1:= x*b1 - b2 + 1.0009030408600137;
(* computation of q1; *)
b1:=0.0; b2:=0.0;
FOR i:=1 TO 11 DO
b0:= x2 * b1 - b2 + ar2[i];
b2:= b1; b1:= b0
END;
q1:=(x*b1 - b2 + 0.46777787069535E-01)*y
END
END besspq1;
PROCEDURE BesselJ1(x : LONGREAL) : LONGREAL;
TYPE Vek15 = ARRAY [1..15] OF LONGREAL;
CONST ar = Vek15{
-0.19554E-15, +0.1138572E-13,
-0.57774042E-12, +0.2528123664E-10,
-0.94242129816E-09, +0.2949707007278E-07,
-0.76175878054003E-06, +0.158870192399321E-04,
-0.260444389348581E-03, +0.324027018268386E-02,
-0.291755248061542E-01, +0.177709117239728E-00,
-0.661443934134543E-00, +0.128799409885768E+01,
-0.119180116054122E+01
};
VAR z,z2,b0,b1,b2,sgnx : LONGREAL;
c,cosx,sinx,p1,q1 : LONGREAL;
i : CARDINAL;
BEGIN
IF debug THEN TIO.WrStr("Aufruf von bessj1"); TIO.WrLn; END;
IF (x = 0.0) THEN
RETURN 0.0;
ELSIF (ABS(x) < 8.0) THEN
x:= x/8.0; z:= 2.0*x*x - 1.0; z2:= z+z;
(* computation of j1 *)
b1:=0.0; b2:=0.0;
FOR i:=1 TO 15 DO
b0 := z2*b1 - b2 + ar[i];
b2 := b1; b1 := b0;
END;
RETURN x*(z*b1 - b2 + 0.648358770605265);
ELSE
IF (x < 0.0) THEN sgnx:=-1.0; ELSE sgnx:=1.0; END;
x := ABS(x);
c:= 0.797884560802865 / sqrt(x);
cosx:= cos(x - 0.706858347057703E+01);
sinx:= sin(x - 0.706858347057703E+01);
besspq1(x,p1,q1);
RETURN sgnx*c*(p1*sinx + q1*cosx)
END
END BesselJ1;
PROCEDURE BesselJn(x : LONGREAL;
n : CARDINAL) : LONGREAL;
VAR xh2ipn,xh,Z,Sum,Fakt,Faktipn : LONGREAL;
res,zdx,BesJm1,BesJ,BesJp1 : LONGREAL;
i,ipn : CARDINAL;
BEGIN
IF (x = 0.0) THEN
res := 0.0;
ELSIF (ABS(x) > 0.95*LFLOAT(n)) THEN
IF (n = 0) THEN
res := BesselJ0(x);
ELSIF (n = 1) THEN
res := BesselJ1(x);
ELSE
(* Aufwaertsrekusion von J0 und J1 *)
x := ABS(x);
zdx := 2.0 / x;
BesJm1 := BesselJ0(x);
BesJ := BesselJ1(x);
FOR i:=1 TO n-1 DO
BesJp1 := LFLOAT(i)*zdx*BesJ - BesJm1;
BesJm1 := BesJ;
BesJ := BesJp1;
END;
res := BesJ;
END;
ELSE
(* Reihenentwicklung *)
xh:=0.5*x;
xh2ipn:=LMathLib.CardPot(xh,n);
xh:=xh*xh;
Sum:=0.0; i:=0; ipn:=n; (* ipn=i+n *)
Fakt:=1.0; Faktipn:=LMathLib.fact(n);
LOOP
Z:=Sum;
IF ODD(i) THEN
Sum:=Sum - xh2ipn/(Fakt*Faktipn);
ELSE
Sum:=Sum + xh2ipn/(Fakt*Faktipn);
END;
IF (Z = Sum) THEN EXIT END;
INC(i); INC(ipn);
Fakt:=Fakt*LFLOAT(i);
Faktipn:=Faktipn*LFLOAT(ipn);
xh2ipn:=xh2ipn*xh;
END; (* LOOP *)
res := Sum;
END;
RETURN res;
END BesselJn;
PROCEDURE start(x : LONGREAL;
n : INTEGER;
t : INTEGER) : INTEGER;
(*----------------------------------------------------------------*)
(* generate a starting value for the miller algorithm for *)
(* computing an array of bessel functions *)
(* *)
(* x : the argument of the bessel functions, x > 0 *)
(* n : the number of bessel functions to be computed, n >= 0 *)
(* t : the type of bessel function in question, *)
(* t = 0 corresponds to ordinary bessel functions, *)
(* t = 1 corresponds to modified bessel functions. *)
(*----------------------------------------------------------------*)
VAR p,q,r,y : LONGREAL;
s : INTEGER;
BEGIN
s := 2*t-1;
p := 36.0/x - LFLOAT(t); r:= LFLOAT(n)/x;
IF (r > 1.0) OR (t = 1) THEN
q := sqrt(r*r + LFLOAT(s)); r:=r*ln(q + r) - q;
ELSE
r:=0.0;
END;
q:= 18.0/x + r;
IF (p > q) THEN r:=p ELSE r:=q; END;
p := sqrt(2.0*(LFLOAT(t) + r));
p := x*((1.0 + r) + p) / (1.0 + p);
q:=0.0;
WHILE (p > q) OR (p < q-1.0) DO
y := p;
p := p/x;
q := sqrt(p*p + LFLOAT(s));
p := x*(r + q) / ln(p + q);
q := y;
END;
IF (t = 1) THEN
RETURN + LMathLib.entier( p+1.0);
ELSE
RETURN - 2*LMathLib.entier(-p/2.0);
END;
END start;
PROCEDURE BesselJN( x : LONGREAL;
n : INTEGER;
VAR Jn : ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* Generates an array of ordinary bessel functions of the first *)
(* kind of order l (l = 0,...,n) with argument x. *)
(* *)
(* x : the argument of the bessel functions *)
(* n : the upper bound of the indices of array j, n >= 0 *)
(* Jn : array j[0:n], on exit j[l] is the ordinary bessel *)
(* function of the first kind of order l and argument x. *)
(* *)
(* [3] Gautschi, W., "Computational aspects of three term *)
(* recurrence relations", SIAM review, 9, pp 24-82 (1967) *)
(*----------------------------------------------------------------*)
VAR x2,r,s : LONGREAL;
l,m,nu : INTEGER;
signx : BOOLEAN;
BEGIN
IF (x = 0.0) THEN
Jn[0]:= 1.0;
FOR l:=n TO 1 BY -1 DO Jn[l]:= 0.0; END;
ELSE
signx := (x < 0.0);
x:= ABS(x);
r:=0.0; s:=0.0; x2:= 2.0/x; l:=0;
nu := start(x,n,0);
FOR m:=nu TO 1 BY -1 DO
r := 1.0 / (x2*LFLOAT(m) - r);
l := 2 - l; s := r*(LFLOAT(l) + s);
IF (m <= n) THEN Jn[m]:= r; END;
END;
r := 1.0 / (1.0 + s);
Jn[0] := r;
FOR m:=1 TO n DO r:=r*Jn[m]; Jn[m]:= r; END;
IF signx THEN
FOR m:=1 TO n BY 2 DO Jn[m] := -Jn[m]; END;
END;
END;
END BesselJN;
PROCEDURE ZeroBesselJ0 (n : INTEGER) : LONGREAL;
VAR x,y,Zero : LONGREAL;
BEGIN
IF (n = 0) THEN
RETURN MAX(LONGREAL);
ELSIF (n = 1) THEN
RETURN 2.404825557695772768622;
ELSIF (n = 2) THEN
RETURN 5.520078110286310649597;
ELSE
x := pi*VAL(REAL,(4*n-1));
y := 1.0 / (x*x);
Zero := x/4.0*(1.0 + y*(2.0 + y*(-62.0/3.0 + y*(15116.0/15.0 +
y*(-12554474.0/105.0 + y*8368654292.0/315.0)))));
(* Newton step to refine the root *)
RETURN Zero + (BesselJ0(Zero) / BesselJ1(Zero));
END;
END ZeroBesselJ0;
PROCEDURE ZeroBesselJ1(n : CARDINAL) : LONGREAL;
VAR x,y,Zero : LONGREAL;
BEGIN
IF (n = 0) THEN
RETURN 0.0;
ELSIF (n = 1) THEN
RETURN 3.831705970207512315614;
ELSIF (n = 2) THEN
RETURN 7.015586669815618753537;
ELSE
x := pi*VAL(LONGREAL,(4*n+1));
y := 1.0 / (x*x);
Zero := x/4.0*(1.0 + y*(-6.0 + y*(6.0 + y*(-4716.0/5.0 +
y*(-3902418.0/35.0 - y*8952167324.0/35.0)))));
(* Newton step to refine the root *)
RETURN Zero - BesselJ1(Zero) / (BesselJ0(Zero) - BesselJ1(Zero)/Zero);
END;
END ZeroBesselJ1;
PROCEDURE BesselY0(x : LONGREAL) : LONGREAL;
CONST C = 0.577215664901532; (* Eulerkonstante. *)
VAR xh2,y,Z,Sum,Fakt,Phi : LONGREAL;
p0,q0 : LONGREAL;
i : CARDINAL;
BEGIN
IF (x < 8.0) THEN
xh2:=0.25*x*x; Sum:=-2.0*xh2;
i:=2; Fakt:=2.0;
Phi:=1.5; y:=xh2*xh2;
LOOP
Z:=Sum;
IF ODD(i) THEN
Sum:=Sum - 2.0*y*Phi/(Fakt*Fakt);
ELSE
Sum:=Sum + 2.0*y*Phi/(Fakt*Fakt);
END; (* IF *)
IF (Z = Sum) THEN EXIT END;
INC(i);
Fakt:=Fakt*LFLOAT(i);
y:=y*xh2;
Phi:=Phi + 1.0 / LFLOAT(i);
END; (* LOOP *)
RETURN (2.0*(C + ln(x/2.0))*BesselJ0(x) - Sum) / pi;
ELSE (* 13.05.18 *)
besspq0(x,p0,q0);
y := x - 0.706858347057703E01;
RETURN (0.797884560802865 / sqrt(x))*(p0*sin(y) + q0*cos(y));
END;
END BesselY0;
PROCEDURE DBesselY1(X : LONGREAL) : LONGREAL;
CONST imag = CMPLX(0.0,1.0);
CE = 0.577215664901532861;
PI1 = 2.0 / pi;
PI3 = 3.0*pi / 4.0;
VAR alfa,b0,b1,b2,h,p,r,v,y : LONGREAL;
i : CARDINAL;
cb0,cb1,cb2 : LONGCOMPLEX;
TYPE RVek17 = ARRAY [1..17] OF LONGREAL;
CVek20 = ARRAY [1..20] OF LONGCOMPLEX;
CONST D = RVek17{
0.648358770605264921, -1.191801160541216873,
1.287994098857677620, -0.661443934134543253,
0.177709117239728283, -0.029175524806154208,
0.003240270182683857, -0.000260444389348581,
0.000015887019239932, -0.000000761758780540,
0.000000029497070073, -0.000000000942421298,
0.000000000025281237, -0.000000000000577740,
0.000000000000011386, -0.000000000000000196,
0.000000000000000003
};
CONST E = RVek17{
-0.040172946544414076, -0.444447147630558063,
-0.022719244428417736, 0.206644541017490520,
-0.086671697056948524, 0.017636703003163134,
-0.002235619294485095, 0.000197062302701541,
-0.000012885853299241, 0.000000652847952359,
-0.000000026450737175, 0.000000000878030117,
-0.000000000024343279, 0.000000000000572612,
-0.000000000000011578, 0.000000000000000203,
-0.000000000000000003
};
CONST F = CVek20{
CMPLX( 1.001702234853820996, 0.037261715000537654),
CMPLX( 0.002255572846561180, 0.037145322479807690),
CMPLX( 0.000543216487508013, -0.000137263238201907),
CMPLX(-0.000011179461895408, -0.000019851294687597),
CMPLX(-0.000000946901382392, 0.000001070014057386),
CMPLX( 0.000000111032677121, 0.000000038305261714),
CMPLX(-0.000000001294398927, -0.000000011628723277),
CMPLX(-0.000000001114905944, 0.000000000759733092),
CMPLX( 0.000000000157637232, 0.000000000075476075),
CMPLX(-0.000000000002830457, -0.000000000024752781),
CMPLX(-0.000000000002932169, 0.000000000002493893),
CMPLX( 0.000000000000617809, 0.000000000000156198),
CMPLX(-0.000000000000043162, -0.000000000000103385),
CMPLX(-0.000000000000010133, 0.000000000000018129),
CMPLX( 0.000000000000003973, -0.000000000000000709),
CMPLX(-0.000000000000000632, -0.000000000000000520),
CMPLX( 0.000000000000000006, 0.000000000000000172),
CMPLX( 0.000000000000000027, -0.000000000000000026),
CMPLX(-0.000000000000000008, -0.000000000000000000),
CMPLX( 0.000000000000000001, 0.000000000000000001)
};
BEGIN
IF (X <= 0.0) THEN
ErrOut('NON-POSITIVE ARGUMENT X (BesselY1)');
RETURN 0.0;
END;
IF (ABS(X) <= 2.0*MachEps) THEN
RETURN 0.0;
END;
v := ABS(X);
IF (v < 8.0) THEN
y := v / 8.0;
h := 2.0*y*y - 1.0;
alfa := h+h;
b1:=0.0; b2:=0.0;
FOR i:=17 TO 1 BY -1 DO
b0 := D[i] + alfa*b1 - b2;
b2 := b1;
b1 := b0;
END;
p := y*(b0-h*b2);
b1:=0.0; b2:=0.0;
FOR i:=17 TO 1 BY -1 DO
b0 := E[i] + alfa*b1 - b2;
b2 := b1;
b1 := b0;
END;
p := PI1*((CE + ln(0.5*X))*p - 1.0/X) + y*(b0-b2);
ELSE
r := 1.0 / v;
h := 10.0*r - 1.0;
alfa := h+h;
cb1 := zero;
cb2 := zero;
FOR i:=20 TO 1 BY -1 DO
cb0 := F[i] + scalarMult(alfa,cb1) - cb2;
cb2 := cb1;
cb1 := cb0;
END;
cb0 := scalarMult(sqrt(PI1*r),CEXP(imag*CMPLX(v-PI3,0.0))) *
(cb0 - scalarMult(h,cb2)); (*red?*)
p := RE(-imag*cb0); (*red?*)
END;
IF (X < 0.0) THEN p := -p; END;
RETURN p;
END DBesselY1;
PROCEDURE BesselY01( x : LONGREAL;
VAR y0,y1 : LONGREAL);
TYPE Vek15 = ARRAY [1..15] OF LONGREAL;
CONST ar1 = Vek15{
+0.164349E-14,
-0.8747341E-13, +0.402633082E-11,
-0.15837552542E-09, +0.524879478733E-08,
-0.14407233274019E-06, +0.32065325376548E-05,
-0.563207914105699E-04, +0.753113593257774E-03,
-0.72879624795521E-02, +0.471966895957634E-01,
-0.177302012781143E-00, +0.261567346255047E-00,
+0.179034314077182E-00, -0.274474305529745
};
CONST ar2 = Vek15{
+0.42773E-15, -0.2440949E-13,
+0.121143321E-11, -0.5172121473E-10,
+0.187547032473E-08, -0.5688440039919E-07,
+0.141662436449235E-05, -0.283046401495148E-04,
+0.440478629867099E-03, -0.51316411610611E-02,
+0.423191803533369E-01, -0.226624991556755E-00,
+0.675615780772188E-00, -0.767296362886646E-00,
-0.128697384381350E-00
};
CONST c0 = 0.636619772367581;
VAR z,z2,c,lnx,b0,b1,b2 : LONGREAL;
cosx,sinx,p0,q0,p1,q1 : LONGREAL;
i : CARDINAL;
BEGIN
IF debug THEN TIO.WrStr("Aufruf von bessy01"); TIO.WrLn; END;
IF (x < 8.0) THEN
lnx:= c0*ln(x);
c:= c0/x; x:= x/8.0; z := 2.0*x*x - 1.0; z2:= z+z;
(* computation of y0 *)
b1:=0.0; b2:= 0.0;
FOR i:=1 TO 15 DO
b0:= z2*b1 - b2 + ar1[i];
b2:= b1; b1:= b0;
END;
y0:= lnx*BesselJ0(8.0*x) + z*b1 - b2 - 0.33146113203285E-01;
(* computation of y1; *)
b1:=0.0; b2:=0.0;
FOR i:=1 TO 15 DO
b0:= z2*b1 - b2 + ar2[i];
b2:= b1; b1:= b0;
END;
y1 := lnx*BesselJ1(x*8.0) - c + x*(z*b1 - b2 + 0.2030410588593425E-01);
ELSE
c := 0.797884560802865 / sqrt(x);
besspq0(x,p0,q0);
besspq1(x,p1,q1);
x:= x - 0.706858347057703E01; cosx := cos(x); sinx := sin(x);
y0 := c*(p0*sinx + q0*cosx);
y1 := c*(q1*sinx - p1*cosx);
END;
END BesselY01;
MODULE CHank;
IMPORT pi,sqrt;
IMPORT one,CABS,CEXP,scalarMult;
IMPORT ErrOut;
EXPORT CHank01Asym;
VAR Pvek,Qvek : ARRAY [1..32] OF LONGREAL;
PROCEDURE CHank01Asym( Z : LONGREAL;
RelErr : LONGREAL;
VAR H01 : LONGCOMPLEX);
(*----------------------------------------------------------------*)
(* This routine uses Hankels asymptotic expansion [2,(9.2.5)] *)
(* to calculate the Hankel functions of orders 0 and argument Z *)
(* within the desired relative error RelErr. *)
(*----------------------------------------------------------------*)
CONST imag = CMPLX(0.0,1.0);
VAR pk0,qk0,zsq,sqtmp : LONGCOMPLEX;
term01,ex1,xi : LONGCOMPLEX;
k : INTEGER;
BEGIN
zsq := CMPLX(Z,0.0)*CMPLX(Z,0.0);
pk0 := one;
qk0 := CMPLX(-1.0/(8.0*Z),0.0);
xi := CMPLX(Z - pi/4.0 ,0.0);
ex1 := CEXP(imag*xi);
H01 := one + imag*qk0;
k:=1;
REPEAT
pk0 := scalarMult(Pvek[k],pk0) / zsq;
qk0 := scalarMult(Qvek[k],qk0) / zsq;
term01 := pk0 + imag*qk0;
H01 := H01 + term01;
INC(k);
UNTIL (CABS(term01/H01) < RelErr) OR (k > 32);
IF (k > 32) THEN
ErrOut('Not enough terms in HankelH01');
END;
sqtmp := CMPLX(sqrt(2.0/(pi*Z)),0.0);
IF (RE(sqtmp) < 0.0) THEN sqtmp := -sqtmp; END;
H01 := sqtmp*H01*ex1;
END CHank01Asym;
VAR k : INTEGER ;
k2,t1,t2,xi,tmp : LONGREAL;
BEGIN
(* Set up the constants that are used for hankel's *)
(* asymptotic expansion (see [3], (9.2.5)) *)
t2 := -1.0;
xi := 1.0;
FOR k:=1 TO 32 DO
k2 := 2.0*LFLOAT(k);
xi := xi + 2.0;
t1 := t2;
t2 := -xi*xi;
tmp := k2*(k2-1.0)*64.0;
(* pqarray(1,k) = (0.2*k)*(-1/4)**k, where hankel's symbol is *)
(* defined as (n,k) = gamma(.5+n+k)/(gamma(k+1)*gamma(.5+n-k) *)
Pvek[k] := (-1.0*t1*t2)/tmp;
xi := xi + 2.0;
t1 := t2;
t2 := -xi*xi;
tmp := k2*(k2+1.0)*64.0;
(* pqarrya(2,k) = 0.5*(0.2*k+1)*(-1/4)**k *)
Qvek[k] := (-1.0*t1*t2)/tmp;
END;
END CHank; (* END MODULE *)
PROCEDURE Hankel01(x : LONGREAL) : LONGCOMPLEX;
CONST C = LMathLib.C; (* Eulerkonstante *)
VAR xh,xh2,Zi,SumR,SumI,Phi,yf2,Fakt,y : LONGREAL;
res : LONGCOMPLEX;
i : CARDINAL;
BEGIN
IF (x <= 0.0) THEN
ErrOut('x = 0.0 (SpezFunkt3.Hankel01)');
RETURN CMPLX(0.0,0.0);
END;
IF (x <= 8.0) THEN
xh := 0.5*x; xh2 := xh*xh;
SumR := 1.0; SumI := 0.0;
Fakt := 1.0; Phi := 1.0;
i:=1; y := xh2;
LOOP
Zi := SumI;
yf2 := y / (Fakt*Fakt);
IF ODD(i) THEN
SumR:=SumR - yf2;
SumI:=SumI - 2.0*yf2*Phi;
ELSE
SumR:=SumR + yf2;
SumI:=SumI + 2.0*yf2*Phi;
END;
IF (Zi = SumI) THEN EXIT; END; (* !!!! *)
INC(i);
Fakt := Fakt*VAL(LONGREAL,i);
Phi := Phi + 1.0/VAL(LONGREAL,i);
y:=y*xh2;
END; (* LOOP *)
res := CMPLX(SumR , ((2.0*SumR*(C + ln(xh)) - SumI) / pi));
ELSIF (x <= 17.5) THEN
res := CMPLX(BesselJ0(x),BesselY0(x));
ELSE
CHank01Asym(x,MachEps,res);
END;
RETURN res;
END Hankel01;
END SpezFunkt3.