Coulomb

Coulomb Wave Functions for the HP-41

Overview

1°) Regular Coulomb Wave Function

   a) Program#1
   b) Program#2
   c) Non-integer values of L

2°) Irregular Coulomb Wave Function

   a) Program#1
   b) Program#2
   c) Non-integer values of L

3°) Asymptotic Expansion

-The regular Coulomb wave function F_L(n,r) and the irregular Coulomb wave function G_L(n,r) are 2 independant solutions of the differential equation:

d²y/dr² + [ 1 - 2n/r - L(L+1)/r² ] y = 0

1°) Regular Coulomb Wave Function

a) Program#1

Formulae: F_L(n,r) = C_L(n) r^L+1 Sum_k>L A_k^L (n) r^k-L-1

with C_L(n) = (1/Gam(2L+2)) 2^L e^-pi.n/2 | Gam(L+1+i.n) |
and A_L+1^L = 1 ; A_L+2^L = n/(L+1) ; (k+L)(k-L-1) A_k^L = 2n A_k-1^L - A_k-2^L ( k > L+2 )

Data Registers: R00 thru R05: temp

R06 = r ; R07 = n ; R08 = L ; R09 = C_L(n)
Flags: /
Subroutine: "GAMZ" ( Gamma Function in the complex plane )

01 LBL "RCWF"
02 STO 06
03 RDN
04 STO 07
05 X<>Y
06 STO 08
07 1
08 +
09 XEQ "GAMZ"
10 LASTX
11 RCL 07
12 PI
13 *
14 2
15 /

16 E^X
17 /
18 2
19 RCL 08
20 Y^X
21 *
22 RCL 08
23 ST+ X
24 1
25 +
26 FACT
27 /
28 STO 09
29 CLX
30 STO 01

31 SIGN
32 STO 02
33 STO 03
34 LBL 01
35 2
36 STO 04
37 X<>Y
38 LBL 02
39 RCL 07
40 ST+ X
41 RCL 02
42 ST* Y
43 X<> 01
44 RCL 06
45 *

46 -
47 RCL 06
48 *
49 RCL 03
50 ST/ Y
51 RCL 08
52 ST+ X
53 +
54 1
55 ST+ 03
56 +
57 /
58 STO 02
59 +
60 X#Y?

61 GTO 01
62 DSE 04
63 GTO 02
64 RCL 09
65 *
66 RCL 06
67 ST* Y
68 RCL 08
69 Y^X
70 *
71 END

( 96 bytes / SIZE 010 )

STACK	INPUTS	OUTPUTS
Z	L	/
Y	n	/
X	r	F_L(n,r)

      L is a non-negative integer,
      n is real,
      r is non-negative.

Example:

     2   ENTER^
   0.7 ENTER^
   1.8 XEQ "RCWF" >>>>   F₂( 0.7 , 1.8 ) = 0.141767744   ( in 34 seconds )

b) Program#2

-This program uses the same basic formulae but | Gam( L+1+i n ) | is calculated without "GAMZ"

-Instead, we have also | Gam( 1+i y ) |² = (Pi.y) / Sinh (Pi y) and

| Gam( 1+L+i y ) |² = [ L² + y² ] [ (L-1)² + y² ] .................. [ 1 + y² ] (Pi.y) / Sinh (Pi y)

Data Registers: R01 = r ; R02 = n ; R03 = L ; R00 & R04 thru R06: temp
Flag: F10
Subroutines: /

-Line 30 is a synthetic TEXT0 or another NOP like LBL 00

01 LBL "RCWF2"
02 STO 01
03 RDN
04 STO 02
05 X<>Y
06 STO 00
07 STO 03
08 X=0?
09 GTO 00
10 SIGN
11 LBL 01
12 RCL 00
13 X^2
14 RCL 02
15 X^2
16 +
17 *
18 DSE 00
19 GTO 01
20 LBL 00

21 X=0?
22 SIGN
23 RCL 02
24 PI
25 *
26 STO 00
27 ENTER^
28 X=0?
29 ISG Y
30 ""
31 E^X-1
32 LASTX
33 CHS
34 E^X-1
35 -
36 2
37 /
38 X=0?
39 SIGN
40 /

