code 35162;
    procedure 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 
    begin real X2, R, S; integer L, M, NU, SIGNX;
        SIGNX:= SIGN(X); X:= ABS(X);
        R:= S:= 0; X2:= 2/X; L:= 0; NU:= START(X,N,0);
        for M:= NU    step -1 until 1 do 
        begin R:= 1/(X2*M-R);
            L:= 2-L; S:= R*(L+S);
            if M<=N then J[M]:= R
        end;
        J[0]:= R:= 1/(1+S);
        for M:= 1 step 1 until N do 
        J[M]:= R:= R*J[M];
        if SIGNX < 0 then 
        for M:= 1 step 2 until N do 
        J[M]:= -J[M];
    end BESSELJ;
        eop