code 35180;
  procedure BESS JAPLUSN(A, X, N, JA); value A, X, N;
  integer N; real X, A; array JA;
  if X = 0 then 
  begin JA[0]:= if A = 0 then 1 else 0;
    for N:= N step -1 until 1 do JA[N]:= 0
  end else 
  if A = 0 then 
  begin 
    BESS J(X, N, JA)
  end else 
  if A = .5 then 
  begin real S;
    S:= SQRT(X) * .797 884 560 802 865; comment S = SQRT(2X / PI);
    SPHER BESS J(X, N, JA);
    for N:= N step - 1 until 0 do JA[N]:= JA[N] * S
  end 

  else 
  begin real A2, X2, R, S, L, LABDA; integer K, M, NU;
    L:= 1; NU:= START(X,N,0);
    for M:= 1 step 1 until NU do 
    L:= L * (M+A) / (M+1); R:= S:= 0; X2:= 2 / X; K:= -1; A2:= A + A;
    for M:= NU+NU step - 1 until 1 do 
    begin R:= 1 / (X2 * (A + M)    - R);
      if K = 1 then LABDA:= 0 else 
      begin L:= L * (M + 2) / (M + A2); LABDA:= L * (M + A) end;
      S:= R * (LABDA + S); K:= -K;
      if M<= N then JA[M]:= R
    end;
    JA[0]:= R:= 1 / GAMMA(1 + A) / (1 + S) / X2 ** A;
    for M:= 1 step 1 until N do JA[M]:= R:= R * JA[M];
  end BESS JAPLUSN;
        eop