41 *
42 SQRT
43 2
44 ST/ 00
45 RCL 03
46 Y^X
47 *
48 RCL 00
49 E^X
50 /
51 RCL 03
52 ST+ X
53 1
54 +
55 FACT
56 /
57 STO 00
58 CLX
59 STO 04
60 SIGN

61 STO 05
62 STO 06
63 LBL 02
64 SF 10
65 LBL 03
66 RCL 02
67 ST+ X
68 RCL 05
69 ST* Y
70 X<> 04
71 RCL 01
72 *
73 -
74 RCL 01
75 *
76 RCL 06
77 ST/ Y
78 RCL 03
79 ST+ X
80 +

81 1
82 ST+ 06
83 +
84 /
85 STO 05
86 +
87 X#Y?
88 GTO 02
89 FS?C 10
90 GTO 03
91 RCL 00
92 *
93 RCL 01
94 ST* Y
95 RCL 03
96 Y^X
97 *
98 END

( 125 bytes / SIZE 007 )

STACK	INPUTS	OUTPUTS
Z	L	/
Y	n	/
X	r	F_L(n,r)

       L is a non-negative integer,
       n is real,
       r is non-negative.

Example:

     2   ENTER^
   0.7 ENTER^
   1.8 XEQ "RCWF2" >>>>   F₂( 0.7 , 1.8 ) = 0.141767746                ---Execution time = 17s---

-"RCWF2" is both faster and more accurate than "RCWF"

c) Non-integer Values of L

Formula:

F_L(n,r) = (1/2) [ 1/(2L+1)! ] | Gam(L+1+i.n) | exp [ -(PI) n/2 - i (L+1)(PI)/2 ] M_i.n,L+1/2 (2i.r)

where M_k,µ = Whittaker's function of the 1st kind

-Since it involves complex numbers, we cannot use the program "WHIM" listed in "Whittaker's Function for the HP-41"

Data Registers: R00 thru R11: temp

R06 = r ; R07 = n ; R08 = L
Flags: /
Subroutines: "GAM3" or another "gamma-function program" ( cf "Gamma Function for the HP-41" )
"GAMZ+" or "GAMZ" ( gamma function in the complex plane )

-The M-Code routine Z*Z may be replaced by XEQ "Z*Z" ( cf "Complex Functions for the HP-41" )
-Lines 12-13-14 & 17-18 are only useful to compute the irregular Coulomb wave function with "ICWF3"
-Otherwise, they may be deleted.

01 LBL "RCWF3"
02 STO 06
03 RDN
04 STO 07
05 X<>Y
06 STO 08
07 1
08 +
09 STO 09
10 XEQ "GAMZ+"
11 DEG
12 STO 10
13 X<>Y
14 STO 11
15 R-P
16 STO 05
17 ST/ 10
18 ST/ 11
19 CLX
20 STO 00

21 STO 02
22 STO 04
23 SIGN
24 STO 01
25 STO 03
26 LBL 01
27 ISG 00
28 CLX
29 RCL 00
30 RCL 08
31 +
32 RCL 06
33 ST+ X
34 ST* Y
35 RCL 07
36 *
37 RCL 08
38 RCL 09
39 +
40 RCL 00

41 ST+ Y
42 *
43 ST/ Z
44 /
45 RCL 04
46 RCL 03
47 Z*Z
48 STO 03
49 RCL 01
50 +
51 STO 01
52 LASTX
53 -
54 ABS
55 X<>Y
56 STO 04
57 RCL 02
58 +
59 STO 02
60 LASTX

61 -
62 ABS
63 +
64 X#0?
65 GTO 01
66 RCL 06
67 R-D
68 CHS
69 LASTX
70 ST+ X
71 RCL 09
72 Y^X
73 PI
74 RCL 07
75 *
76 2
77 /
78 E^X
79 /
80 P-R

81 RCL 02
82 RCL 01
83 Z*Z
84 STO 03
85 X<>Y
86 STO 04
87 RCL 09
88 ST+ X
89 XEQ "GAM3"
90 RCL 05
91 X<>Y
92 ST+ X
93 /
94 ST* 03
95 ST* 04
96 RCL 04
97 RCL 03
98 END

