SphHarm

# 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:

   Abramowitz & Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
   http://functions.wolfram.com/