hp41programs

Bessel

Bessel Functions for the HP-41


Overview
 

 1°) Bessel & modified Bessel Functions of the first kind  Jn(x) &  In(x)
 2°) Bessel Functions Jn(x) of integer order ( n >= 0 )
 3°) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( non-integer order )
 4°) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( integer order )
 5°) Continued Fractions for Jn(x) & Yn(x)
 6°) An Asymptotic Expansion for  Jn(x) and Yn(x)
 7°) An Asymptotic Expansion for  Kn(x)
 8°) Integrals of  Jn(x) &  In(x)
 9°) Integral of  Jn(x)  ( integer order )
10°) Spherical Bessel Functions

   a) Spherical Bessel Functions of the first kind
   b) Spherical Bessel Functions of the second ( and first ) kind

11°) Zeros of Bessel Functions, Mac Mahon's Asymptotic Expansions

   a) Zeros of  Jn(x) & Yn(x)
   b) Zeros of  the first derivatives  J'n(x) & Y'n(x)

12°)  Bessel Functions - Complex order & argument

   a) Bessel Functions of the 1st kind
   b) Bessel Functions of the 2nd kind
   c) Bessel Functions of the 2nd kind ( real integer order )
 
 

1°) Bessel & modified Bessel Functions of the first kind  Jn(x) &  In(x)
 

Formulae:

       Jn(x) = (x/2)n [ 1/Gam(n+1)  +  (-x2/4)1/ (1! Gam(n+2) )  + .... + (-x2/4)k/ (k! Gam(n+k+1) ) + ....  ]      n # -1 ; -2 ; -3 ; ....

       In(x) = (x/2)n [ 1/Gam(n+1)  +  (x2/4)1/ (1! Gam(n+2) )  + .... + (x2/4)k/ (k! Gam(n+k+1) ) + ....  ]      n # -1 ; -2 ; -3 ; ....
 

Data Registers:  /
Flag:  F01
Subroutine:  "GAM" or "GAM+"  ( cf "Gamma Function for the HP-41" )

-Synthetic registers N & O may be replaced by any unused data registers.
 
 

01  LBL "JNX" 
02  2
03  /
04  STO N 
05  X<>Y
06  STO O        
07  1
08  ENTER^
09  STO M
10  ENTER^
11  LBL 01
12  X<> Z
13  FC? 01
14  CHS
15  RCL N         
16  X^2
17  *
18  RCL M
19  ST/ Y
20  RCL O        
21  +
22  /
23  ISG M 
24  CLX
25  ST+ Y
26  X<> Z
27  X#Y?
28  GTO 01
29  RCL N
30  RCL O
31  Y^X
32  *
33  X<> O        
34  1
35  +
36  XEQ "GAM"
37  RCL O        
38  X<>Y
39  /
40  CLA
41  END

 
    ( 70 bytes / SIZE 000 )
 
 

      STACK        INPUTS     OUTPUTS
           Y             n          f(x)
           X             x          f(x)

  where  f(x) = Jn(x)  if flag F01 is clear  and  f(x) = In(x)  if flag F01 is set.

Example:    Calculate   J0.7(1.9)  and  I0.7(1.9)

   CF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   J0.7(1.9) = 0.584978102   (  in 12 seconds )

   SF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   I0.7(1.9) = 1.727630598   (  in 12 seconds )

Notes:

-If  n  is a positive integer,  J-n(x) = (-1)n Jn(x)   and    I-n(x) = In(x)
 J0(0) = 1 but this program doesn't work if x = n = 0 ( DATA ERROR line 31 )
-Add  X=Y?  X#0?  GTO 00  SIGN  X<>Y  LBL 00  after line 30 and it will work in this case too.

-If you have the M-Code routine "HGF+" in your HP-41 ( cf "Hypergeometric Functions" ), the routine may be simplified as follows:
-Line 21 is - if possible - an M-Code routine that gives  0^0 = 1, see for instance "A few M-Code Routines for the HP-41"
 
 

 01  LBL "JNX"
 02  X<>Y
 03  STO 02
 04  1
 05  +
 06  STO 01
 07  CLX
 08  2
 09  /
 10  STO 00
 11  CLST
 12  SIGN
 13  CHS
 14  RCL 00
 15  X^2
 16  FC? 01
 17  CHS
 18  HGF+
 19  RCL 00
 20  RCL 02 
 21  Y^X
 22  *
 23  END

 
( 33 or 34 bytes / SIZE 003 )
 
 

      STACK        INPUTS     OUTPUTS
           Y             n            /
           X             x          f(x)

  where  f(x) = Jn(x)  if flag F01 is clear  and  f(x) = In(x)  if flag F01 is set.

Example:

   CF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   J0.7(1.9) = 0.584978103   (  in 4 seconds )

   SF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   I0.7(1.9) = 1.727630603   (  in 4 seconds )
 

-However, with or without "HGF+",  Jn(x) is not accurately computed for large arguments, "JNX1" is better:
 

2°) Bessel Functions Jn(x) of integer order ( n >= 0 )
 

-The following program is only a modification of a program given by Keith Jarett in "HP-41 Extended Functions made easy"
-It uses the relations:   J0(x) + 2 J2(x) + 2 J4(x) + 2 J6(x) + ...... = 1  and   Jn-1(x) + Jn+1(x) = (2n/x) Jn(x)
 

Data Registers:  R00 and R02 to R05: temp   R01 = x
Flags: /
Subroutines: /
 
 

01  LBL "JNX1"
02  X#0?
03  GTO 00      
04  X#Y?
05  RTN
06  SIGN
07  RTN
08  LBL 00
09  STO 01
10  ABS
11  5
12  +
13  X<>Y
14  STO 00      
15  X<Y?
16  X<>Y
17  INT
18  ST+ X
19  STO 03
20  ST- 00
21  CLST
22  STO 04
23  STO 05      
24  SIGN
25  LBL 01
26  STO Z
27  RCL 03 
28  ST+ X
29  *
30  RCL 01
31  /
32  X<>Y
33  - 
34  ISG 00
35  STO 02      
36  ST+ 04
37  RCL 05 
38  X<> 04
39  STO 05
40  RDN
41  DSE 03
42  GTO 01
43  RCL 05      
44  ST+ X
45  X<>Y
46  -
47  RCL 02 
48  X<>Y
49  /
50  END

 
    ( 70 bytes / SIZE 006 )
 
 

      STACK        INPUTS     OUTPUTS
           Y             n           /
           X             x         Jn(x)

 
