hp41programs

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:

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