( 134 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
Z	L	/
Y	n	Im(F_L) ~ 0
X	r > 0	F_L(n,r)

Example:

   1.4   ENTER^
-2.1   ENTER^
   1.6   XEQ "RCWF3" >>>>   F_1.4( -2.1 , 1.6 ) = -0.779739876                ---Execution time = 43s---

-If there were no roundoff errors, Y-output = 0
-Actually, the small number in Y-register gives an idea of the accuracy of the results:
-If, for instance, Y = E-4 , there are probably no more than 4 exact decimals.
-Here, Y = -2 E-9 , so the result is probably correct to 8D ( an HP-48 returns -0.779739873 )

-This program also works if L is an integer.

2°) Irregular ( Logarithmic ) Coulomb Wave Function

a) Program#1

Formulae: G_L(n;r) = (1/Pi).(exp(2.Pi.n) - 1) F_L(n;r) [ Ln(2r) + q_L(n)/p_L(n) ] + (1/C_L(n)/r^L/(2L+1)) Sum_k>-L-1 a_k^L (n) r^k+L

with a_-L-1^L = 0 ; a_-L^L = 1 ; (k+L)(k-L-1) a_k^L = 2n a_k-1^L - a_k-2^L - (2k-1).p_L(n).A_k^L and a_L+1^L = 0

p_L(n) = 2n ( 1+n² )( 4+n² ) ........... ( L²+n² ) 2^2L / ( 2L+1) / [(2L)!]²

