program besy uses sampler var ordr res real integer const PName String

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
program besy;
uses sampler;
var
x, ordr, res: real;
i: integer;
const
PName: String[15] = 'orig_f.pas';
function bessy(x, n: real): real;
const
small = 0.00000001;
euler = 0.57721566;
pi = 3.1415926;
pi2 = 0.63661977;
var
j: integer;
x2, sum, sum2, t, t2,
ts, term, xx, y0, y1,
ya, yb, yc, ans, a, b,
sina, cosa: real;
begin
if x < 12 then
begin
xx := 0.5 * x;
x2 := xx * xx;
t := ln(xx) + euler;
sum := 0.0;
j := 0;
term := t;
y0 := t;
repeat
j := j + 1;
if j <> 1 then sum := sum + 1 / (j - 1);
ts := t - sum;
term := -x2 * term / (j * j) * (1 - 1 / (j * ts));
y0 := y0 + term;
until abs(term) < small;
term := xx * (t - 0.5);
sum := 0.0;
y1 := term;
j := 1;
repeat
j := j + 1;
sum := sum + 1 / (j - 1);
ts := t - sum;
term := (-x2 * term) / (j * (j - 1)) * ((ts - 0.5 / j) / (ts + 0.5 / (j - 1)));
y1 := y1 + term;
until abs(term) < small;
y0 := pi2 * y0;
y1 := pi2 * (y1 - 1 / x);
if n = 0.0 then ans := y0
else if n = 1.0 then ans := y1
else
begin
ts := 2.0 / x;
ya := y0;
yb := y1;
for j := 2 to trunc(n + 0.01) do
begin
yc := ts * (j - 1) * yb - ya;
ya := yb;
yb := yc
end;
ans := yc
end;
bessy := ans;
end
else
bessy := sqrt(2 / (pi * x)) * sin(x - pi / 4 - n * pi / 2);
end;
begin
Sample(PName, 88);
ordr := 15;
x := 5;
res := bessy(x, ordr);
Sample(PName, 92);
end.