Examples:

     2    ENTER^
    10   XEQ "JNX1"   produces   J2(10) =  0.254630314  (  in 17 seconds )

     3    ENTER^
   100    R/S   >>>>    J3(100)  =  7.628420178 10 -2  ( in 2mn01s )

-Unlike "JNX" , "JNX1" yields accurate results for large x-values

-If you replace lines 45 thru 49 by  RCL Y   -   ST/ 02   ST/ Z   /   RCL 02
 you'll get  J1(x)  in Z-register , J0(x)  in Y-register and  Jn(x)  in X-register.
 

3°) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( non-integer order )
 

Formulae:     Yn(x) = ( Jn(x) cos(n(pi)) - J-n(x) ) / sin(n(pi))      ;      Kn(x) = (pi/2) ( I-n(x) - In(x) ) / sin(n(pi))      n # .... -3 ; -2 ; -1 ; 0 ; 1 ; 2 ; 3 ....
 

Data Registers:   R00 = x , R01 = n , R02 = temp

Flag: F01      If  F01 is clear,  "YNX" calculates  Yn(x)
                      If  F01 is set,     "YNX" ---------   Kn(x)
Subroutine: "JNX"
 
 

01  LBL "YNX"
02  DEG
03  STO 00      
04  X<>Y
05  STO 01
06  CHS
07  X<>Y
08  XEQ "JNX"
09  STO 02
10  RCL 01      
11  RCL 00
12  XEQ "JNX"
13  PI
14  R-D
15  RCL 01      
16  *
17  STO Z
18  FS? 01
19  CLX
20  COS
21  *
22  RCL 02      
23  -
24  X<>Y
25  SIN
26  /
27  PI
28  2
29  /
30  CHS
31  FC? 01       
32  ST/ X
33  *
34  END

 
    ( 54 bytes / SIZE 003 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x           f(x)

  with  f(x) =  Yn(x)  if  CF 01  ;   f(x) =  Kn(x)  if SF 01

Example:  Evaluate  Y1.4(3)  and   K1.4(3)

    CF 01   1.4  ENTER^   3   XEQ "YNX"   >>>>   Y1.4(3) = 0.137821836    ( in 31 seconds )
    SF 01    1.4  ENTER^   3          R/S          >>>>   K1.4(3) = 0.046088036      --------------
 

4°) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( integer order )
 

Formulae:       with  psi(x) = the digamma function, we have ( positive n ):

     Yn(x) = -(1/pi) (x/2)-n SUMk=0,1,....,n-1  (n-k-1)!/(k!) (x2/4)k + (2/pi) ln(x/2) Jn(x) - (1/pi) (x/2)n SUMk=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (-x2/4)k / (k!(n+k)!)

Kn(x) =  (1/2) (x/2)-n SUMk=0,1,....,n-1  (n-k-1)!/(k!) (-x2/4)k - (-1)n ln(x/2) In(x) + (1/2) (-1)n (x/2)n SUMk=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (x2/4)k / (k!(n+k)!)
 
 

Data Registers:   R00 = x/2 ; R01 = n ; R02 to R05: temp

Flag: F01      If  F01 is clear,  "YNX1" calculates  Yn(x)
                      If  F01 is set,     "YNX1" ---------   Kn(x)

Subroutine: "JNX"
 
 

  01  LBL "YNX1"
  02  STO 00        
  03  2
  04  ST/ 00
  05  X<> Z
  06  STO 01
  07  X<>Y
  08  XEQ "JNX"
  09  RCL 00
  10  X^2
  11  LN
  12  *
  13  1
  14  FS? 01
  15  CHS
  16  RCL 01
  17  Y^X
  18  *
  19  STO 02
  20  RCL 01
  21  STO 03
  22  CLX
  23  LBL 01
  24  RCL 01        
  25  RCL 03
  26  -
  27  FACT
  28  RCL 03 
  29  1
  30  -
  31  X<0?
  32  CLST
  33  FACT
  34  ST/ Y
  35  CLX
  36  RCL 00
  37  ST* X
  38  FS? 01
  39  CHS
  40  LASTX
  41  Y^X
  42  *
  43  +
  44  DSE 03
  45  GTO 01
  46  RCL 00
  47  RCL 01
  48  Y^X
  49  /
  50  ST- 02
  51  1.15443133 
  52  STO 05        
  53  RCL 01
  54  LBL 02
  55  X#0?
  56  1/X
  57  ST- 05
  58  X<> L
  59  DSE X
  60  GTO 02
  61  1
  62  STO 03
  63  STO 04
  64  RCL 05
  65  LBL 03
  66  RCL 04
  67  RCL 00
  68  X^2
  69  FC? 01
  70  CHS
  71  *
  72  RCL 03 
  73  ST/ Y
  74  RCL 01        
  75  +
  76  /
  77  STO 04
  78  RCL 05
  79  LASTX
  80  1/X
  81  -
  82  RCL 03
  83  1/X
  84  -
  85  STO 05
  86  *
  87  +
  88  ISG 03
  89  CLX
  90  X#Y?
  91  GTO 03
  92  RCL 01 
  93  FACT
  94  /
  95  RCL 00        
  96  FS? 01
  97  CHS
  98  RCL 01
  99  Y^X
100  *
101  RCL 02
102  +
103  FC? 01
104  PI
105  FS? 01
106  -2
107  /
108  END

 
    ( 151 bytes / SIZE 006 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x           f(x)

  with  f(x) =  Yn(x)  if  CF 01  ;   f(x) =  Kn(x)  if SF 01   ;   n = 0 , 1 , 2 , ........ , 69

Example:  Evaluate  Y2(3)  and   K2(3)

    CF 01    2  ENTER^   3   XEQ "YNX1"   >>>>   Y2(3) = -0.160400393    ( in 23 seconds )
    SF 01    2  ENTER^   3          R/S            >>>>    K2(3) =  0.061510458      --------------

Note:

-If n is a positive integer, Y-n(x) = (-1)n Yn(x)  and  K-n(x) = Kn(x)
 

5°) Continued Fractions for Jn(x) & Yn(x)
 

-Unlike   "JNX" &  "YNX"  ( with CF 01 )  and  "YNX1"  , the following program produces accurate results for large arguments.

Formulae:

 With   p + i.q = -1/(2x) + i + (i/x) [ ( 0.52 - n2 )/( 2x + 2i + ( 1.52 - n2 )/( 2x + 4i + ..... ) ) ]
  and      gn =  -1/(((2n + 2)/x) - 1/(((2n + 4)/x) - ..... ))

   Jn(x)  =  sign(D) [ ( q2 + ( p - gn - n/x )2 ) ( 2q/(x.Pi) ) ] -1/2    where  D = the denominator of the second continued fraction
   Yn(x) =  [ ( p - gn - n/x )/q ] Jn(x)
 

