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
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall