home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / BESY.PAS next >
Encoding:
Pascal/Delphi Source File  |  1985-07-18  |  1.7 KB  |  90 lines
  1. program besy;  { -> 348 }
  2. { evaluation of Bessel function of the second kind }
  3.  
  4. var x,ordr : real;
  5.  done : boolean;
  6.  
  7.  
  8. function bessy(x,n: real): real;
  9. { cylindical bessel function of the second kind }
  10. const small = 1.0E-8;
  11.  euler = 0.57721566;
  12.  pi = 3.1415926;
  13.  pi2 = 0.63661977; { 2/pi }
  14. var j : integer;
  15.  
  16.  x2,sum,sum2,t,t2,
  17.  ts,term,xx,y0,y1,
  18.  ya,yb,yc,ans,a,b,
  19.  sina,cosa  : real;
  20.  
  21. begin   { function bessy }
  22.   if x<12 then
  23.     begin
  24.       xx:=0.5*x;
  25.       x2:=xx*xx;
  26.       t:=ln(xx)+euler;
  27.       sum:=0.0;
  28.       term:=t;
  29.       y0:=t;
  30.       j:=0;
  31.       repeat
  32.  j:=j+1;
  33.  if j<>1 then sum:=sum+1/(j-1);
  34.  ts:=t-sum;
  35.  term:=-x2*term/(j*j)*(1-1/(j*ts));
  36.  y0:=y0+term
  37.       until abs(term)<small;
  38.       term:=xx*(t-0.5);
  39.       sum:=0.0;
  40.       y1:=term;
  41.       j:=1;
  42.       repeat
  43.  j:=j+1;
  44.  sum:=sum+1/(j-1);
  45.  ts:=t-sum;
  46.  term:=(-x2*term)/(j*(j-1))*((ts-0.5/j)/(ts+0.5/(j-1)));
  47.  y1:=y1+term
  48.       until abs(term)<small;
  49.       y0:=pi2*y0;
  50.       y1:=pi2*(y1-1/x);
  51.       if n=0.0 then ans:=y0
  52.       else if n=1.0 then ans:=y1
  53.       else
  54.  begin  { find y by recursion }
  55.    ts:=2.0/x;
  56.    ya:=y0;
  57.    yb:=y1;
  58.    for j:=2 to trunc(n+0.01) do
  59.      begin
  60.        yc:=ts*(j-1)*yb-ya;
  61.        ya:=yb;
  62.        yb:=yc
  63.      end;
  64.    ans:=yc
  65.  end;
  66.       bessy:=ans;
  67.     end  { x<12 }
  68.   else  { x>11, asymtotic expansion }
  69.     bessy:=sqrt(2/(pi*x))*sin(x-pi/4-n*pi/2)
  70. end; { function bessy }
  71.  
  72. begin
  73.   clrscr;
  74.   done:=false;
  75.   writeln;
  76.   repeat
  77.     write('Order? ');
  78.     readln(ordr);
  79.     if ordr<0.0 then done:=true
  80.     else
  81.       begin
  82.  repeat
  83.    write('Arg? ');
  84.    readln(x)
  85.  until x>=0.0;
  86.       writeln('Y Bessel is ',bessy(x,ordr))
  87.     end  { if }
  88.   until done
  89. end.
  90.