Data Registers:   R00 thru R15 are used by "CFRZ"  R01 = x ,  R16 = n
Flags: /
Subroutines: "CFR" & "CFRZ"  ( cf "Continued Fractions for the HP-41" )
 
 

01  LBL "JYNX"
02  STO 01
03  X<>Y
04  STO 16
05  "T"
06  ASTO 00
07  CLST
08  RCL 01
09  XEQ "CFRZ"
10  RCL 01        
11  ST/ Z
12  /
13  1
14  +
15  X<>Y
16  CHS
17  RCL 01
18  ST+ X
19  1/X
20  -
21  STO 09
22  X<>Y
23  STO 10
24  "U"
25  ASTO 00
26  CLX
27  CLA
28  RCL 01        
29  XEQ "CFR"
30  CHS
31  STO 08
32  RCL 09
33  +
34  RCL 16
35  RCL 01
36  /
37  -
38  STO 11
39  RCL 10 
40  R-P
41  LAST X
42  ST+ X
43  PI
44  RCL 01        
45  *
46  /
47  SQRT
48  X<>Y
49  /
50  RCL 04
51  SIGN
52  *
53  STO 12
54  RCL 11
55  *
56  RCL 10 
57  /
58  RCL 12
59  RTN
60  LBL "T"
61  RCL 03        
62  ST+ X
63  RCL 01
64  ST+ X
65  RCL 03
66  .5
67  -
68  X^2
69  RCL 16 
70  X^2
71  -
72  0
73  X<>Y
74  RTN
75  LBL "U"
76  RCL 02
77  RCL 16        
78  +
79  ST+ X
80  RCL 01
81  /
82  1
83  CHS
84  END

 
    ( 125 bytes / SIZE 017 )
 
 

      STACK        INPUTS      OUTPUTS
           Y         n >= 0          Yn(x)
           X          x > 0          Jn(x)

 
Examples:

    10   ENTER^  XEQ "JYNX"   >>>>   J10(10) =  0.207486107   X<>Y   Y10(10) =  -0.359814151    ( in 2mn27s )

  3.14  ENTER^
   100  XEQ "JYNX"  >>>>  J3.14(100) =  0.079535723     X<>Y    Y3.14(100) =  0.006582327    ( in 4mn14s )

-The method doesn't work if n is a negative integer.
-If n < 0 we can use the relations  Jn = J-n cos n.Pi + Y-n sin n.Pi  and  Yn = -J-n sin n.Pi + Y-n cos n.Pi
-If x < 0 the results are generally complex.
 

6°) An Asymptotic Expansion for  Jn(x) and Yn(x)
 

Formulae:

    Jn(x) = (2/(pi*x))1/2 ( P.cos(x-(2n+1)pi/4) - Q.sin(x-(2n+1)pi/4) )  and   Yn(x) = (2/(pi*x))1/2 ( P.sin(x-(2n+1)pi/4) + Q.cos(x-(2n+1)pi/4) )
 

  where   P ~  1 - (4n2-1)(4n2-9)/(2!(8x)2) + (4n2-1)(4n2-9)(4n2-25)(4n2-49)/(4!(8x)4) - ......

     and   Q ~  (4n2-1)/(8x) - (4n2-1)(4n2-9)(4n2-25)/(3!(8x)3) + ......
 

Data Registers:  R00 = x ; R01 = n ; R02 to R06: temp
Flags: /
Subroutines: /
 
 

