The Museum of HP Calculators

# Spheroidal Harmonics for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) Eigenvalues
2°) Spheroidal Wave Functions ( | x | <= 1 )

-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  ------------------  oblate   -------------------------           ( 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

-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"                               ( the same character as the LBL line 27 )
11  1
12  +
13  *
14  RCL 10
15  2
16  /
17  +
18  ENTER^
19  ASTO 20
20  1
21  +
22  "U"                               ( the same character as the LBL line 43 )
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 )

 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 ........

2°) Spheroidal Wave Functions    ( -1 <= x <= +1 )

-The following program computes 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"

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          R06 thru R10: temp                         ( Registers R01 thru R04 are to be initialized before executing "SMN" )

•  R01 = m      •  R03 = c       R05 = partial sums and finally  Smn(x)
•  R02 = n       •  R04 = Lmn
Flag:   F01
Subroutines: /

01  LBL "SMN"
02  STO 00
03  SF 01
04  RCL 02
05  RCL 01
06  -
07  2
08  MOD
09  X=0?
10  CF 01
11  ST+ 01                                R01 is modified if n-m is odd, but its original content is restored line 51
12  STO 06
13  RCL 00
14  1
15  FS? 01
16  X<>Y
17  STO 05
18  STO 07
19  1
20  CHS
21  RCL 02
22  RCL 01
23  -
24  2
25  /
26  STO 08
27  Y^X
28  RCL 01
29  RCL 02
30  +
31  FACT
32  *
33  RCL 08
34  FACT
35  RCL 02
36  LASTX
37  -
38  FACT
39  *
40  2
41  RCL 02
42  Y^X
43  *
44  /
45  ST* 05
46  STO 10
47  CLX
48  STO 08
49  STO 09
50  FS?C 01
51  DSE 01
52  LBL 01
53  RCL 03
54  RCL 08
55  *
56  RCL 06
57  3
58  -
59  RCL 06
60  *
61  2
62  +
63  RCL 03
64  -
65  RCL 04
66  -
67  RCL 09
68  STO 08
69  *
70  +
71  RCL 04
72  RCL 06
73  X^2
74  ST+ X
75  -
76  RCL 01
77  X^2
78  -
79  RCL 10
80  STO 09
81  *
82  +
83  CHS
84  RCL 06
85  1
86  +
87  ST/ Y
88  LASTX
89  +
90  STO 06
91  /
92  STO 10
93  RCL 07
94  RCL 00
95  X^2
96  *
97  STO 07
98  *                                              The program stops when 2 consecutive sums are equal.
99  RCL 05                                    To avoid a premature interruption, add  DSE 11   GTO 02  after line 104
100  +                                              and  3  STO 11   LBL 02  after line 52
101  STO 05                                    in case 1 or 2 successive ak = 0  and the next one is not.
102  LASTX                                    Practically, this seems highly improbable if c#0
103  X#Y?
104  GTO 01
105  END

( 123 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 19 seconds )
0.9         R/S          >>>>     Smn(0.9) = 3.188333453   ( in 24 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.557370593   ( in 19 seconds )
1           R/S          >>>>       Smn(1)  = 12.41705452   ( in 20 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 19 seconds )

Note:

-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   after line 02

Reference:      Abramowitz & Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4