NU123952@VM1.NoDak.EDU (Mark A. Ordal) (04/13/91)
Here's my latest set of Bessel function programs for the HP-48. The
only programs you need to call directly are Ynx and Jnx. However, if
you are plotting something (e.g., the intensity of the diffraction pattern
of a circular aperature -- which involves J1) you might get a slight
speed improvement by directly calling the J1 program (named ASJ1).
The comments are pretty complete so just read through the listing
to see how to call the programs and see where I got the algorithms.
Note that what follows is the ASCII (translate 3) listing of a directory.
The separate programs can be adapted for use on a 28S. In fact, most run
unchanged on the 28S. For Ynx and Jnx, you have to replace the CASE
statement with nested IF statements. The DECR command used in NRJN and
NRYN has to be replaced with 28S equivalents (the comments explain what to do).
Dr. Mark A. Ordal
Physics Department
North Dakota State University
Fargo, ND 58105
NU123952@NDSUVM1
%%HP: T(3)A(R)F(.);
@ --------------------------------------------------------------------
@
@ J (x) and Y (x)
@ n n
@
@ --------------------------------------------------------------------
@ Bessel Functions of integer order and the first and second kinds
@ --------------------------------------------------------------------
@ HP48 programs by Dr. Mark A. Ordal, North Dakota State University,
@ Physics Department, Fargo, ND, 58105
@ NU123952@NDSUVM1
@ revision of April 11, 1991
@ --------------------------------------------------------------------
@ This is the ASCII listing (translate 3 mode) of a directory that
@ contains the following:
@
@ Programs: Ynx, Jnx, ASY1, ASY0, ASJ1, ASJ0, NRYN, and NRJN
@ Subroutines: JY.SUB, JY1.SUB, and JY0.SUB
@
@ When named ASJY.NX, running the Bytes command on the name of this
@ directory returns: Checksum #36510 and 2651.5 bytes.
@ --------------------------------------------------------------------
DIR
Ynx
@ --------------------------------------------------------------------
@ Bessel Functions of the SECOND kind and integer order
@
@ To use:
@ Level 2 of stack: order (a nonnegative integer)
@ Level 1 of stack: argument (any positive real number)
@
@ Remember that Yn(x) is infinite for x=0 (These programs do not
@ prevent you from entering an invalid order or argument.)
@
@ When named Ynx, the BYTES command returns #18975 and 123.5 bytes
@ --------------------------------------------------------------------
\<< \-> n x
\<<
CASE n 0
SAME
THEN x
ASY0
END n 1
SAME
THEN x
ASY1
END n x
NRYN
END
\>>
\>>
Jnx
@ --------------------------------------------------------------------
@ Bessel Functions of the FIRST kind and integer order
@
@ To use:
@ Level 2 of stack: order (a nonnegative integer)
@ Level 1 of stack: argument (any nonnegative real number)
@
@ Remember: these programs do not prevent you from entering an
@ invalid order or argument.
@
@ When named Jnx, the BYTES command returns #51151 and 123.5 bytes
@ --------------------------------------------------------------------
\<< \-> n x
\<<
CASE n 0
SAME
THEN x
ASJ0
END n 1
SAME
THEN x
ASJ1
END n x
NRJN
END
\>>
\>>
ASY1
@ --------------------------------------------------------------------
@ Bessel Functions of the SECOND kind and order one calculated
@ using Abramowitz and Stegun Eqs 9.4.4 and 9.4.6
@
@ This program is slightly shorter and faster than the RPL equivalent
@ of the "Numerical Recipes" program BESSY1.
@ This program program calls programs JY.SUB, JY1.SUB, and ASJ1.
@
@ To use:
@ Level 1 of stack: argument (any positive real number)
@
@ When named ASY1, the BYTES command returns #36309 and 257.5 bytes
@ --------------------------------------------------------------------
\<< 0 \-> x a
\<<
IF x 3 <
THEN
-.6366198 .2212091
2.1682709
-1.3164827 .3123951
-.0400976 .0027873
x 3 / SQ 6 JY.SUB x
ASJ1 x .5 * LN * x
* 2 * \pi / \->NUM + x
/
ELSE x
JY1.SUB SIN * x \v/ /
END
\>>
\>>
ASY0
@ --------------------------------------------------------------------
@ Bessel Functions of the SECOND kind and order zero calculated
@ using Abramowitz and Stegun Eqs 9.4.1 and 9.4.3
@
@ This program is slightly shorter and faster than the RPL equivalent
@ of the "Numerical Recipes" program BESSY0.
@ This program program calls programs JY.SUB, JY0.SUB, and ASJ0.
@
@ To use:
@ Level 1 of stack: argument (any positive real number)
@
@ When named ASY0, the BYTES command returns #36663 and 243.5 bytes
@ --------------------------------------------------------------------
\<< 0 \-> x a
\<<
IF x 3 <
THEN
.36746691 .60559366
-.74350384
.25300117
-.04261214
.00427916
-.00024846 x 3 / SQ
6 JY.SUB x ASJ0 x
.5 * LN * 2 * \pi /
\->NUM +
ELSE x
JY0.SUB SIN * x \v/ /
END
\>>
\>>
ASJ1
@ --------------------------------------------------------------------
@ Bessel Functions of the FIRST kind and order one calculated
@ using Abramowitz and Stegun Eqs 9.4.4 and 9.4.6
@
@ This program is slightly shorter and faster than the RPL equivalent
@ of the "Numerical Recipes" program BESSJ1.
@ This program program calls programs JY.SUB and JY1.SUB.
@
@ To use:
@ Level 1 of stack: argument (any nonnegative real number)
@
@ When named ASJ1, the BYTES command returns #24497 and 201 bytes
@ --------------------------------------------------------------------
\<< 0 \-> x a
\<<
IF x 3 <
THEN .5
-.56249985
.21093573
-.03954289
.00443319
-.00031761
.00001109 x 3 / SQ
6 JY.SUB x *
ELSE x
JY1.SUB COS * x \v/ /
END
\>>
\>>
ASJ0
@ --------------------------------------------------------------------
@ Bessel Functions of the FIRST kind and order zero calculated
@ using Abramowitz and Stegun Eqs 9.4.1 and 9.4.3
@
@ This program is slightly shorter and faster than the RPL equivalent
@ of the "Numerical Recipes" program BESSJ0.
@ This program program calls programs JY.SUB and JY0.SUB.
@
@ To use:
@ Level 1 of stack: argument (any nonnegative real number)
@
@ When named ASJ0, the BYTES command returns #26811 and 186 bytes
@ --------------------------------------------------------------------
\<< 0 \-> x a
\<<
IF x 3 <
THEN 1
-2.2499997
1.2656208 -.3163866
.0444479 -.0039444
.00021 x 3 / SQ 6
JY.SUB
ELSE x
JY0.SUB COS * x \v/ /
END
\>>
\>>
NRYN
@ --------------------------------------------------------------------
@ Bessel Functions of the SECOND kind and integer order n > 1
@ calculated using the RPL equivalent of the "Numerical Recipes"
@ program BESSY. Needs the equivalent of the "Numerical Recipes"
@ programs BESSY0 and BESSY1. In this case, programs ASY0 and ASY1
@ are used instead -- based on Abramowitz and Stegun Eqs 9.4.1 and
@ 9.4.3 for order zero and Eqs 9.4.4 and 9.4.6 for order one.
@
@ The programs ASY0 and ASY1 based on the Abramowitz and Stegun
@ equations are shorter and faster than RPL equivalents of BESSY0
@ and BESSY1.)
@
@ To use:
@ Level 2 of stack: order (a positive integer > 1)
@ Level 1 of stack: argument (any positive real number)
@
@ When named NRYN, the BYTES command returns #8897 and 149 bytes
@ --------------------------------------------------------------------
\<< 0 \-> n x tox
\<< 2 x / 'tox'
STO x ASY0 x ASY1 1
n 1 -
FOR j j tox @ Keep 'bjm' value in top level of stack and
* OVER * 3 PICK - @ 'bj' value in next lower level. Variable
ROT DROP @ 'bjp' in FORTRAN program BESSJ is just used
NEXT SWAP @ to swap values between 'bjm' and 'bj'. (The
DROP @ old 'bj' value becomes the new 'bjm' value.)
\>>
\>>
NRJN
@ --------------------------------------------------------------------
@ Bessel Functions of the FIRST kind and integer order n > 1
@ calculated using the RPL equivalent of the "Numerical Recipes"
@ program BESSJ. Needs the equivalent of the "Numerical Recipes"
@ programs BESSJ0 and BESSJ1. In this case, programs ASJ0 and ASJ1
@ are used instead -- based on Abramowitz and Stegun Eqs 9.4.1 and
@ 9.4.3 for order zero and Eqs 9.4.4 and 9.4.6 for order one.
@
@ The programs ASJ0 and ASJ1 based on the Abramowitz and Stegun
@ equations are shorter and faster than RPL equivalents of BESSJ0
@ and BESSJ1.)
@
@ To use:
@ Level 2 of stack: order (a positive integer > 1)
@ Level 1 of stack: argument (any nonnegative real number)
@
@ When named NRJN, the BYTES command returns #20392 and 783.5 bytes
@
@ To adapt for use on the HP28S, replace the commands 'm' DECR with
@ the equivalent sequence m 1 - 'm' STO m
@ --------------------------------------------------------------------
\<< \-> n x
\<< 40
10000000000
.0000000001 0 1 0 0
0 0 0 \-> iacc bigno
bigni tox bj bjm
bjp m sum bessj
\<< 2 x /
'tox' STO
IF x n >
THEN x
ASJ0 x ASJ1 1 n 1 -
FOR j j @ Keep 'bjm' value in top level of stack and
tox * OVER * 3 PICK @ 'bj' value in next lower level. Variable
- ROT DROP @ 'bjp' in FORTRAN program BESSJ is just used
NEXT @ to swap values between 'bjm' and 'bj'. (The
SWAP DROP @ old 'bj' value becomes the new 'bjm' value.)
ELSE iacc
n * \v/ IP n + 2 / IP
2 * 'm' STO 0 1 m
DO tox @ Keep 'bjp' value in top level of stack and
* OVER * 3 PICK - @ 'bj' value in next lower level. Variable
ROT DROP @ 'bjm' in FORTRAN program BESSJ is just used
IF @ to swap values between 'bjp' and 'bj'. (The
DUP ABS bigno > @ old 'bj' value becomes the new 'bjp' value.)
THEN
bigni * SWAP bigni
* SWAP bigni DUP
'bessj' STO* 'sum'
STO*
END
IF m
n SAME
THEN
OVER 'bessj' STO
END
'm' DECR tox * OVER
* 3 PICK - ROT DROP
IF
DUP ABS bigno >
THEN
bigni * SWAP bigni
* SWAP bigni DUP
'bessj' STO* 'sum'
STO*
END
IF m
n SAME
THEN
OVER 'bessj' STO
END
DUP 'sum' STO+ 'm'
DECR
UNTIL m
1 <
END
DROP SWAP DROP NEG
sum 2 * + bessj
SWAP /
END
\>>
\>>
\>>
JY.SUB
@ --------------------------------------------------------------------
@ Subprogram for use by ASJ0, ASJ1, ASY0, and ASY1
@
@ When named JY.SUB, the BYTES command returns #10082 and 63 bytes
@ --------------------------------------------------------------------
\<< \-> t j
\<< 1 j
FOR k t * +
NEXT
\>>
\>>
JY1.SUB
@ --------------------------------------------------------------------
@ Subprogram for use by ASJ1 and ASY1
@
@ When named JY1.SUB, the BYTES command returns #4833 and 251.5 bytes
@ --------------------------------------------------------------------
\<< 0 \-> x a
\<< 3 x / 'a'
STO .79788456
.00000156 .01659667
.00017105
-.00249511
.00113653
-.00020033 a 6
JY.SUB -2.35619449
.12499612 .0000565
-.00637879
.00074348 .00079824
-.00029166 a 6
JY.SUB x + RAD
\>>
\>>
JY0.SUB
@ --------------------------------------------------------------------
@ Subprogram for use by ASJ0 and ASY0
@
@ When named JY0.SUB, the BYTES command returns #54944 and 251.5 bytes
@ --------------------------------------------------------------------
\<< 0 \-> x a
\<< 3 x / 'a'
STO .79788456
-.00000077
-.0055274
-.00009512
.00137237
-.00072805
.00014476 a 6
JY.SUB -.78539816
-.04166397
-.00003954
.00262573
-.00054125
-.00029333
.00013558 a 6
JY.SUB x + RAD
\>>
\>>
ENDakcs.joey@hpcvbbs.UUCP (joe montanino) (05/14/91)
Thanks for a program which is very useful to me. How about I(n,x) and K(n,x) for the modified bessel equation?