01  LBL "AEJY"
02  DEG
03  STO 00       
04  X<>Y
05  STO 01
06  ST+ X
07  X^2
08  STO 04
09  1
10  STO 02
11  -
12  RCL 00
13  8
14  *
15  STO 05
16  ST* 05
17  /
18  STO 03
19  XEQ 01
20  STO 06
21  CLX
22  STO 02 
23  SIGN
24  STO 03       
25  XEQ 01
26  GTO 02
27  LBL 01
28  RCL 03
29  RCL 02
30  2
31  +
32  STO 02
33  ST+ X
34  3
35  -
36  X^2
37  RCL 04
38  -
39  *
40  RCL 04 
41  RCL 02       
42  ST+ X
43  DSE X
44  X^2
45  -
46  *
47  RCL 02
48  ST/ Y
49  DSE X
50  /
51  RCL 05
52  /
53  STO 03
54  +
55  X#Y?
56  GTO 01
57  RTN
58  LBL 02 
59  RCL 00       
60  RCL 01
61  ST+ X
62  1
63  +
64  4
65  /
66  PI
67  *
68  -
69  R-D
70  STO 04
71  X<>Y
72  P-R
73  RCL 04
74  RCL 06       
75  P-R
76  ST+ T
77  RDN
78  -
79  2
80  PI
81  /
82  RCL 00
83  /
84  SQRT
85  ST* Z
86  *
87  END

 
    ( 112 bytes / SIZE 007 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             n         Yn(x)
           X          x > 0          Jn(x)

 
Example:       4   ENTER^  100   XEQ "AEJY"  >>>>   J4(100)  =  0.026105809
                                                             X<>Y                Y4(100)  = -0.075430120    ( in 10 seconds )
 

-The program will not work if  x is too "small" or if n is too "large"
-If you have the M-Code routine HGF+ in your HP-41, you can use the following variant:
 
 

01  LBL "AEJY"
02  DEG
03  STO 00       
04  X<>Y
05  ST+ X
06  STO 06
07  1
08  +
09  STO 01
10  LASTX
11  RCL 06
12  -
13  STO 02
14  LASTX
15  3
16  +
17  STO 03
18  LASTX
19  RCL 06       
20  -
21  STO 04
22  .5
23  STO 05
24  X^2
25  ST* 01
26  ST* 02
27  ST* 03
28  ST* 04
29  1/X
30  1
31  RCL 00
32  X^2
33  1/X
34  CHS
35  STO 06       
36  HGF+
37  STO 07
38  RCL 01
39  STO 08
40  RCL 02
41  STO 09
42  1
43  ST+ 01
44  ST+ 02
45  ST+ 05
46  4
47  X<>Y
48  RCL 06
49  HGF+
50  RCL 08 
51  RCL 09       
52  *
53  ST+ X
54  *
55  RCL 00
56  ST/ Y
57  PI
58  RCL 08
59  *
60  -
61  R-D
62  STO 06
63  X<>Y
64  P-R
65  RCL 06
66  RCL 07       
67  P-R
68  X<> Z
69  -
70  X<> Z
71  +
72  PI
73  RCL 00
74  *
75  2
76  /
77  SQRT
78  ST/ Z
79  /
80  END

 
    ( 105 bytes / SIZE 010 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             n         Yn(x)
           X           x > 0          Jn(x)

 
Examples:

    4     ENTER^
  100   XEQ "AEJY"  >>>>      J4(100)  =  0.026105809
                                 X<>Y     Y4(100)  = -0.075430120    ( in 6 seconds )

-Likewise, we find

     JPI(11.6)  =  0.238578119
    YPI(11.6)  =  0.002890137      ( in  13 seconds )

-But for  n = PI  and  x = 11.5 , the series diverge too soon and we have an infinite loop ( press any key to stop the loop - or  ENTER^  ON )
 

7°) An Asymptotic Expansion for  Kn(x)
 

Formula:     Kn(x) ~  (pi/(2x))1/2 e-x  ( 1 + (4n2-1)/(8x) + (4n2-1)(4n2-9)/(2!(8x)2) + (4n2-1)(4n2-9)(4n2-25)/(3!(8x)3) + ......  )
 

Data Registers:  R00 = x ; R01 = 4n2 ; R02 , R03: temp
Flags: /
Subroutines: /
 
 

01  LBL "AEK"
02  STO 00     
03  X<>Y
04  ST+ X
05  X^2
06  STO 01
07  SIGN
08  STO 02
09  STO 03
10  LBL 01
11  RCL 03     
12  RCL 01
13  RCL 02
14  ST+ X
15  DSE X
16  X^2
17  -
18  *
19  RCL 00 
20  RCL 02      
21  *
22  8
23  *
24  /
25  STO 03
26  +
27  ISG 02
28  CLX
29  X#Y?
30  GTO 01      
31  PI
32  RCL 00
33  ST+ X
34  /
35  SQRT
36  *
37  RCL 00      
38  E^X
39  /
40  END

 
    ( 54 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x          Kn(x)

 
Examples:

    2   ENTER^  10   XEQ "AEK"  >>>>   K2(10)  = 0.000,021,509,817,05    ( in 15 seconds )
  1.4  ENTER^  19        R/S           >>>>   K1.4(19) = 1.683198846 10-9          ( in 6 seconds )

-The program will not work if  x is too "small" or if n is too "large".
-For instance, n = 2 , x = 7 dont work.
-However, if we stop the iterations before the terms begin to increase, we find  K2(7) = 0.0005545622
-Divergence is more typical than convergence for asymptotic series...

Remark:     In(x) is accurately computed by "JNX" with flag F01 set, but if you need an asymptotic expansion for this function:

  replace line 39 with *
              line 36 with /
              line 34 with *  and add CHS after line 24   ( for instance,   I1.4(19) = 15,597,340  ( in 8 seconds ) )

-The M-Code routine HGF+ may also be used to simplify this program ( cf "Hypergeometric Functions for the HP-41" )
 
 

 01  LBL "AEK"
 02  STO 00
 03  X<>Y
 04  STO 01
 05  CHS
 06  .5
 07  ST+ 01
 08  +
 09  STO 02
 10  LASTX
 11  1/X
 12  0
 13  RCL 00
 14  ST+ X
 15  STO 03
 16  1/X
 17  CHS
 18  HGF+
 19  PI
 20  RCL 03
 21  /
 22  SQRT
 23  *
 24  RCL 00
 25  CHS
 26  E^X
 27  *
 28  END

 
   ( 40 bytes )
 
 

      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x          Kn(x)

 
Example:

     PI   10.1    XEQ "AEK"   >>>>   KPi(10.1) = 0.00002545492107       ---Execution time = 6s---

-If  x is not large enough, the divergence occurs too soon and we fall in an infinite loop
-Press any key ( or ENTER^  ON ) to stop the loop.

-For instance, with n = PI and x = 10 , this program doesn't work.
 

8°) Integrals of   Jn(x)  &  In(x)
 

-Here, we integrate the series expansions used by "JNX"     ( n # -1 ; -2 ; -3 ; .... )
 

Data Registers:  /

Flag:  F01      CF 01  for  §0x Jn(t).dt
                       SF 01 for  §0x In(t).dt

Subroutine:  "1/G+"  ( Reciprocal of the Gamma Function )

-Synthetic registers M , N & O may be replaced by any unused data registers.
-Lines 46-47-48 may be replaced by  XEQ "GAM+"  RCL O  X<>Y  /
 
 

01  LBL "ITIJ" 
02  STO N     
03  X<>Y
04  STO O 
05  2
06  ST/ N
07  SIGN
08  STO M
09  +
10  /
11  ENTER^
12  ENTER^
13  LBL 01     
14  X<> Z
15  FC? 01
16  CHS
17  RCL N 
18  X^2
19  *
20  RCL M
21  ST/ Y
22  RCL O
23  +
24  ST/ Y
25  RCL M     
26  +
27  1
28  -
29  ST* Y
30  2
31  +
32  /
33  ISG M
34  CLX
35  ST+ Y
36  X<> Z 
37  X#Y?
38  GTO 01     
39  RCL N
40  RCL O
41  Y^X
42  *
43  X<> O
44  1
45  +
46  XEQ "1/G+"
47  RCL O       
48  *
49  CLA
50  END

 
    ( 85 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             x       §0x fn(t).dt 

  where   fn(t) = Jn(t)  if flag F01 is clear
    and    fn(t) =  In(t)  if flag F01 is set

Examples:

          CF 01
  1.4   ENTER^
   3     XEQ"ITIJ"  >>>>   §03  J1.4(x).dx =  1.049262785  ( in 14s )

          SF 01
  1.4   ENTER^
   3     XEQ"ITIJ"  >>>>   §03  I1.4(x).dx =  2.918753200  ( in 13s )

-If flag F01 is clear, this program cannot be used for large arguments ( like "JNX" )
 

9°) Integral of   Jn(x) ( integer order )
 

Formula:       §0x Jn(t).dt = 2 ( Jn+1(x) + Jn+3(x) + Jn+5(x) + ........ )
 

Data Registers:   R00 thru R05 ( or R06 ) are used by "JNX1" ( R01 = x )  R07: temp R08: partial sums and eventually,  (1/2) §0x Jn(t).dt
Flag:  F07
Subroutine: "JNX1"
 
 

01  LBL "ITJ"
02  STO 01       
03  X<>Y
04  STO 07 
05  CLX
06  STO 08
07  DSE 07
08  LBL 01
09  SF 07
10  LBL 02 
11  RCL 07       
12  2
13  +
14  STO 07 
15  RCL 01
16  XEQ "JNX1"
17  RCL 08       
18  +
19  STO 08 
20  LASTX
21  X#Y?
22  GTO 01       
23  FS?C 07
24  GTO 02
25  ST+ X        
26  END

 
    ( 45 bytes / SIZE 009 )
 
 

      STACK        INPUTS      OUTPUTS
           Y              n             /
           X             x        §0x Jn(t).dt 

 
Examples:

     1    ENTER^
     3    XEQ "ITJ"  >>>>  §03  J1(x).dx = 1.260051955  ( in 2mn01s )                 - "ITIJ" yields the same result in 14 seconds -

     0    ENTER^
    10      R/S         >>>>  §010  J0(x).dx = 1.067011304  ( in 5mn54s )                 - "ITIJ" yields 1.067011425 ( error ~ 10 -7 )  in 23 seconds -

    50    ENTER^
    30      R/S         >>>>  §030  J50(x).dx = 1.478729947 10 -8  ( in 12mn29s )

-For small arguments, "ITJ" is much slower than "ITIJ"
-There is also a risk of  OUT OF RANGE  for large x-values,
  but if you modify "JNX1" as explained on the right of its listing ( cf paragraph 2°) )
  you'll get for instance §0100 J50(x).dx = 1.088806747 ( in about 10 seconds with a V41 emulator in turbo mode )
 

10°) Spherical Bessel Functions
 

     a) Spherical Bessel Functions of the first kind
 

