code 35150;
  comment SPHERICAL BESSEL FUNCTIONS J[.5](X),  , J[N+.5](X);
  procedure SPHER BESS J(X, N, J); value X, N;
  real X; integer N; array J;
  if X = 0 then 
  begin J[0]:= 1;
    for N:= N step -1 until 1 do J[N]:=0
  end else if N = 0 then 
  begin real X2;
    if ABS(X) < .015 then 
    begin X2:= X * X / 6; J[0]:= 1 + X2 * (X2 * .3 - 1) end else 
    J[0]:= SIN(X)/X
  end else 
  begin integer M; real R, S;
    R:= 0; M:= START(X,N,0);
    for M:= M step - 1 until 1 do 
    begin R:= 1 / ((M + M + 1) / X - R); if M <= N then J[M]:= R
    end; if X < .015  then 
    begin S:= X * X / 6;
      J[0]:= R:= S * (S * .3 - 1) + 1 end else 
    J[0]:= R:= SIN(X) / X;
    for M:= 1 step 1 until N do J[M]:= R:= J[M] * R;
  end SPHER BESS J;
        eop