code 35184;
procedure BESS ZEROS(A,N,Z,D); value A,N,D; real A;array Z;
                                 integer N,D;
comment COMPUTES Z[1],...Z[N],THE FIRST N ZEROS OF A BESSEL FUNCTION.
   THE CHOICE OF D DETERMINES THE TYPE OF THE BESSEL FUNCTION :
   IF D=1 THEN JA       ELSE
   IF D=2 THEN YA       ELSE
   IF D=3 THEN JA-PRIME ELSE
   IF D=4 THEN YA-PRIME.
   A IS THE ORDER OF THE BESSEL FUNCTION, IT MUST BE NON-NEGATIVE.;
begin real AA,A2,B,BB,C,CHI,CO,MU,MU2,MU3,MU4,P,PI,PA,PA1,P0,P1,PP1,
             Q,QA,QA1,Q1,QQ1,RO,SI,T,TT,U,V,W,X,XX,X4,Y; integer J,S;

   real procedure FI(Y); value Y; real Y;
   comment COMPUTES FI FROM THE EQUATION
          TAN(FI)-FI=Y , WHERE Y>=0.
   THE RELATIVE ACCURACY IS AT LEAST 5 DIGITS;
   if Y=0  then FI:=0         else 
   if Y>"5 then FI:=1.570796  else 
   begin real R,P,PP;
       if Y<1 then 
       begin P:=(3*Y)**(1/3); PP:=P*P;
           P:=P*(1+PP*(-210+PP*(27-2*PP))/1575)
       end else 
       begin P:=1/(Y+1.570796); PP:=P*P;
           P:= 1.570796-P*(1+PP*(2310+PP*(3003+PP*(4818+PP*
                   (8591+PP*16328))))/3465)
       end;
       PP:=(Y+P)*(Y+P); R:=(P-ARCTAN(P+Y))/PP;
       FI:=P-(1+PP)*R*(1+R/(P+Y))
   end FI;

   real procedure R;
   begin BESS PQA01(A,X,PA,QA,PA1,QA1);
       CHI:=X-PI*(A/2+0.25);
       SI :=SIN(CHI); CO:=COS(CHI);
       R:= if D=1 then (PA*CO-QA*SI)/(PA1*SI+QA1*CO) else 
           if D=2 then (PA*SI+QA*CO)/(QA1*SI-PA1*CO) else 
           if D=3 then A/X-(PA1*SI+QA1*CO)/(PA*CO-QA*SI) else 
                           A/X-(QA1*SI-PA1*CO)/(PA*SI+QA*CO)
   end R;
   PI:=4*ARCTAN(1); AA:=A*A; MU:=4*AA; MU2:=MU*MU;
   MU3:=MU*MU2; MU4:=MU2*MU2;
   if D<3 then 
   begin P:=7*MU-31; P0:=MU-1;
       P1:=4*(253*MU2-3722*MU+17869)/15/P*P0;
       Q1:=8*( 83*MU2- 982*MU+ 3779)/ 5/P
   end else 
   begin P:=7*MU2+82*MU-9; P0:=MU+3;
       P1:=(4048*MU4+131264*MU3-221984*MU2-417600*MU+1012176)/60/P;
       Q1:=1.6*(83*MU3+2075*MU2-3039*MU+3537)/P
   end;
   T:=if D=1or D=4 then 0.25 else 0.75; TT:=4*T;
   if D<3 then 
   begin PP1:= 5/48; QQ1:= -5/36  end else 
   begin PP1:=-7/48; QQ1:= 35/288 end;
   Y:= 3*PI/8; BB:= if A>=3 then A **(-2/3) else 0.0 ;
   for S:=1 step 1 until N do 
   begin if A=0and S=1and D=3 then 
       begin X:=0; J:=0 end else 
       begin if S >= 3*A -8  then 
           begin B:=(S+A/2-T)*PI; C:=1/B/B/64;
               X:=B-1/B/8*(P0-P1*C)/(1-Q1*C)
           end else 
           begin if S=1 then 
               begin X:= if D=1 then -2.33811 else 
                           if D=2 then -1.17371 else 
                           if D=3 then -1.01879 else -2.29444
               end else 
               begin X:= Y*(4*S-TT); V:= 1/X/X;
                   X:= -X**(2/3)*(1+V*(PP1+QQ1*V))
               end;
               U:=X*BB; V:=FI(2/3*(-U)**1.5);
               W:=1/COS(V); XX:=1-W*W; C:=SQRT(U/XX);
               X:=W*(A+C/A/U*
               (if D<3 then -5/48/U-C*(-5/24/XX+1/8)
                         else  7/48/U+C*(-7/24/XX+3/8)))
           end;  J:=0;
           L1: XX:=X*X; X4:=XX*XX; A2:=AA-XX; RO:=R; J:=J+1;
               if D<3 then 
               begin U:=RO; P:=(1-4*A2)/6/X/(2*A+1);
                   Q:=(2*(XX-MU)-1-6*A)/3/X/(2*A+1)
               end else 
               begin U:=-XX*RO/A2; V:=2*X*A2/(AA+XX)/3;
                   W:=A2*A2*A2;
                   Q:=V*(1+( MU2+32*MU*XX+48*X4)/32/W);
                   P:=V*(1+(-MU2+40*MU*XX+48*X4)/64/W)
               end;
           W:=U*(1+P*RO)/(1+Q*RO); X:=X+W;
           if ABS(W/X)>"-13 and J<5 then goto L1
       end; Z[S]:=X
   end 
end BESS ZEROS

        eop