-Spherical Bessel functions of the first kind are defined by   jn(x) = (Pi/(2x))1/2 Jn+1/2(x)   where n is an integer
-So, we can compute these functions by one of the programs above.
-However, another approach is used hereafter:
 

Formulae:             jn-1(x) =  jn(x) (2n+1)/x - jn+1(x)   this recurrence relation is used , starting with  jm(x) = 1 , jm+1(x) = 0  for some large m
 -Then, the results are normalized by the sum     1 j02(x) + 3 j12(x) + 5 j22(x) + ....... = 1
 -We could also normalize by calculating  j0(x) = (sin x)/x  but  this would be inaccurate if  sin x ~ 0
 

Data Registers:   R00 , R02 , R03: temp    R01 = x , R04 = Sum
Flags: /
Subroutines: /

-The sum in register R04 may easily exceed  9.999999999 1099
-So, in order to avoid an OUT OF RANGE,
 add    ENTER^   ABS   RCL 05   X>Y?   GTO 01   ST/ 02   ST/ Z   ST/ T   X^2   ST/ 04   LBL 01   R^  R^   after line 36
 and add   E40   STO 05  after line 23
-Thus you'll get, for instance,   j100(50) = 1.019012264 10 -22  ( in 3mn22s )
 
 

01  LBL "SB1"
02  X#0?
03  GTO 00   
04  X#Y?
05  RTN
06  SIGN
07  RTN
08  LBL 00
09  STO 01
10  ABS
11  5
12  +
13  X<>Y
14  STO 00   
15  X<Y?
16  X<>Y
17  INT
18  ST+ X
19  STO 03
20  ST- 00
21  ST+ X
22  STO 04
23  CLST
24  SIGN
25  ST+ 04
26  LBL 01
27  RCL 03   
28  ST+ X
29  1
30  +
31  RCL Y
32  *
33  RCL 01
34  /
35  R^
36  - 
37  ISG 00 
38  STO 02   
39  ENTER^
40  X^2 
41  RCL 03
42  ST+ X
43  DSE X
44  *
45  ST+ 04
46  RDN
47  DSE 03
48  GTO 01
49  RCL 02   
50  RCL 04
51  SQRT
52  /
53  END

 
    ( 74 bytes / SIZE 005 )
 
 

      STACK        INPUTS      OUTPUTS
           Y          n >= 0             /
           X             x           jn(x) 

 
Examples:

  2  PI  XEQ "SB1"  >>>>  j2(Pi) =  0.303963551   ( execution time = 14 seconds )

  10  ENTER^
   2      R/S        >>>>  j10(2) =  6.825300865 10 -8  ( in 17 seconds )

  100  ENTER^  R/S   >>>>  j100(100) =  0.01088047702  ( in 3mn )
 

     b) Spherical Bessel Functions of the second ( and first ) kind
 

-Spherical Bessel functions of the second kind are defined by    yn(x) = (Pi/(2x))1/2 Yn+1/2(x)   where n is an integer
-But the following program uses another method:
 

Formulae:

      yn+1(x) =  yn(x) (2n+1)/x - yn-1(x)                                  ,        jn+1(x) =  jn(x) (2n+1)/x - jn-1(x)
       y0(x) = -(cos x)/x  ,  y1(x) = -(cos x)/x2 - (sin x)/x          ,        j0(x) = (sin x)/x  ,   j1(x) = (sin x)/x2 - (cos x)/x
 

Data Registers:   R00  = k.nnn  ; R01 = x
Flag: F01    CF 01 to calculate  yn(x)
                    SF 01 to calculate   jn(x)    Set flag F01 only if  x is not significantly smaller than  n  -  the reccurence relation is unstable otherwise.
Subroutines: /

-Line 21 may be replaced by RTN and line 33 may be deleted.
 
 

01  LBL "SB2"
02  RAD
03  STO 00   
04  X<>Y
05   E3
06  /
07  STO 01
08  SIGN
09  P-R
10  RCL 00   
11  ST/ Z
12  /
13  FC? 01
14  CHS
15  DEG
16  FS? 01
17  X<>Y
18  LBL 01   
19  ISG 01
20  FS? 30
21  GTO 02
22  STO Z
23  RCL 01 
24  INT
25  ST+ X
26  DSE X
27  *
28  RCL 00   
29  /
30  X<>Y
31  -
32  GTO 01   
33  LBL 02
34  END

 
   ( 53 bytes / SIZE 002 )
 
 

      STACK        INPUTS      OUTPUTS
           Y   0 <= n <  1000             /
           X         x # 0           fn(x) 

 where  fn(x) = yn(x)  if  CF 01
   and   fn(x) = jn(x)   if  SF 01     ( only if  x is not significantly smaller than  n )

Examples:
 

a)         CF 01
     2     ENTER^
   3.14  XEQ "SB2"  >>>>   y2(3.14) = -0.222053752  ( in 2 seconds )

           CF 01
    30   ENTER^
     5   XEQ "SB2"  >>>>   y5(30) = -7.760717577 1018  ( in 15 seconds )