q_L(n)/p_L(n) = Sum_s=1,...,L s/(s²+n²) - 1/1 - 1/2 - 1/3 - ...... - 1/(2L+1) + Ré ( Psi(1+i.n) ) + 2.gamma + r_L(n)/p_L(n) ( gamma = Euler's Constant )

where r_L(n) = ( (-1)^L+1 / (2L)! ) Im [ 1/(2L+1) + 2(i.n-L)/1!/(2L) + 2²(i.n-L)(i.n-L+1)/2!/(2L-1)
+ ...................................... + 2^2L(i.n-L)(i.n-L+1) ....... (i.n+L-1)/(2L)! ]

Data Registers: R00 thru R05 & R11 thru R15: temp

R06 = r ; R07 = n ; R08 = L ; R09 = C_L(n) ; R10 = F_L(n;r)
Flags: /
Subroutines: "RCWF" ( regular Coulomb wave function )
"PSIZ" ( Digamma function in the complex plane )

-Line 106 is 2 x Euler's Constant
-Line 178 is TEXT 0 or another NOP instruction like STO X

01 LBL "ICWF"
02 XEQ "RCWF"
03 STO 10
04 RCL 07
05 ST+ X
06 STO 15
07 PI
08 *
09 E^X-1
10 *
11 PI
12 /
13 STO 11
14 RCL 07
15 1
16 XEQ "PSIZ"
17 STO 12
18 RCL 08
19 ST+ X
20 STO 03
21 1
22 +
23 LBL 01
24 1/X
25 ST- 12
26 LASTX
27 DSE X
28 GTO 01
29 STO 02
30 STO 04
31 RCL 15
32 4
33 RCL 08
34 Y^X
35 *
36 1
37 STO 13
38 RCL 03
39 ST+ Y
40 FACT

41 STO 05
42 X^2
43 *
44 /
45 STO 00
46 RCL 08
47 STO 01
48 X=0?
49 GTO 00
50 LBL 02
51 RCL 01
52 ENTER^
53 X^2
54 RCL 07
55 X^2
56 +
57 ST* 00
58 /
59 ST+ 12
60 DSE 01
61 GTO 02
62 LBL 03
63 RCL 07
64 RCL 08
65 RCL 03
66 -
67 R-P
68 X<>Y
69 RCL 02
70 +
71 STO 02
72 X<>Y
73 RCL 13
74 *
75 ST+ X
76 RCL 08
77 ST+ X
78 RCL 03
79 -
80 1

81 +
82 /
83 STO 13
84 P-R
85 CLX
86 RCL 03
87 /
88 ST+ 04
89 DSE 03
90 GTO 03
91 LBL 00
92 RCL 12
93 1
94 CHS
95 RCL 08
96 Y^X
97 RCL 04
98 *
99 RCL 05
100 /
101 RCL 00
102 X=0?
103 SIGN
104 /
105 -
106 1.15443133
107 +
108 RCL 06
109 ST+ X
110 LN
111 +
112 ST* 11
113 CLX
114 STO 02
115 STO 03
116 SIGN
117 STO 04
118 STO 05
119 STO 13
120 LBL 04

121 2
122 STO 14
123 LBL 05
124 RCL 15
125 RCL 04
126 ST* Y
127 X<> 03
128 RCL 06
129 *
130 -
131 RCL 06
132 *
133 RCL 15
134 RCL 02
135 ST* Y
136 X<> 01
137 RCL 06
138 *
139 -
140 RCL 06
141 *
142 RCL 05
143 ST/ Y
144 RCL 08
145 ST+ X
146 -
147 1
148 -
149 STO 12
150 X#0?
151 GTO 00
152 RDN
153 CLX
154 RCL 06
155 RCL 05
156 Y^X
157 1
158 LBL 00
159 /
160 STO 02

161 RCL 00
162 *
163 RCL 05
164 RCL 08
165 -
166 ST+ X
167 1
168 -
169 *
170 -
171 RCL 05
172 /
173 RCL 12
174 X#0?
175 /
176 STO 04
177 ISG 05
178 ""
179 RCL 13
180 +
181 STO 13
182 LASTX
183 X#Y?
184 GTO 04
185 DSE 14
186 GTO 05
187 RCL 06
188 RCL 08
189 Y^X
190 RCL 09
191 *
192 RCL08
193 ST+ X
194 1
195 +
196 *
197 /
198 RCL 11
199 +
200 END

( 259 bytes / SIZE 016 )

STACK	INPUTS	OUTPUTS
Z	L	/
Y	n	/
X	r	G_L(n,r)

L is a non-negative integer, n is real and r is positive. F_L(n,r) is stored in R10

Example:

            3   ENTER^
         -0.4 ENTER^
          1.2 XEQ "ICWF" >>>> G₃( -0.4 , 1.2 ) = 6.563265438        ( in 102 seconds )
                               and we have F₃( -0.4 , 1.2 ) = 0.029547268          in register R10.

-Unfortunately, several significant digits are sometimes cancelled when the last operation ( line 199 ) is performed.
( Check the value in LASTX or in R11 )

b) Program#2

-Here, Ré [ Psi( 1+i.x) ] is computed without calling "PSIZ"

-We use: Ré [ Psi( 1+i.x) ] = -gamma + x² SUM_{k=1,2,.........} 1/[k(k²+x²)] ( gamma = Euler's Constant )

-This infinite sum is actually calculated with the help of Euler-Mac Laurin formula, and after some algebra it yields:

Ré [ Psi( 1+i.x) ] = - SUM_{k=1,2,.....15} k/(k²+x²) + (1/2) Ln ( 16²+x² ) - (16/2)( 16²+x² )^-1 + (1/12)( 16²+x² )^-2 [ 0.1 + x² - 16² - (24/30) 16² x² ( 16²+x² )^-2 ]
+ eps

where eps is a very small number of the order of -2.9 10^-10

Data Registers:   R01 = r , R02 = n , R03 = L , R04 = F_L(n,r)    R00 & R05 thru R15: temp
Flag:  F10
Subroutine: "RCWF2"    ( § 1°) b) )

-You could perhaps add 3 E-10 - after line 58 though it dosn't seem to change the results significantly.

01 LBL "ICWF2"
02 XEQ "RCWF2"
03 STO 04
04 RCL 02
05 ST+ X
06 STO 05
07 PI
08 *
09 E^X-1
10 *
11 PI
12 /
13 STO 06
14 RCL 02
15 X^2
16 STO 15
17 15
18 STO 14
19 CLX
20 LBL 01
21 RCL 14
22 ENTER^
23 X^2
24 RCL 15
25 +
26 /
27 -
28 DSE 14
29 GTO 01
30 RCL 15
31 256
32 STO 14
33 +
34 STO 08
35 SQRT
36 LN
37 +
38 8
39 RCL 08
40 ST* 08
41 /

