--- a/SpezFunkt3.mod
+++ b/SpezFunkt3.mod
@@ -22,6 +22,7 @@
(* 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 *)
@@ -39,7 +40,7 @@
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
- (* $Id: SpezFunkt3.mod,v 1.5 2018/05/16 07:08:36 mriedl Exp mriedl $ *)
+ (* $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;
@@ -54,6 +55,14 @@
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);
@@ -64,10 +73,10 @@
(* Fortran bearbeitet : MRi, 07.05.18 *)
(*----------------------------------------------------------------*)
- CONST R32 = 1.0/32.0;
+ CONST R32 = 1.0 / 32.0;
PI1 = 2.0 / pi;
PI4 = pi / 4.0;
- I = CMPLX(0.0,1.0);
+ imag = CMPLX(0.0,1.0);
TYPE RVek17 = ARRAY [1..17] OF LONGREAL;
@@ -127,15 +136,15 @@
r := 1.0 / v;
h := 10.0*r - 1.0;
alfa := h + h;
- cb1 := CMPLX(0.0,0.0);
- cb2 := CMPLX(0.0,0.0);
+ 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 := LongComplexMath.exp(I*CMPLX((v - PI4),0.0));
- cb0 := scalarMult(sqrt(PI1*r),ctmp)*(cb0 - scalarMult(h,cb2));
+ 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;
@@ -183,14 +192,13 @@
x2,b0,b1,b2,y : LONGREAL;
i : CARDINAL;
BEGIN
- IF debug THEN TIO.WrStr("Aufruf von besspq0"); END;
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);
+ 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;
@@ -208,7 +216,6 @@
END;
q0:=(x*b1 - b2 -0.015555854605337)*y
END;
- IF debug THEN TIO.WrLn; END;
END besspq0;
PROCEDURE BesselJ0(x : LONGREAL) : LONGREAL;
@@ -244,7 +251,7 @@
END;
RETURN z*b1 - b2 + 0.15772797147489;
ELSE
- x := ABS(x); (* Codefragment korrigiert, nicht original NuMaL *)
+ 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);
@@ -291,7 +298,7 @@
x2,b0,b1,b2,y : LONGREAL;
i : CARDINAL;
BEGIN
- IF debug THEN TIO.WrStr("Aufruf bon besspq1"); TIO.WrLn; END;
+ 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);
@@ -300,7 +307,7 @@
p1:= b*(j1x*sinx - y1*cosx);
q1:= b*(j1x*cosx + y1*sinx)
ELSE
- y:= 8 / x; x:= 2 * y * y - 1; x2 := x + x;
+ 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
@@ -341,7 +348,7 @@
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;
+ 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
@@ -575,10 +582,10 @@
PROCEDURE DBesselY1(X : LONGREAL) : LONGREAL;
- CONST I = CMPLX(0.0,1.0);
- CE = 0.577215664901532861;
- PI1 = 2.0 / pi;
- PI3 = 3.0*pi / 4.0;
+ 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;
@@ -634,54 +641,53 @@
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(I*CMPLX(v-PI3,0.0))) *
- (cb0 - scalarMult(h,cb2));
- p := RE(-I*cb0);
-
- END;
-
- IF (X < 0.0) THEN p := -p; END;
-
- RETURN p;
+ 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;
@@ -721,21 +727,21 @@
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;
+ 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*x) + z*b1 - b2 - 0.33146113203285E-01;
+ 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) - c + x*(z*b1 - b2 + 0.2030410588593425E-01);
+ y1 := lnx*BesselJ1(x*8.0) - c + x*(z*b1 - b2 + 0.2030410588593425E-01);
ELSE
c := 0.797884560802865 / sqrt(x);
besspq0(x,p0,q0);
@@ -749,7 +755,7 @@
MODULE CHank;
IMPORT pi,sqrt;
-IMPORT CABS,CEXP,scalarMult;
+IMPORT one,CABS,CEXP,scalarMult;
IMPORT ErrOut;
EXPORT CHank01Asym;
@@ -765,8 +771,7 @@
(* within the desired relative error RelErr. *)
(*----------------------------------------------------------------*)
- CONST IM = CMPLX(0.0,1.0);
- one = CMPLX(1.0,0.0);
+ CONST imag = CMPLX(0.0,1.0);
VAR pk0,qk0,zsq,sqtmp : LONGCOMPLEX;
term01,ex1,xi : LONGCOMPLEX;
@@ -775,15 +780,15 @@
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(IM*xi);
- H01 := one + IM*qk0;
+ 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 + IM*qk0;
+ term01 := pk0 + imag*qk0;
H01 := H01 + term01;
INC(k);
UNTIL (CABS(term01/H01) < RelErr) OR (k > 32);
@@ -808,17 +813,17 @@
t1 := t2;
t2 := -xi*xi;
tmp := k2*(k2-1.0)*64.0;
- (* pqarry(1,k) = (0,2*k)*(-1/4)**k, where hankel's symbol is *)
+ (* 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;
- (* pqarry(2,k) = 0.5*(0,2*k+1)*(-1/4)**k *)
+ (* pqarrya(2,k) = 0.5*(0.2*k+1)*(-1/4)**k *)
Qvek[k] := (-1.0*t1*t2)/tmp;
END;
-END CHank;
+END CHank; (* END MODULE *)
PROCEDURE Hankel01(x : LONGREAL) : LONGCOMPLEX;
@@ -833,10 +838,10 @@
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;
+ 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);
@@ -847,7 +852,7 @@
SumR:=SumR + yf2;
SumI:=SumI + 2.0*yf2*Phi;
END;
- IF (Zi = SumI) THEN EXIT END; (* !!!! *)
+ IF (Zi = SumI) THEN EXIT; END; (* !!!! *)
INC(i);
Fakt := Fakt*VAL(LONGREAL,i);
Phi := Phi + 1.0/VAL(LONGREAL,i);