b)         SF 01
      4    ENTER^
    100  XEQ "SB2"  >>>>   j4(100) = -4.179461836 10 -3  ( in 3 seconds )

            SF 01
    100  ENTER^   XEQ "SB2"  >>>>   j100(100) = 1.088047702 10 -2  ( in 47 seconds )
 

-If  n < 0 , we can use the relation    yn(x) = (-1)n+1 j -n-1(x)   which actually holds for all integers n  ( n < 0 , n = 0 , n > 0 )
 

11°) Zeros of Bessel Functions:  Mac Mahon's Asymptotic Expansions
 

   a) Zeros of  Jn(x) & Yn(x)
 

-The s-th positive zero of  Jn(x) is denoted  jn,s and the s-th positive zero of  Yn(x) is denoted  yn,s

-This routine employs the first 5 terms of Mac Mahon's asymptotic expansions ( s >> n ):

    jn,s ~  b - (µ-1) / (8.b) - (4/3) (µ-1) (7µ-31) / (8.b)3 - (32/15) (µ-1) (83µ2-982µ+3779) / (8.b)5
                                      - (64/105) (µ-1) (6949µ3-153855µ2+1585743µ-6277237) / (8.b)7 - ................

       where   b = ( s + n/2 - 1/4 ) PI   and   µ = 4 n2
 

   yn,s ~  the same formula  with  b' = ( s + n/2 - 3/4 ) PI  instead of b
 

Data Registers:    R00  temp

                               R01 = b         R03 = 8.b       R05 = (8.b)2      R07 = last term of the series for J
                               R02 = b'        R04 = 8.b'      R06 = (8.b')2      R08 = last term of the series for Y
Flags: /
Subroutines: /
 
 

  01  LBL "JYZ"
  02  X<>Y
  03  STO 02    
  04  2
  05  /
  06  +
  07  4
  08  1/X
  09  -
  10  STO 01
  11  .5
  12  -
  13  PI
  14  ST* 01
  15  *
  16  X<> 02
  17  ST+ X
  18  X^2
  19  ENTER^
  20  ENTER^
  21  STO 00
  22  6949
  23  *
  24  153855
  25  -
  26  *
  27  1585743
  28  +
  29  *
  30  6277237
  31  -
  32  64
  33  *
  34  105
  35  /
  36  STO 07    
  37  STO 08 
  38  RCL 01
  39  RCL 02
  40  8
  41  ST* Z
  42  *
  43  STO 04
  44  X^2
  45  STO 06
  46  ST/ T
  47  RDN
  48  STO 03
  49  X^2
  50  STO 05
  51  /
  52  RCL 00 
  53  83
  54  *
  55  982
  56  -
  57  RCL 00    
  58  *
  59  3779
  60  +
  61  32
  62  *
  63  15
  64  /
  65  ST+ Z
  66  +
  67  RCL 05
  68  /
  69  X<>Y
  70  RCL 06
  71  /
  72  RCL 00 
  73  7
  74  *
  75  31
  76  -
  77  .75
  78  /
  79  ST+ Z
  80  +
  81  RCL 06    
  82  /
  83  X<>Y
  84  RCL 05
  85  /
  86  1
  87  ST+ Y
  88  ST+ Z
  89  RCL 00
  90  -
  91  ST* 07
  92  ST* 08
  93  ST* Z
  94  *
  95  RCL 03 
  96  /
  97  X<>Y
  98  RCL 04
  99  /
100  RCL 03    
101  7
102  Y^X
103  ST/ 07
104  CLX
105  RCL 04
106  LASTX
107  Y^X
108  ST/ 08
109  CLX
110  RCL 02
111  +
112  X<>Y
113  RCL 01
114  +
115  END

 
    ( 172 bytes / SIZE 009 )
 
 

      STACK        INPUTS      OUTPUTS
           Y          n >= 0          yn,s
           X              s           jn,s

 
Examples:

  PI   ENTER^
 12    XEQ "JYZ"   >>>>   jPi,12 = 41.73324496            R07 = -12 E-9
                              RDN   yPi,12 = 40.15792512            R08 = -16 E-9

  0    ENTER^
  4       R/S            >>>>    j0,4  = 11.79153443            R07 =  -58 E-9
                              RDN   y0,4  = 10.22234503            R08 = -158 E-9

  0    ENTER^
  3       R/S            >>>>    j0,3  = 8.653728            R07 =  -508 E-9
                              RDN   y0,3  = 7.086051            R08 = -2069 E-9

Notes:

-Registers R07 & R08 contain the 5th terms of the expansions: this may be used to give an idea of the accuracy.
-In the first 2 examples, almost all the decimals are correct. In the 3rd example, the results are exact to 6D
-Otherwise, you can delete lines 100 to 109 and lines  92-91-37-36
 

   b) Zeros of  the first derivatives  J'n(x) & Y'n(x)
 

-The s-th positive zero of  J'n(x) is denoted  j'n,s and the s-th positive zero of  Y'n(x) is denoted  y'n,s

-This routine also uses the first 5 terms of Mac Mahon's asymptotic expansions ( s >> n ):

    j'n,s ~  b - (µ+3) / (8.b) - (4/3) (7µ2+82µ-9) / (8.b)3 - (32/15) (83µ3+2075µ2-3039µ+3537) / (8.b)5
                                         - (64/105) (6949µ4+296492µ3-1248002µ2+7414380µ-5853627) / (8.b)7 - ................

       where   b = ( s + n/2 - 3/4 ) PI   and   µ = 4 n2
 

   y'n,s ~  the same formula  with  b' = ( s + n/2 - 1/4 ) PI  instead of b
 

