DIRECTORY RealFns, RealHypBessel; RealHypBesselImpl: CEDAR PROGRAM IMPORTS RealFns EXPORTS RealHypBessel = BEGIN SinH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN ex: REAL = RealFns.Exp[x]; emx: REAL = 1/ex; f _ (ex-emx)/2.0 END; CosH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN ex: REAL = RealFns.Exp[x]; emx: REAL = 1/ex; f _ (ex+emx)/2.0 END; TanH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF ABS[x] > 20. THEN f _ (IF x<0. THEN -1. ELSE 1.) ELSE BEGIN ex: REAL = RealFns.Exp[x]; emx: REAL = 1/ex; f _ (ex-emx)/(ex+emx); END; END; CotH: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF ABS[x] > 20. THEN f _ (IF x<0. THEN -1. ELSE 1.) ELSE BEGIN ex: REAL = RealFns.Exp[x]; emx: REAL = 1/ex; f _ (ex+emx)/(ex-emx); END; END; pi: REAL _ 3.1415926536; pi4: REAL _ 0.7853981633; J0: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN IF x < 0.0 THEN x _ -x; IF x < 8.0 THEN BEGIN x2: REAL = x*x; P: REAL = (((( (-0.1849052456E3)*x2+ (0.7739233017E5))*x2+ (-0.1121442418E8))*x2+ (0.6516196407E9))*x2+ (-0.1336259035E11))*x2+ (0.5756849057E11); Q: REAL = (((( x2+ (0.2678532712E3))*x2+ (0.5927264853E5))*x2+ (0.9494680718E7))*x2+ (0.1029532985E10))*x2+ (0.5756849041E11); f _ P/Q; END ELSE BEGIN X0: REAL = x-pi4; p0, q0: REAL; [p0, q0] _ PQ0[x]; f _ RealFns.SqRt[2.0/(pi*x)]*(p0*RealFns.Cos[radians: X0]-q0*RealFns.Sin[radians: X0]); END; END; PQ0: PROC [x: REAL] RETURNS [P0, Q0: REAL] = BEGIN z: REAL = 8.0/x; z2: REAL = z*z; P0 _ (( -- PZERO 6502 - 8.78 digits (-0.165711693E-5)*z2+ (0.270871648E-4))*z2+ (-0.1098577711E-2))*z2+ (0.9999999983); Q0 _ ((( -- QZERO 6901 - 9.13 digits (0.576476505E-6)*z2+ (-0.6796319617E-5))*z2+ (0.1430262718E-3))*z2+ (-0.1562499927E-1))*z; END; pi34: REAL _ 3*pi4; J1: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = BEGIN sign: REAL _ 1.0; IF x < 0.0 THEN {sign _ -1.0; x _ -x}; IF x < 8.0 THEN BEGIN x2: REAL = x*x; P: REAL = ((( (0.6543001487E1)*x2+ (-0.1988185996E4))*x2+ (0.1978615055E6))*x2+ (-0.7064812917E7))*x2+ (0.6706259146E8); Q: REAL = ((( x2+ (0.1766217266E3))*x2+ (0.2668508067E5))*x2+ (0.2635953855E7))*x2+ (0.1341252501E9); f _ sign*x*P/Q; END ELSE BEGIN X1: REAL = x-pi34; p1, q1: REAL; [p1, q1] _ PQ1[x]; f _ sign*RealFns.SqRt[2.0/(pi*x)]*(p1*RealFns.Cos[radians: X1]-q1*RealFns.Sin[radians: X1]); END; END; PQ1: PROC [x: REAL] RETURNS [P1, Q1: REAL] = BEGIN z: REAL = 8.0/x; z2: REAL = z*z; P1 _ (( -- PONE 6702 - 8.72 digits (0.19796752E-5)*z2+ (-0.348677998E-4))*z2+ (0.1830991520E-2))*z2+ (0.1000000001E1); Q1 _ ((( -- QONE 7102 - 9.08 digits (-0.67222066E-6)*z2+ (0.831923005E-5))*z2+ (-0.2002434943E-3))*z2+ (0.4687499917E-1))*z; END; Jn: PUBLIC PROC [n: INT, x: REAL] RETURNS [f: REAL] = BEGIN sign: REAL _ 1.0; IF n<0 THEN {n _ -n; sign _ IF n MOD 2 # 0 THEN -1.0 ELSE 1.0}; SELECT n FROM 0 => f _ J0[x]; 1 => f _ sign*J1[x]; ENDCASE => BEGIN v0: REAL _ J0[x]; v1: REAL _ J1[x]; FOR m: INT IN [2..n] DO v2: REAL _ 2*(m-1)*v1/x-v0; v0 _ v1; v1 _ v2; ENDLOOP; f _ sign*v1; END; END; Y0: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = {ERROR}; Y1: PUBLIC PROC [x: REAL] RETURNS [f: REAL] = {ERROR}; Yn: PUBLIC PROC [n: INT, x: REAL] RETURNS [f: REAL] = {ERROR}; END. ╚RealHypBesselImpl.mesa Last Update by E. McCreight, November 13, 1984 12:48:32 pm PST Bessel function approximations are taken from the book "Computer Approximations", by Hart, Cheney, Lawson, Maehly, Mesztenyi, Rice, Thacher, and Witzgall, published by Wiley. Function names and indexes refer to that book. Bessel function of first kind of order 0 JZERO 5838 - 8.54 digits Bessel function of first kind of order 1 JONE 6036 - 7.08 digits ..actually, JONE 6038 - 8.19 digits would be better, but that page is missing from my copy. Bessel function of first kind of order n Bessel function of second kind of order 0 Bessel function of second kind of order 1 Bessel function of second kind of order n Э╩0ШJЪЬЩJЪЬ>Щ>JШЪ╧k Ш JЪЬШJЪЬШJШЧЪ ЬЭЬЭЬЭЬ ЭЬШHJЪЭШJШЪ╧nЬЭЬЭЬЭЬЭЬЭЬШ/JЪЭШJЪЬЭЬШJЪЬЭЬШJЪЬШJЪЭЬШJШЧЪЮЬЭЬЭЬЭЬЭЬЭЬШ/JЪЭШJЪЬЭЬШJЪЬЭЬШJЪЬШJЪЭЬШJШЧЪЮЬЭЬЭЬЭЬЭЬЭЬШ/JЪЭШJЪЭЬЭЬ ЭЬЭЬЭЬЭЬШ3ЪЭШJЪЭШJЪЬЭЬШJЪЬЭЬШJЪЬШJЪЭЬШЧJЪЭЬШJШЧЪЮЬЭЬЭЬЭЬЭЬЭЬШ/JЪЭШJЪЭЬЭЬ ЭЬЭЬЭЬЭЬШ3ЪЭШJЪЭШJЪЬЭЬШJЪЬЭЬШJЪЬШJЪЭЬШЧJЪЭЬШJШJШЧJШJЪЬ▐Щ▐JШJЪЬЭЬШJЪЬЭЬШJШЪЮЬЭЬЭЬЭЬЭЬЭЬШ-JЪЬ(Щ(JЪЭШJЪЭЬ ЭЬШЪЭЬ ЭШJЪЭШJЪЬЩJЪЬЭЬШЪЭЬЭЬШJЪЬШJЪЬШJЪЬШJЪЬШJЪЬШJЪЬШЧЪЭЬЭЬШJЪЬШJЪЬШJЪЬШJЪЬШJЪЬШJЪЬШЧJЪЬЭЬЭЬШJЪЭШЧЪЭШJЪЭШJЪЬЭЬ ШJЪЬЭЬШ JЪЬШJЪЬWШWJЪЭЬШЧJЪЭЬШJЩЧЪ ЮЬЭЬЭЬЭЬ ЭЬШ,JЪЭШJЪЬЭЬ ШJЪЬЭЬШЪЬ╧cШ#JЪЬШJЪЬШJЪЬШJЪЬШЧЪЬ ЯШ$JЪЬШJЪЬШJЪЬШJЪЬШЧJЪЭЬШJШЧJЪЬЭЬ ШJШЪЮЬЭЬЭЬЭЬЭЬЭЬШ-JЪЬ(Щ(JЪЭШJЪЬЭЬШJЪЭЬ ЭЬШ&ЪЭЬ ЭШJЪЭШЪЬЩJЩ[ЧJЪЬЭЬШЪЭЬЭЬШ JЪЬШJЪЬШJЪЬШJЪЬШJЪЬШЧЪЭЬЭЬШ JЪЬШJЪЬШJЪЬШJЪЬШJЪЬШЧJЪЬЭЬЭЬШJЪЭШЧЪЭШJЪЭШJЪЬЭЬ ШJЪЬЭЬШ JЪЬШJЪЬ\Ш\JЪЭЬШЧJЪЭЬШJЩЧЪ ЮЬЭЬЭЬЭЬ ЭЬШ,JЪЭШJЪЬЭЬ ШJЪЬЭЬШЪЬЯШ"JЪЬШJЪЬШJЪЬШJЪЬШЧЪЬ ЯШ#JЪЬШJЪЬШJЪЬШJЪЬШЧJЪЭЬШJЩЧЪЮЬЭЬЭЬЭЬЭЬЭЬЭЬШ5JЪЬ(Щ(JЪЭШJЪЬЭЬШJЪЭЬЭЬЭЬЭЬЭЬЭЬШ?ЪЭЬЭШ JШJШЪЭЬШ JЪЭШJЪЬЭЬ ШJЪЬЭЬ ШЪЭЬЭЬЭЬЭШJЪЬЭЬШJЪЬШJЪЬШJЪЭЬШЧJШJЪЭЬШЧЧJЪЭЬШJШJЩЧЪЮЬЭЬЭЬЭЬЭЬЭЬЭЬШ6JЪЬ)Щ)JЩЧЪЮЬЭЬЭЬЭЬЭЬЭЬЭЬШ6JЪЬ)Щ)JЩЧЪЮЬЭЬЭЬЭЬЭЬЭЬЭЬЭЬШ>JЪЬ)Щ)JЩЧJЪЭЬШЧJШЧЕЧJB