42 -
43 204.8
44 RCL 15
45 *
46 RCL 08
47 /
48 RCL 14
49 .1
50 -
51 RCL 15
52 -
53 +
54 12
55 /
56 RCL 08
57 /
58 -
59 STO 07
60 RCL 03
61 ST+ X
62 STO 08
63 1
64 +
65 LBL 02
66 1/X
67 ST- 07
68 LASTX
69 DSE X
70 GTO 02
71 STO 09
72 STO 10
73 RCL 05
74 4
75 RCL 03
76 Y^X
77 *
78 1
79 STO 11
80 RCL 08
81 ST+ Y
82 FACT

83 STO 12
84 X^2
85 *
86 /
87 STO 14
88 RCL 03
89 STO 13
90 X=0?
91 GTO 00
92 LBL 03
93 RCL 13
94 ENTER^
95 X^2
96 RCL 15
97 +
98 ST* 14
99 /
100 ST+ 07
101 DSE 13
102 GTO 03
103 LBL 04
104 RCL 03
105 RCL 08
106 -
107 STO 15
108 RCL 11
109 *
110 RCL 02
111 RCL 09
112 *
113 -
114 X<> 11
115 RCL 02
116 *
117 RCL 15
118 RCL 09
119 *
120 +
121 2
122 RCL 03
123 RCL 15

124 +
125 1
126 +
127 /
128 ST* 11
129 *
130 STO 09
131 RCL 08
132 /
133 ST+ 10
134 DSE 08
135 GTO 04
136 LBL 00
137 RCL 07
138 1
139 CHS
140 RCL 03
141 Y^X
142 RCL 10
143 *
144 RCL 12
145 /
146 RCL 14
147 X=0?
148 SIGN
149 /
150 -
151 1.15443133
152 +
153 RCL 01
154 ST+ X
155 LN
156 +
157 ST* 06
158 CLX
159 STO 08
160 STO 09
161 SIGN
162 STO 10
163 STO 11
164 STO 12

165 LBL 05
166 SF 10
167 LBL 06
168 RCL 05
169 RCL 10
170 ST* Y
171 X<> 08
172 RCL 01
173 *
174 -
175 RCL 01
176 *
177 RCL 05
178 RCL 09
179 ST* Y
180 X<> 13
181 RCL 01
182 *
183 -
184 RCL 01
185 *
186 RCL 12
187 ST/ Y
188 RCL 03
189 ST+ X
190 -
191 1
192 -
193 STO 07
194 X#0?
195 GTO 00
196 RDN
197 CLX
198 RCL 01
199 RCL 12
200 Y^X
201 1
202 LBL 00
203 /
204 STO 09
205 RCL 14

206 *
207 RCL 12
208 RCL 03
209 -
210 ST+ X
211 1
212 -
213 *
214 -
215 RCL 12
216 /
217 RCL 07
218 X#0?
219 /
220 STO 10
221 ISG 12
222 ""
223 RCL 11
224 +
225 STO 11
226 LASTX
227 X#Y?
228 GTO 05
229 FS?C 10
230 GTO 06
231 RCL 01
232 RCL 03
233 Y^X
234 RCL 00
235 *
236 RCL 03
237 ST+ X
238 1
239 +
240 *
241 /
242 RCL 06
243 +
244 END

( 313 bytes / SIZE 016 )

STACK	INPUTS	OUTPUTS
Z	L	/
Y	n	/
X	r	G_L(n,r)

    L is a non-negative integer,
    n is real,
    r is positive    &    F_L(n,r) is stored in R04

Example:

            3   ENTER^
         -0.4 ENTER^
          1.2 XEQ "ICWF" >>>> G₃( -0.4 , 1.2 ) = 6.563265348                                         ---Execution time = 66s---
                               and we have F₃( -0.4 , 1.2 ) = 0.029547268          in register R04.