Data Registers:    R00  temp

                               R01 = b         R03 = 8.b       R05 = (8.b)2      R07 = last term of the series for J'
                               R02 = b'        R04 = 8.b'      R06 = (8.b')2      R08 = last term of the series for Y'
Flags: /
Subroutines: /
 
 

  01  LBL "dJYZ"
  02  X<>Y
  03  STO 01     
  04  2
  05  /
  06  +
  07  4
  08  1/X
  09  -
  10  STO 02
  11  .5
  12  -
  13  PI
  14  ST* 02
  15  *
  16  X<> 01
  17  ST+ X
  18  X^2
  19  ENTER^
  20  ENTER^
  21  STO 00
  22  6949
  23  *
  24  296492
  25  +
  26  *
  27  1248002
  28  -
  29  *
  30  7414380
  31  +
  32  *
  33  5853627
  34  -
  35  64
  36  *
  37  105
  38  /
  39  CHS
  40  STO 07     
  41  STO 08
  42  RCL 01 
  43  RCL 02
  44  8
  45  ST* Z
  46  *
  47  STO 04
  48  X^2
  49  STO 06
  50  ST/ T
  51  RDN
  52  STO 03
  53  X^2
  54  STO 05
  55  /
  56  RCL 00
  57  83
  58  *
  59  2075
  60  +
  61  RCL 00 
  62  *
  63  3039
  64  -
  65  RCL 00     
  66  *
  67  3537
  68  +
  69  32
  70  *
  71  15
  72  /
  73  ST- Z
  74  -
  75  RCL 05
  76  /
  77  X<>Y
  78  RCL 06
  79  /
  80  RCL 00
  81  7
  82  *
  83  82
  84  +
  85  RCL 00     
  86  *
  87  9
  88  -
  89  .75
  90  /
  91  ST- Z
  92  -
  93  RCL 06
  94  /
  95  X<>Y
  96  RCL 05
  97  /
  98  RCL 00
  99  3
100  +
101  ST- Z
102  -
103  RCL 03
104  /
105  X<>Y
106  RCL 04
107  /
108  RCL 03     
109  7
110  Y^X
111  ST/ 07
112  CLX
113  RCL 04
114  LASTX
115  Y^X
116  ST/ 08
117  CLX
118  RCL 02
119  +
120  X<>Y
121  RCL 01
122  +
123  END

 
     ( 187 bytes / SIZE 009 )
 
 

      STACK        INPUTS      OUTPUTS
           Y          n >= 0          y'n,s
           X              s           j'n,s

 
Examples:

  PI   ENTER^
 12    XEQ "dJYZ"   >>>>   j'Pi,12 = 40.14532124            R07 = -57 E-9
                                RDN   y'Pi,12 = 41.72112783            R08 = -43 E-9

  0    ENTER^
  4       R/S            >>>>    j'0,4  = 10.17346816            R07 =  147 E-9
                              RDN   y'0,4  = 11.74915483           R08 =   54 E-9

  0    ENTER^
  3       R/S            >>>>    j'0,3  = 7.015587            R07 =  1930 E-9
                              RDN   y'0,3  = 8.596006            R08 =   474 E-9

Notes:

-Registers R07 & R08 contain the 5th terms of the expansions: this may be used to give an idea of the accuracy.
-In the first 2 examples, almost all the decimals are correct. In the 3rd example, the results are exact to 6D
-Otherwise, delete lines 108 to 117 and lines  41-40.
 

12°)  Bessel Functions - Complex order & arguments
 

    a) Bessel Functions of the 1st kind
 

Data Registers:  R00 thru R10: temp

Flag:  F01     CF 01  to compute  Jn(z)
                        SF 01  to compute   In(z)

Subroutines:   "GAMZ+"                          ( cf "Gamma Function for the HP-41" § 1°)-h) )
                         "Z^2"  "Z*Z"  "Z^Z"  "Z/Z"  ( cf "Complex Functions for the HP-41" )

     "Z^2"  "Z*Z" ... and so on ...  may replaced by M-Code routines ( cf "A few M-Code Routines for the HP-41" )
 
 

01  LBL "JNZ"
02  STO 01           
03  RDN
04  STO 02
05  RDN
06  STO 05
07  X<>Y
08  STO 06 
09  2
10  ST/ 01
11  ST/ 02
12  RCL 02
13  RCL 01
14  XEQ "Z^2"
15  STO 09
16  X<>Y
17  STO 10
18  CLX
19  STO 04           
20  STO 08
21  SIGN
22  STO 00
23  STO 03
24  STO 07 
25  LBL 01
26  RCL 10
27  RCL 09
28  RCL 04
29  RCL 03
30  XEQ "Z*Z"
31  RCL 00
32  FC? 01
33  CHS
34  ST/ Z
35  /
36  RCL 05           
37  RCL 00
38  +
39  RCL 06 
40  X<>Y
41  XEQ "Z/Z"
42  STO 03
43  RCL 07
44  +
45  STO 07
46  LATSX
47  -
48  X<>Y
49  STO 04
50  RCL 08
51  +
52  STO 08
53  LASTX
54  -
55  R-P
56  ISG 00
57  CLX
58  X#0?
59  GTO 01
60  RCL 02          
61  RCL 01 
62  RCL 06
63  RCL 05
64  XEQ "Z^Z"
65  RCL 08
66  RCL 07
67  XEQ "Z*Z"
68  STO 07
69  X<>Y
70  STO 08
71  RCL 06
72  RCL 05
73  1
74  +
75  XEQ "GAMZ+"
76  RCL 08           
77  RCL 07
78  R^
79  R^
80  XEQ "Z/Z"
81  END

 
    ( 125 bytes / SIZE 011 )
 
 

      STACK        INPUTS     OUTPUTS
           T             b            /
           Z             a            /
           Y             y      Im [ f(z) ]
           X             x      Re [ f(z) ]

     n = a + i b               f(z) = Jn(z)  if flag F01 is clear
     z = x + i y                f(z) = In(z)  if flag F01 is set.

Examples:       n = 1 + 4 i  ,   z = 2 + 3 i

   •   CF 01

   4   ENTER^
   1   ENTER^
   3   ENTER^
   2   XEQ "JNZ"   >>>>   J1+4i ( 2+3.i ) = 0.342267397 - 0.424349408 i        ---Execution time = 33s---

   •   SF 01

   4   ENTER^
   1   ENTER^
   3   ENTER^
   2   XEQ "JNZ"   >>>>   I1+4i ( 2+3.i ) = 1.391608538 + 0.315041667 i        ---Execution time = 33s---
 

     b) Bessel Functions of the 2nd kind
 

Data Registers:  R00 thru R14: temp

Flag:  F01     CF 01  to compute  Yn(z)
                        SF 01  to compute   Kn(z)

Subroutines:   "JNZ"   ( cf §12-a) )
                         "Z*Z"  "Z/Z"  "COSZ"  "SINZ"  ( cf "Complex Functions for the HP-41" )
 
 

