Spheroidal Harmonics for the HP-41
Overview
1°) Eigenvalues
a) Program#1
b) Program#2
2°) Spheroidal Wave Functions ( | x | <= 1 )
a) Program#1
b) Program#2
c) A Generalization
3°) Other Normalizations
a) Normalization factors
b) Normalized Functions ( including
"Meixner-Shaefke" Scheme )
-We want to solve the differential equation ( 1 - x2 ) d2S/dx2 - 2.x dS/dx + [ Lmn - c2 x2 - m2/( 1 - x2 ) S ] = 0
where m , n are non-negative integers: m =
0 , 1 , 2 , 3 , ....... n = m , m + 1 , m + 2 , ......
label the successive eigenvalues Lmn
and c2 may be positive or negative ( i-e
c may be purely imaginary )
c2 > 0 corresponds to the prolate spheroidal
wave functions
c2 < 0 corresponds to the oblate
spheroidal wave functions
( c = 0 leads to the associated Legendre functions )
-But first, we have to find the numbers Lmn(c2)
so that the solutions are regular at x = +/- 1
( if c = 0 , the eigenvalues are simply n ( n + 1 ) )
1°) Eigenvalues
a) Program#1
-The program below finds Lmn by solving the transcendental equation U1(Lmn) + U2(Lmn) = 0 where U1 & U2 are 2 continued fractions:
-brm
-br-2m
U1(Lmn) = grm
- Lmn + ---------------
--------------- ..............
gr-2m - Lmn +
gr-4m - Lmn +
r = n - m
-br+2m
-br+4m
U2(Lmn) =
--------------- ----------------
..............
gr+2m - Lmn +
gr+4m - Lmn +
with brm
= [ r.(r-1).(2m+r).(2m+r-1).c4 ] / [ (2m+2r-1)2.(2m+2r+1).(2m+2r-3)
]
and grm
= (m+r).(m+r+1) + (c2/2).[ 1 - (4m2-1)/((2m+2r-1).(2m+2r+3))
]
Data Registers: R00 to R08 are used by "CFR"
R09 = n-m R11 = m R15 thru R19: unused
R10 = c2 R12 to R14: temp
R20 thru R23 are used by "SOLVE"
Flag: F01
Subroutines: "CFR"
( cf "Continued Fractions for the HP-41" )
"SOLVE" ( cf "Linear and non-linear Systems for the HP-41" ) after
replacing registers R00 thru R03 by R20 thru R23 ( for instance )
-In other words:
LBL "SOLVE"
STO 23
ENTER^ /
-
GTO 01
STO 21
LBL 01
ENTER^ RCL 22
*
END
X<>Y
VIEW 21
X<> 23 RCL 21
+
STO 22
RCL 21
-
STO 22
STO 21 ( 27 lines / 51 bytes
)
XEQ IND 20
XEQ IND 20 X#0?
STO T
X#Y?
| 01 LBL "LMN"
02 STO 10 03 RDN 04 STO 09 05 X<>Y 06 STO 11 07 - 08 X<> 09 09 ENTER^ 10 "T" 11 1 12 + 13 * 14 RCL 10 15 2 16 / 17 + 18 ENTER^ 19 ASTO 20 20 1 21 + 22 "U" 23 ASTO 00 24 CLA |
25 XEQ "SOLVE"
26 RTN 27 LBL "T" 28 STO 01 29 RCL 09 30 STO 12 31 XEQ 01 32 X<> L 33 SF 01 34 XEQ "CFR" 35 STO 14 36 CLX 37 CF 01 38 RCL 01 39 XEQ "CFR" 40 RCL 14 41 + 42 RTN 43 LBL "U" 44 RCL 09 45 RCL 02 46 1 47 FC? 01 48 CLX |
49 -
50 ST+ X 51 FS? 01 52 CHS 53 + 54 STO 12 55 1 56 RCL 12 57 - 58 * 59 RCL 11 60 ST+ X 61 RCL 12 62 + 63 ST* Y 64 1 65 - 66 * 67 RCL 10 68 X^2 69 * 70 RCL 11 71 RCL 12 72 + |
73 ST+ X
74 3 75 - 76 ST/ Y 77 4 78 + 79 ST/ Y 80 2 81 FS? 01 82 ST- 12 83 - 84 X^2 85 / 86 STO 13 87 LBL 01 88 1 89 RCL 11 90 ST+ 12 91 ST+ X 92 X^2 93 - 94 RCL 12 95 ST+ X 96 1 |
97 -
98 ST/ Y 99 4 100 + 101 / 102 1 103 + 104 RCL 10 105 * 106 2 107 / 108 RCL 12 109 RCL 12 110 1 111 + 112 * 113 + 114 RCL 01 115 - 116 RCL 13 117 RTN 118 END |
( 172 bytes / SIZE 024 )
| STACK | INPUTS | OUTPUTS |
| Z | m | Lmn |
| Y | n | Lmn |
| X | c2 | Lmn |
Example1: m = 4 ,
n = 11 , c2 = -1
4 ENTER^
11 ENTER^
-1 XEQ "LMN" >>>> Lmn
= 131.5600809 ( in 1mn27s )
( the HP-41 displays the successive approximations )
Example2: m = 0 , n = 0 , c2 = 3
0 ENTER^ ENTER^
3 XEQ "LMN" >>>> Lmn
= 0.879933945 ( in 2mn30s )
Notes:
-This program uses 2 initial guesses: Li = n(n+1)
+ c2/2 and L'i = 1 + Li
( lines 02 to 21 )
-A better approximation is Li = n(n+1)
+ ( c2/2 ) [ 1 - (2m-1).(2m+1)/((2n-1).(2n+3)) ]
-More accurate L-approximations may be found in the "Handbook of Mathematical
Functions" but the formulas are quite complex.
-Too bad guesses can produce correct eigenvalues, but corresponding
to another value of n.
-For example, if m = n = 0 and c2 = -16 , Li
= 0 and L'i = 1 would yield 0.221407906
which is actually Lmn for m = 0 , n = 2 , c2
= -16
-So, it could be preferable to improve the 2 guesses, especially for
large c-values.
-Here are a few results obtained by "LMN" and Warren Furlow's excellent
V41 emulator ( in "turbo" mode )
( of course - due to roundoff errors - the last decimal may be
doubtful )
***** TABLE OF SPHEROIDAL EIGENVALUES *****
| m n c2 | Lmn | m n c2 | Lmn | m n c2 | Lmn | m n c2 | Lmn |
| 0 0 1 | 0.319000055 | 0 0 -1 | -0.348602400 | 0 1 1 | 2.593084580 | 0 1 -1 | 1.393206310 |
| 4 | 1.127734066 | -4 | -1.594493214 | 4 | 4.287128544 | -4 | -0.505243981 |
| 9 | 2.136732226 | -9 | -4.343292705 | 9 | 6.820888328 | -9 | -3.899400290 |
| 16 | 3.172067425 | -16 | -9.150793382 | 16 | 9.805943842 | -16 | -9.025710674 |
| 25 | 4.195128875 | -25 | -16.07904274 | 25 | 12.91170325 | -25 | -16.05041268 |
| 0 2 1 | 6.533471801 | 0 2 -1 | 5.486800054 | 0 3 1 | 12.51446214 | 0 3 -1 | 11.49212090 |
| 4 | 8.225713002 | -4 | 4.091509101 | 4 | 14.10020387 | -4 | 10.00386381 |
| 9 | 11.19293865 | -9 | 2.251269209 | 9 | 16.88903022 | -9 | 7.611465213 |
| 16 | 15.30630000 | -16 | 0.221407908 | 16 | 21.04896065 | -16 | 4.339082933 |
| 25 | 20.17691472 | -25 | -2.448598906 | 25 | 26.58735961 | -25 | 0.060929887 |
| 1 1 1 | 2.195548355 | 1 1 -1 | 1.795304587 | 1 2 1 | 6.424699144 | 1 2 -1 | 5.567527453 |
| 4 | 2.734111026 | -4 | 1.118553391 | 4 | 7.653149562 | -4 | 4.222747334 |
| 9 | 3.504795867 | -9 | -0.269420805 | 9 | 9.555998398 | -9 | 1.821541436 |
| 16 | 4.399593067 | -16 | -2.909200759 | 16 | 11.94871939 | -16 | -1.869861296 |
| 25 | 5.350422298 | -25 | -7.493388287 | 25 | 14.64295625 | -25 | -7.127837527 |
| 1 3 1 | 12.46791533 | 1 3 -1 | 11.53481845 | 1 4 1 | 20.48169601 | 1 4 -1 | 19.52068334 |
| 4 | 13.88149342 | -4 | 10.16324505 | 4 | 21.94014372 | -4 | 18.09765523 |
| 9 | 16.23846623 | -9 | 8.007074153 | 9 | 24.40831218 | -9 | 15.77725227 |
| 16 | 19.46526054 | -16 | 5.405903145 | 16 | 27.91188168 | -16 | 12.62871852 |
| 25 | 23.39761313 | -25 | 2.750367212 | 25 | 32.42194360 | -25 | 8.694959247 |
| 2 2 1 | 6.140948992 | 2 2 -1 | 5.855162574 | 2 3 1 | 12.33110151 | 2 3 -1 | 11.66440927 |
| 4 | 6.542495275 | -4 | 5.395010784 | 4 | 13.29825047 | -4 | 10.62995129 |
| 9 | 7.151100524 | -9 | 4.526460462 | 9 | 14.82777821 | -9 | 8.809393921 |
| 16 | 7.903860949 | -16 | 3.027624006 | 16 | 16.81295852 | -16 | 6.046076609 |
| 25 | 8.747674251 | -25 | 0.428957106 | 25 | 19.13598192 | -25 | 2.109805814 |
| 2 4 1 | 20.40235305 | 2 4 -1 | 19.59722019 | 2 5 1 | 30.43614539 | 2 5 -1 | 29.56437148 |
| 4 | 21.60513304 | -4 | 18.38832874 | 4 | 31.74704320 | -4 | 28.26120113 |
| 9 | 23.58709333 | -9 | 16.38610602 | 9 | 33.93615107 | -9 | 26.10501394 |
| 16 | 26.29348662 | -16 | 13.67312110 | 16 | 36.99626751 | -16 | 23.12875359 |
| 25 | 29.62878614 | -25 | 10.51236733 | 25 | 40.89293269 | -25 | 19.38439052 |
........ and so on
........
b) Program#2
-The variant listed hereunder doesn't call any subroutine:
-The equation is re-written f(Lmn) = Lmn
and f(f(f(...(f(x))...)) is computed until the result is approximately
constant.
-Moreover, the number of terms in the continued fractions is fixed
to 5, which facilitates the calculations.
-On the other hand, the program is slower.
Data Registers: R00 = 10
R01 = m R02 = n R03 = c2
R04 = Lmn R05 to R09: temp
Flag: F01
Subroutines: /
-Line 09, 10 may be replaced by another even integer
-Line 113, this M-Code routine ( cf "A Few M-Code Routines for
the HP-41" ) may be replaced by X#Y?
-Line 114 is a three-byte GTO 01
| 01 LBL "LMN"
02 STO 03 03 RDN 04 STO 02 05 X<>Y 06 STO 01 07 - 08 STO 05 09 10 10 STO 00 11 SIGN 12 RCL 02 13 ST+ Y 14 * 15 RCL 03 16 2 17 / 18 + 19 STO 04 20 LBL 01 21 VIEW 04 22 RCL 05 23 RCL 00 |
24 STO 09
25 + 26 STO 06 27 SF 01 28 CLX 29 GTO 04 30 LBL 02 31 1 32 RCL 01 33 ST+ X 34 X^2 35 - 36 RCL 01 37 RCL 06 38 + 39 STO 08 40 ST* 08 41 ST+ 08 42 ST+ X 43 1 44 - 45 ST/ Y 46 4 |
47 +
48 / 49 RCL 03 50 2 51 / 52 ST* Y 53 + 54 RCL 08 55 + 56 + 57 RTN 58 LBL 03 59 XEQ 02 60 STO 07 61 RCL 05 62 RCL 00 63 STO 09 64 - 65 STO 06 66 CLX 67 LBL 04 68 XEQ 02 69 RCL 04 |
70 -
71 2 72 FC? 01 73 ST+ 06 74 SIGN 75 ST- 09 76 RCL 06 77 ST- Y 78 * 79 RCL 01 80 ST+ X 81 LASTX 82 + 83 ST* Y 84 1 85 - 86 ST* Y 87 RCL 06 88 + 89 ST/ Y 90 ST/ Y 91 2 92 FS? 01 |
93 ST- 06
94 + 95 ST/ Y 96 4 97 - 98 / 99 RCL 03 100 X^2 101 * 102 X<>Y 103 / 104 DSE 09 105 GTO 04 106 FS?C 01 107 GTO 03 108 RCL 04 109 X<>Y 110 RCL 07 111 + 112 STO 04 113 X#Y?? 114 GTO 01 115 END |
( 157 bytes / SIZE 010 )
| STACK | INPUTS | OUTPUTS |
| Z | m | / |
| Y | n | / |
| X | c2 | Lmn |
Example1: m = 4 ,
n = 11 , c2 = -1
4 ENTER^
11 ENTER^
-1 XEQ "LMN" >>>> Lmn
= 131.5600809
--- Execution time = 82s ---
-In this example, the 2nd routine is faster but it's seldom the case!
Example2: m = 0.2 , n = 0.6 , c2 = 1.7 ( actually, m & n are not necessarily integers )
0.2 ENTER^
0.6 ENTER^
1.7 R/S
>>>> Lmn = 2.246866650
---Execution time = 5mn46s ---
Note:
-Though R00 = 10 seems enough, the number of terms in the continued
fraction ( half the value of line 09 ) is not known in advance,
so you could press XEQ 01 after storing - say 12 - in
R00
-If the last decimal is ( almost ) the same, the result should be correct.
2°) Spheroidal Wave Functions (
-1 <= x <= +1 )
-The following programs compute the angular spheroidal wave function
of the
first kind by a series expansion in powers of x
-But first, given m , n and c2 , the corresponding
eigenvalue Lmn must be calculated by "LMN"
a) Program#1
Formulae:
Smn(x) = a0
+ a2 x2 + a4 x4 + ........
with a0 = [ (-1)(n+m)/2 (n+m)! ] / [ 2n
((n-m)/2)! ((n+m)/2)! ]
if n - m is even
Smn(x) = a1x
+ a3 x3 + a5 x5 + ........
with a1 = [ (-1)(n+m-1)/2 (n+m+1)! ] / [ 2n
((n-m-1)/2)! ((n+m+1)/2)! ]
if
n - m is odd
-In both cases, (k+1)(k+2) ak+2 + ( Lmn
- 2.k2 - m2 ) ak + ( k2 - 3.k
+ 2 - Lmn - c2 ) ak-2 + c2
ak-4 = 0
Data Registers: R00 = x R05 thru R09: temp ( Registers R01 thru R04 are to be initialized before executing "SMN" )
• R01 = m • R03 = c2
• R02 = n • R04 = Lmn
Flags: /
Subroutines: /
-Line 10 Y^X , an M-Code routine that gives 0^0 =
1 is preferable to avoid a DATA ERROR if x = (n-m) mod 2= 0.
-Otherwise, add X=Y? X#0? GTO 00 SIGN
X<>Y LBL 00 before this line
| 01 LBL "SMN"
02 STO 00 03 RCL 02 04 RCL 01 05 - 06 STO 07 07 2 08 MOD 09 STO 05 10 Y^X 11 STO 06 12 1 13 CHS 14 RCL 01 15 RCL 02 16 + 17 STO 09 18 RCL 05 19 ST- 07 |
20 ST+ 09
21 - 22 2 23 / 24 Y^X 25 RCL 09 26 FACT 27 * 28 2 29 ST/ 07 30 ST/ 09 31 RCL 02 32 Y^X 33 / 34 0 35 STO 08 36 X<> 07 37 FACT 38 RCL 09 |
39 FACT
40 * 41 / 42 STO 09 43 RCL 06 44 * 45 LBL 01 46 3 47 RCL 05 48 ST* Y 49 X^2 50 - 51 2 52 - 53 RCL 03 54 + 55 RCL 04 56 + 57 RCL 08 |
58 ST* Y
59 X<> 07 60 RCL 03 61 * 62 - 63 RCL 05 64 X^2 65 ST+ X 66 RCL 01 67 X^2 68 + 69 RCL 04 70 - 71 RCL 09 72 STO 08 73 * 74 + 75 RCL 05 76 1 |
77 +
78 ST/ Y 79 LASTX 80 + 81 STO 05 82 / 83 STO 09 84 RCL 00 85 X^2 86 RCL 06 87 * 88 STO 06 89 * 90 + 91 X#Y? 92 GTO 01 93 END |
( 113 bytes / SIZE 011 )
| STACK | INPUTS | OUTPUTS |
| Y | / | Smn(x) |
| X | x | Smn(x) |
Example1: m = n = 2 and
c2 = -25 "LMN" gives 0.4289571064
We store these 4 numbers into R01 thru R04 and
0.6 XEQ "SMN" >>>>
Smn(0.6) = 4.564797329 ( in 18 seconds
)
0.9
R/S >>>>
Smn(0.9) = 3.188333453 ( in 23 seconds
)
Example2: m = n = 0 and c2 = -16 "LMN" gives -9.150793382 We store these 4 numbers into R01 thru R04 and
0.7 XEQ "SMN" >>>>
Smn(0.7) = 4.557370657 ( in 18 seconds
)
1
R/S >>>>
Smn(1) = 12.41705490 ( in 19 seconds
)
Example3: m = 2 , n = 5 and c2 = 16 "LMN" gives 36.99626751 We store these 4 numbers into R01 thru R04 and
0.3 XEQ "SMN" >>>>
Smn(0.3) = -9.214845515 ( in 15 seconds
)
0.7
R/S >>>>
Smn(0.7) = 10.51929252
( in 18 seconds )
Notes:
-If m > 0 Smn(1) = Smn(-1) = 0 therefore, you could add ABS 1 X#Y? GTO 00 RCL 01 0 X<Y? RTN LBL 00 RCL 00 after line 02
-To avoid a premature stop, add FS?C 10
GTO 01 after line 92 and SF 10 after line 44.
-Practically, this seems highly improbable if c#0
b) Program#2
-Here, Smn(x) is computed under the form Smn(x)
= ( 1 - x2 )m/2 f(x)
where f(x) = Sumk=0,1,.... ak xk
with a0 = [ (-1)(n+m)/2
(n+m)! ] / [ 2n ((n-m)/2)! ((n+m)/2)! ]
a1 = 0
if n - m is even
and a1 = [ (-1)(n+m-1)/2
(n+m+1)! ] / [ 2n ((n-m-1)/2)! ((n+m+1)/2)! ]
a0 = 0
if n - m is odd
-In both cases, (k+1)(k+2) ak+2 - [ k ( k + 2m
+ 1 ) - Lmn + m ( m + 1 ) ] ak - c2
ak-2 = 0
Data Registers: R00 = x R05 thru R08: temp ( Registers R01 thru R04 are to be initialized before executing "SMN" )
• R01 = m • R03 = c2
• R02 = n • R04 = Lmn
Flags: /
Subroutines: /
-Lines 10 & 86 Y^X , an M-Code routine that gives
0^0 = 1 is preferable to avoid a DATA ERROR.
-Otherwise, add X=Y? X#0? GTO 00 SIGN
X<>Y LBL 00 before these lines
| 01 LBL "SMN"
02 STO 00 03 RCL 02 04 RCL 01 05 - 06 STO 07 07 2 08 MOD 09 STO 05 10 Y^X 11 STO 06 12 1 13 CHS 14 RCL 01 15 RCL 02 16 + 17 STO 08 18 RCL 05 |
19 ST- 07
20 ST+ 08 21 - 22 2 23 / 24 Y^X 25 RCL 08 26 FACT 27 * 28 2 29 ST/ 07 30 ST/ 08 31 RCL 02 32 Y^X 33 / 34 0 35 X<> 07 36 FACT |
37 RCL 08
38 FACT 39 * 40 / 41 STO 08 42 RCL 06 43 * 44 LBL 01 45 RCL 03 46 RCL 07 47 * 48 RCL 01 49 ST+ X 50 RCL 05 51 ST+ Y 52 ST* Y 53 + 54 RCL 01 |
55 ST+ Y
56 X^2 57 + 58 RCL 04 59 - 60 RCL 08 61 STO 07 62 * 63 + 64 RCL 05 65 1 66 + 67 ST/ Y 68 LASTX 69 + 70 STO 05 71 / 72 STO 08 |
73 RCL 00
74 X^2 75 RCL 06 76 * 77 STO 06 78 * 79 + 80 X#Y? 81 GTO 01 82 RCL 00 83 ASIN 84 COS 85 RCL 01 86 Y^X 87 * 88 END |
( 109 bytes / SIZE 009 )
| STACK | INPUT | OUTPUT |
| X | x | Smn(x) |
Example: m = n = 2 and
c2 = -25 "LMN" gives 0.428957107
We store these 4 numbers into R01 thru R04 and
0.6 XEQ "SMN" >>>>
Smn(0.6) = 4.564797329 ( in 16 seconds
)
0.9
R/S >>>>
Smn(0.9) = 3.188333454 ( in 17 seconds
)
Note:
-This variant is often faster than the first version.
c) A Generalization
-Many normalization schemes are possible when m and n are not integers.
-The routine listed hereafter employs:
Smn(x) = ( 1 - x2 )m/2 f(x) where f(x) = Sumk=0,1,.... ak xk
with a0 = Pnm(0)
= 2m sqrt(PI) / [ Gam((1-m-n)/2) Gam((2-m+n)/2 ]
( Flammer's scheme )
and a1 = P'nm(0)
= ( m + n ) 2m sqrt(PI) / [ Gam((2-m-n)/2) Gam((1-m+n)/2 ]
-In both cases, (k+1)(k+2) ak+2 - [ k ( k + 2m + 1 ) - Lmn + m ( m + 1 ) ] ak - c2 ak-2 = 0
-It gives the same results as above when m & n are non-negative
integers, with n = m , m+1 , m+2 , ......
Data Registers: R00 = x R05 thru R12: temp ( Registers R01 thru R04 are to be initialized before executing "SWF" )
• R01 = m • R03 = c2
• R02 = n • R04 = Lmn
Flags: /
Subroutine: "GAM" or "GAM+"
... ( cf "Gamma Function for the HP-41" )
-Line 110 Y^X , an M-Code routine that gives 0^0 =
1 is preferable to avoid a DATA ERROR if x = 1 and m = 0
-Otherwise, add X=Y? X#0? GTO 02 SIGN
X<>Y LBL 02 after line 109
| 01 LBL "SWF"
02 STO 00 03 2 04 RCL 01 05 Y^X 06 PI 07 SQRT 08 * 09 STO 09 10 STO 11 11 SIGN 12 STO 07 13 RCL 01 14 RCL 02 15 + 16 ST* 11 17 - 18 2 19 / 20 STO 05 21 XEQ "GAM" 22 ST/ 09 23 RCL 02 |
24 RCL 01
25 - 26 2 27 ST+ Y 28 / 29 STO 06 30 XEQ "GAM" 31 ST/ 09 32 RCL 05 33 .5 34 ST- 06 35 + 36 XEQ "GAM" 37 ST/ 11 38 RCL 06 39 XEQ "GAM" 40 ST/ 11 41 CLX 42 STO 06 43 STO 08 44 STO 10 45 RCL 09 46 RCL 00 |
47 RCL 11
48 * 49 + 50 STO 05 51 GTO 01 52 LBL 00 53 RCL 11 54 X<> 09 55 STO 11 56 RCL 10 57 X<> 08 58 STO 10 59 RCL 03 60 * 61 RCL 06 62 RCL 01 63 ST+ X 64 + 65 RCL 06 66 ST* Y 67 + 68 RCL 04 69 - |
70 RCL 01
71 ST+ Y 72 X^2 73 + 74 R^ 75 STO 10 76 * 77 + 78 RCL 06 79 1 80 + 81 STO 06 82 ST/ Y 83 LASTX 84 + 85 / 86 STO 11 87 RCL 07 88 RCL 00 89 ST* 07 90 X^2 91 * 92 * |
93 RTN
94 LBL 01 95 XEQ 00 96 STO 12 97 XEQ 00 98 RCL 12 99 + 100 RCL 05 101 + 102 STO 05 103 LASTX 104 X#Y? 105 GTO 01 106 RCL 00 107 ASIN 108 COS 109 RCL 01 110 Y^X 111 * 112 END |
( 158 bytes / SIZE 013 )
| STACK | INPUT | OUTPUT |
| X | x | Smn(x) |
Example: m = 0.2
n = 0.6 c2 = 1.7
"LMN" gives 2.246866650
-Store these 4 numbers into R01 thru R04 and
0.7 XEQ "SWF" >>>>
Smn(0.7) = 0.682645661
---Execution time = 91s---
3°) Other Normalizations
a) Normalization Factors
"SHNF" calculates f = 1 / sqrt [ §-1+1 { S(x) }2 dx ] where S(x) is non-normalized,
i-e S(0) = 1 , S'(0) = 0
if n-m is even
and S(0) = 0 , S'(0) = 1 if n-m is odd
-It may use many registers if the convergence is slow ( the series described
in paragraph 2°) a) are integrated )
Data Registers: R00 = partial sums R05 thru ............ : temp ( Registers R01 thru R04 are to be initialized before executing "SHNF" )
• R01 = m • R03 = c2
• R02 = n • R04 = Lmn
at the end, R05 = f = normalization factor.
Flags: /
Subroutines: /
| 01 LBL "SHNF"
02 RCL 02 03 RCL 01 04 - 05 2 06 MOD 07 STO 05 08 STO N 09 ST+ X 10 1 11 + 12 1/X 13 STO 00 14 CLX 15 STO 07 16 STO 08 17 SIGN 18 STO 09 19 9 20 STO 06 21 .1 22 % |
23 +
24 STO O 25 LBL 01 26 RCL 06 27 2 28 - 29 SIGN 30 ST- 06 31 E-3 32 ST+ O 33 RCL 03 34 RCL IND L 35 * 36 RCL 05 37 X^2 38 LASTX 39 3 40 * 41 - 42 2 43 + 44 RCL 03 |
45 -
46 RCL 04 47 - 48 RCL IND 06 49 * 50 + 51 RCL 04 52 RCL 05 53 X^2 54 ST+ X 55 - 56 RCL 01 57 X^2 58 - 59 ISG 06 60 CLX 61 RCL IND 06 62 * 63 + 64 CHS 65 RCL 05 66 1 |
67 +
68 ST/ Y 69 LASTX 70 ST+ 06 71 + 72 STO 05 73 / 74 STO IND 06 75 RCL O 76 RCL 06 77 STO M 78 CLX 79 LBL 02 80 RCL IND Y 81 RCL IND M 82 * 83 + 84 DSE M 85 ISG Y 86 GTO 02 87 RCL 05 88 1 |
89 +
90 RCL N 91 + 92 / 93 RCL 00 94 + 95 STO 00 96 VIEW X 97 LASTX 98 X#Y? 99 GTO 01 100 + 101 SQRT 102 1/X 103 STO 05 104 CLA 105 CLD 106 END |
( 141 bytes / SIZE ? )
| STACK | INPUTS | OUTPUTS |
| X | / | f |
Example1: m = 0 , n = 1
, c2 = 2 "LMN" gives
3.172127920 We store these 4
numbers into R01 thru R04 and
XEQ "SHNF" >>>> f = 1.376455541 = R05 ( in 34 seconds )
Example2: m = n = 2 and c2 = 3 "LMN" gives 6.412006610 We store these 4 numbers into R01 thru R04 and
XEQ "SHNF" >>>> f = 0.9962441297 ( in 51 seconds )
Notes:
-The successive partial sums are displayed ( line 96 ).
-If they oscillate between large numbers of opposite signs,
the last decimals may be meaningless !
-For instance, with m = 4 , n = 11 , c2 = -1 , Lmn
= 131.5600809 in registers R01 thru R04, "SHNF"
displays numbers of the order of +/- 3000
-So the final result f = 8.879731920 probably has
no more than 6 correct decimals...
b) Normalized Functions
-The results obtained just above may now be used to compute the angular
spheroidal wave function of the first kind, so that §-1+1
[ S(x) ]2 dx = 1
-We simply have to multiply the non-normalized function by the factor
f given by "SHNF"
-Here, Smn(x) is computed under the form Smn(x)
= ( 1 - x2 )m/2 f(x)
where f(x) = Sumk=0,1,.... ak xk
with a0 = 1
, a1 = 0
if n - m is even
and a0 = 0
, a1 = 1
if n - m is odd
-In both cases, (k+1)(k+2) ak+2 - [ k ( k + 2m + 1 ) - Lmn + m ( m + 1 ) ] ak - c2 ak-2 = 0
-Finally, the result is multiplied by f.
Data Registers: R00 = x R06 thru R10: temp ( Registers R01 thru R05 are to be initialized before executing "SWF2" )
• R01 = m • R03 = c2
• R05 = f
• R02 = n • R04 = Lmn
Flags: /
Subroutines: /
-An M-Code routine Y^X that gives 0^0 = 1 is preferable
to avoid a DATA ERROR if X = Y = 0 ( lines 09 & 61 )
| 01 LBL "SWF2"
02 STO 00 03 RCL 02 04 RCL 01 05 - 06 2 07 MOD 08 STO 07 09 Y^X 10 STO 06 11 STO 08 12 CLX 13 STO 09 |
14 SIGN
15 STO 10 16 LBL 01 17 RCL 03 18 RCL 09 19 * 20 RCL 01 21 ST+ X 22 RCL 07 23 ST+ Y 24 ST* Y 25 + 26 RCL 01 |
27 ST+ Y
28 X^2 29 + 30 RCL 04 31 - 32 RCL 10 33 STO 09 34 * 35 + 36 RCL 07 37 1 38 + 39 ST/ Y |
40 LASTX
41 + 42 STO 07 43 / 44 STO 10 45 RCL 08 46 RCL 00 47 X^2 48 * 49 STO 08 50 * 51 RCL 06 52 + |
53 STO 06
54 LASTX 55 X#Y? 56 GTO 01 57 RCL 00 58 ASIN 59 COS 60 RCL 01 61 Y^X 62 * 63 RCL 05 64 * 65 END |
( 80 bytes / SIZE 011 )
| STACK | INPUT | OUTPUT |
| X | x | Smn(x) |
Example1: m = 0 , n = 1
, c2 = 2 , Lmn = 3.172127920
, f = 1.376455541 These
5 numbers are in R01 thru R05 after executing "SHNF"
0.4 XEQ "SWF2" >>>> S(0.4) = 0.533565783 ---Execution time = 8s---
Example2: m = n = 2 , c2 = 3 , Lmn = 6.412006610 , f = 0.9962441297 these 5 numbers must be in R01 thru R05
0.4 XEQ "SWF2" >>>>
S(0.4) = 0.809618196
---Execution time = 8s---
Note:
-If m > 0 Smn(1) = Smn(-1) = 0
so you could add ABS 1 X#Y? GTO 00 RCL 01
0 X<Y? RTN LBL 00 RCL 00 after line
02
"Meixner-Shaefke" Scheme:
-It requires that §-1+1 [ S(x) ]2 dx = [ 2/(2n+1) ] (m+n)! / (n-m)!
-Simply add RCL 01 RCL 02 + FACT RCL 02 RCL 01 - FACT / ST+ X RCL 02 ST+ X 1 + / SQRT * after line 64
-In the last example, it yields S(0.4) = 2.508510232
References:
[1] Abramowitz & Stegun , "Handbook of Mathematical
Functions" - Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/