-Like with "ICWF", several significant digits may be cancelled when the last operation is performed.
( Check the value in LASTX or in R06 )

c) Non-integer Values of L

Formula:

G_L(n,r) = i F_L(n,r) + | Gam(L+1+i.n) | [ 1/Gam(L+1+i.n) ] exp [ (PI) n/2 + i L(PI)/2 ] W_i.n,L+1/2 (2i.r)

where W_k,µ = Whittaker's function of the 2nd kind

Data Registers:   R00 thru R15: temp
                              R06 = r , R07 = n , R08 = L
Flags: /
Subroutines: "GAM3" or another "gamma-function program" ( cf "Gamma Function for the HP-41" )
                         "GAMZ+" or "GAMZ"    ( gamma function in the complex plane )
                         "RCWF3" cf § 1-c)

-The M-Code routines Z*Z & Z/Z may be replaced by XEQ "Z*Z" & XEQ "Z/Z" ( cf "Complex Functions for the HP-41" )

01 LBL "ICWF3"
02 XEQ "RCWF3"
03 RCL 08
04 RCL 09
05 +
06 90
07 *
08 RCL 06
09 R-D
10 -
11 RCL 06
12 ST+ X
13 RCL 09
14 Y^X
15 PI
16 RCL 07
17 *
18 2
19 /
20 E^X
21 *
22 P-R
23 RCL 02
24 RCL 01
25 Z*Z
26 STO 12
27 X<>Y
28 STO 13
29 RCL 08

30 RCL 09
31 +
32 CHS
33 XEQ "GAM3"
34 ST* 12
35 ST* 13
36 CLX
37 STO 00
38 STO 02
39 STO 04
40 SIGN
41 STO 01
42 STO 03
43 LBL 01
44 RCL 00
45 RCL 08
46 -
47 RCL 06
48 ST+ X
49 ST* Y
50 RCL 07
51 *
52 RCL 00
53 RCL 08
54 ST+ X
55 -
56 ISG 00
57 CLX
58 RCL 00

59 *
60 ST/ Z
61 /
62 RCL 04
63 RCL 03
64 Z*Z
65 STO 03
66 RCL 01
67 +
68 STO 01
69 LASTX
70 -
71 ABS
72 X<>Y
73 STO 04
74 RCL 02
75 +
76 STO 02
77 LASTX
78 -
79 ABS
80 +
81 X#0?
82 GTO 01
83 RCL 06
84 R-D
85 CHS
86 LASTX
87 ST+ X

88 RCL 08
89 CHS
90 Y^X
91 PI
92 RCL 07
93 *
94 2
95 /
96 E^X
97 *
98 P-R
99 RCL 02
100 RCL 01
101 Z*Z
102 STO 14
103 X<>Y
104 STO 15
105 RCL 08
106 RCL 09
107 +
108 XEQ "GAM3"
109 ST* 14
110 ST* 15
111 RCL 07
112 CHS
113 RCL 09
114 XEQ "GAMZ+"
115 RCL 15
116 RCL 14

117 R^
118 R^
119 Z/Z
120 STO 14
121 X<>Y
122 STO 15
123 RCL 07
124 CHS
125 RCL 08
126 CHS
127 XEQ "GAMZ+"
128 RCL 13
129 RCL 12
130 R^
131 R^
132 Z/Z
133 X<>Y
134 RCL 15
135 +
136 X<>Y
137 RCL 14
138 +
139 RCL 11
140 RCL 10
141 Z/Z
142 END

( 191 bytes / SIZE 016 )

STACK	INPUTS	OUTPUTS
Z	L	/
Y	n	/
X	r > 0	G_L(n,r)

Example:

   1.4   ENTER^
-2.1   ENTER^
   1.6   XEQ "ICWF3" >>>>   G_1.4( -2.1 , 1.6 ) = -0.206368268                         ---Execution time = 107s---

-The first term: i F_L(n,r) is not taken into account since this complex number is purely imaginary.
-We can, however, add this term to the final result to get an idea of the accuracy:
-Add X<>Y RCL 16 + X<>Y after line 141 and add STO 16 after line 02
-Here, it yields Y-output = 2 E-9 suggesting 8 or 9 correct decimals ( it should be 0 if there were no roundoff errors ).
-The HP-48 returns -0.206368268 so the result is exact to 9D.

-This program does not work if L is an integer because we cannot use the same formulae to compute Whittaker's function of the 2nd kind.

3°) Asymptotic Expansion

-This routine is useful to compute the Coulomb Wave Functions for large value of r.
-It also works if L is not an integer.

Formulae:

F_L(n,r) = g cos µ + f sin µ where µ = r - n Ln (2r) - L.(Pi)/2 + Arg Gam ( 1 + L + i n )
G_L(n,r) = f cos µ - g sin µ and f + i.g ~ u₀ + u₁ + u₂ + ............ + u_k + ......................

with u₀ = 1 and u_k / u_k-1 = ( i n - L + k - 1 ) ( i n + L + k ) / ( 2k i r ) ( k > 0 )

Data Registers: R00 thru R08: temp
Flags: /
Subroutines: "GAMZ+" ( Gamma Function in the Complex Plane - cf "Gamma Function for the HP-41" )

-The M-Code routine Z*Z may be replaced by XEQ "Z*Z" ( cf "Complex Functions for the HP-41" )

01 LBL "AECWF"
02 STO 04
03 RDN
04 STO 05
05 X<>Y
06 STO 06
07 1
08 +
09 XEQ "GAMZ+"
10 DEG
11 R-P
12 CLX
13 RCL 06
14 90
15 *
16 -
17 RCL 04

18 ENTER^
19 ST+ X
20 LN
21 RCL 05
22 *
23 -
24 R-D
25 +
26 STO 00
27 CLX
28 STO 02
29 STO 03
30 STO 08
31 SIGN
32 STO 01
33 STO 07
34 LBL 01

35 RCL 05
36 RCL 03
37 RCL 06
38 -
39 ISG 03
40 CLX
41 RCL 03
42 RCL 06
43 +
44 CHS
45 RCL 05
46 Z*Z
47 RCL 08
48 RCL 07
49 Z*Z
50 RCL 03
51 ST+ X

52 RCL 04
53 *
54 ST/ Z
55 /
56 STO 07
57 RCL 01
58 +
59 STO 01
60 LASTX
61 -
62 ABS
63 X<>Y
64 STO 08
65 RCL 02
66 +
67 STO 02
68 LASTX

69 -
70 ABS
71 +
72 RND
73 X#0?
74 GTO 01
75 RCL 00
76 RCL 02
77 P-R
78 RCL 00
79 RCL 01
80 P-R
81 R^
82 -
83 X<> Z
84 +
85 END

( 110 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
Z	L	/
Y	n	G_L(n,r)
X	r	F_L(n,r)

Example: L = 3 , n = 4 , r = 16

FIX 8

   3   ENTER^
   4   ENTER^
16 XEQ "AECWF" >>>>    F₃(4,16) = -0.997776721                              ---Execution time = 38s---
                                   RDN     G₃(4,16) = -0.694942604

Notes:

-In fact, these series are divergent and line 72 ( RND ) is used to stop the loop when 2 rounded consecutive sums are equal
-Check registers R07 & R08 to see the last u_k = -2 E-11 - 4 E-9 i in the example above.

-If r is even larger, the routine works in FIX 9 or even in SCI 9
-For instance, "AECWF" returns F₀(1,20) = -0.3292255394 and G₀(1,20) = -0.9724283978

-This program can be modified so that the loop stops when 2 successive u_k start to increase:

Replace line 83-84 by RCL 09 X<>Y RDN RDN +
Replace line 72-73-74 by X=0? GTO 00 ENTER^ X<> 09 X>Y? GTO 01 LBL 00
and add STO 09 after line 14

-Register Z and R09 will give an idea of the accuracy.

-It works often ... but not always! A premature stop is possible.
-But you can press XEQ 01 to check if the loop is still going on.

Reference:

[1] Abramowitz and Stegun , "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4

hp41programs

Coulomb Wave Functions for the HP-41