01  LBL "YNZ"
02  STO 11         
03  X<>Y
04  STO 12 
05  X<>Y
06  XEQ "JNZ"
07  FS? 01
08  GTO 00
09  STO 13
10  X<>Y
11  STO 14
12  RCL 06
13  RCL 05 
14  PI
15  ST* Z
16  *
17  XEQ "COSZ"
18  RCL 14         
19  RCL 13
20  XEQ "Z*Z"
21  LBL 00
22  STO 13
23  X<>Y
24  STO 14
25  RCL 06 
26  CHS
27  RCL 05         
28  CHS
29  RCL 12
30  RCL 11
31  XEQ "JNZ"
32  ST- 13
33  X<>Y
34  ST- 14
35  RCL 06
36  RCL 05         
37  PI
38  ST* Z
39  *
40  XEQ "SINZ"
41  RCL 14
42  RCL 13
43  R^
44  R^
45  XEQ "Z/Z"
46  PI
47  2
48  /
49  FC? 01
50  SIGN
51  FC? 01          
52  CHS
53  ST* Z
54  *
55  END

 
   ( 98 bytes / SIZE 015 )
 
 

      STACK        INPUTS     OUTPUTS
           T             b            /
           Z             a            /
           Y             y      Im [ f(z) ]
           X             x      Re [ f(z) ]

     n = a + i b               f(z) = Jn(z)  if flag F01 is clear
     z = x + i y                f(z) = In(z)  if flag F01 is set.

Examples:       n = 1 + 4 i  ,   z = 2 + 3 i

   •   CF 01

   4   ENTER^
   1   ENTER^
   3   ENTER^
   2   XEQ "YNZ"   >>>>   Y1+4i ( 2+3.i ) = -0.552969731 - 0.212458380 i        ---Execution time = 73s---

   •   SF 01

   4   ENTER^
   1   ENTER^
   3   ENTER^
   2   XEQ "YNZ"   >>>>   K1+4i ( 2+3.i ) = 0.011395936 - 0.065079638 i        ---Execution time = 71s---

Note:

-This program does not work if  n is a real integer.
 

     c) Bessel Functions of the 2nd kind ( real integer order )
 
 

Data Registers:  R00 thru R12: temp   R01-R02 = the result ,  R05 = n

Flag:  F01     CF 01  to compute  Yn(z)
                        SF 01  to compute   Kn(z)

Subroutines:   "JNZ"   ( cf §12-a) )
                         "Z*Z"  "Z^2"  "Z^X"  "LNZ"  ( cf "Complex Functions for the HP-41" )
 
 

  01  LBL "YNZ1"
  02  X<>Y
  03  STO 12        
  04  X<>Y
  05  STO 11
  06  0
  07  RDN
  08  XEQ "JNZ"
  09  STO 01
  10  X<>Y
  11  STO 02
  12  2
  13  ST/ 11
  14  ST/ 12
  15  RCL 12
  16  RCL 11
  17  XEQ "LNZ"
  18  1
  19  FS? 01
  20  CHS
  21  STO 00
  22  RCL 05
  23  STO 10
  24  Y^X
  25  ST+ X
  26  ST* Z
  27  *
  28  RCL 02
  29  RCL 01
  30  XEQ "Z*Z"
  31  STO 01
  32  X<>Y
  33  STO 02
  34  RCL 12
  35  RCL 11
  36  XEQ "Z^2"
  37  STO 03
  38  X<>Y
  39  STO 04        
  40  CLX
  41  STO 07
  42  STO 08
  43  LBL 01
  44  RCL 05
  45  RCL 10
  46  -
  47  FACT
  48  RCL 10
  49  1
  50  -
  51  X<0?
  52  CLST
  53  STO 09
  54  FACT
  55  /
  56  RCL 00
  57  RCL 09
  58  Y^X
  59  *
  60  RCL 04
  61  RCL 03
  62  RCL 09
  63  XEQ "Z^X"
  64  R^
  65  ST* Z
  66  *
  67  ST+ 07
  68  X<>Y
  69  ST+ 08
  70  DSE 10
  71  GTO 01
  72  RCL 12
  73  RCL 11
  74  RCL 05        
  75  CHS
  76  XEQ "Z^X"
  77  RCL 08
  78  RCL 07
  79  XEQ "Z*Z"
  80  ST- 01
  81  X<>Y
  82  ST- 02
  83  1.15443133
  84  STO 09
  85  RCL 05
  86  LBL 02
  87  X#0?
  88  1/X
  89  ST- 09
  90  X<> L
  91  DSE X
  92  GTO 02
  93  RCL 00
  94  ST* 11
  95  ST* 12
  96  CHS
  97  ST* 03
  98  ST* 04
  99  CLX
100  STO 08
101  STO 06
102  SIGN
103  STO 10
104  RCL 09
105  RCL 05
106  FACT
107  1/X
108  STO 07
109  *
110  STO 00        
111  LBL 03
112  RCL 04
113  RCL 03
114  RCL 08
115  RCL 07
116  XEQ "Z*Z"
117  STO 07
118  X<>Y
119  STO 08
120  RCL 10
121  RCL 10
122  RCL 05
123  +
124  *
125  ST/ 07
126  ST/ 08
127  RCL 09
128  LASTX
129  1/X
130  -
131  RCL 10
132  1/X
133  -
134  STO 09
135  ISG 10
136  CLX
137  RCL 07
138  *
139  RCL 00
140  +
141  STO 00
142  LASTX
143  -
144  ABS
145  RCL 08        
146  RCL 09
147  *
148  RCL 06
149  +
150  STO 06
151  LASTX
152  -
153  ABS
154  +
155  X#0?
156  GTO 03
157  RCL 12
158  RCL 11
159  RCL 05
160  XEQ "Z^X"
161  RCL 06
162  RCL 00
163  XEQ "Z*Z"
164  ST+ 01
165  X<>Y
166  ST+ 02
167  PI
168  FS? 01
169  -2
170  ST/ 01
171  ST/ 02
172  RCL 02
173  RCL 01
174  END

 
    ( 262 bytes / SIZE 013 )
 
 

      STACK        INPUTS     OUTPUTS
           Z             n            /
           Y             y      Im [ f(z) ]
           X             x      Re [ f(z) ]

  where        z = x + i y  ,  n is a non-negative integer

    and      f(z) = Yn(z)  if flag F01 is clear     ;      f(z) = Kn(z)  if flag F01 is set.

Examples:      n = 3  ,   z = 1 + 2 i

   •   CF 01

   3   ENTER^
   2   ENTER^
   1   XEQ "YNZ1"   >>>>   Y3 ( 1+2.i ) = 0.290153294 - 0.212118771 i        ---Execution time = 54s---

   •   SF 01

   3   ENTER^
   2   ENTER^
   1   XEQ "YNZ1"   >>>>   K3 ( 1+2.i ) = -0.681436426 + 0.625154655 i        ---Execution time = 54s---

-Likewise,   Y0 ( 1+2.i ) = 1.367418719 + 1.521506580 i    &   K0 ( 1+2.i ) = -0.242345103 - 0.176267190 i

Note:

-If n is a negative integer, one can use:   Yn(x) = (-1)n Y -n(x)  and  Kn(x) = K -n(x)
 

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4