home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / TSTBES.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-18  |  1.6 KB  |  83 lines

  1. program tstbes;  { -> 344 }
  2. { test the bessel function }
  3. { the Gamma function is included }
  4.  
  5. var done :boolean;
  6.  x,ordr : real;
  7.  
  8.  
  9. function gamma(x: real): real;
  10. const pi = 3.1415926;
  11.  
  12. var i,j : integer;
  13.  y,gam : real;
  14.  
  15. begin  { gamma function }
  16.   if x>=0.0 then
  17.     begin
  18.       y:=x+2.0;
  19.       gam:=sqrt(2*pi/y)*exp(y*ln(y)+(1-1/(30*y*y))/(12*y)-y);
  20.       gamma:=gam/(x*(x+1))
  21.     end
  22.   else  { x<0 }
  23.     begin
  24.       j:=0;
  25.       y:=x;
  26.       repeat
  27.  j:=j+1;
  28.  y:=y+1.0
  29.       until y>0.0;
  30.       gam:=gamma(y);  { recursive call }
  31.       for i:=0 to j-1 do
  32.  gam:=gam/(x+1);
  33.       gamma:=gam
  34.     end { x<0 }
  35. end;  { gamma function }
  36.  
  37. function bessj(x,n: real): real;
  38. { cylindrical Bessel function of the first kind }
  39. { the gamma function is required }
  40.  
  41. const tol = 1.0E-4;
  42.  pi = 3.1415926;
  43.  
  44. var i  : integer;
  45.  term,new_term,
  46.  sum,x2  : real;
  47.  
  48. begin { bessj }
  49.   x2:=x*x;
  50.   if (x=0.0)and(N=1.0) then bessj:=0.0
  51.   else if x>15 then { asymptotic expansion }
  52.     bessj:=sqrt(2/(pi*x))*cos(x-pi/4-n*pi/2)
  53.   else
  54.     begin
  55.       if n=0.0 then sum:=1.0
  56.       else sum:=exp(n*ln(x/2))/gamma(n+1.0);
  57.       new_term:=sum;
  58.       i:=0;
  59.   repeat
  60.     i:=i+1;
  61.     term:=new_term;
  62.     new_term:=-term*x2*0.25/(i*(n+1));
  63.     sum:=sum+new_term
  64.   until abs(new_term)<=abs(sum*tol);
  65.   bessj:=sum
  66.  end { if}
  67. end; { bessj }
  68.  
  69. begin  { main }
  70.   done:=false;
  71.   repeat
  72.     write('Order: ');
  73.     readln(ordr);
  74.     if ordr<-25.0 then done:=true
  75.     else
  76.       begin
  77.  write('X: ');
  78.  readln(x);
  79.  writeln('J Bessel is ',bessj(x,ordr))
  80.      end
  81.   until done
  82. end.
  83.