# hp41programs

Ellipsoidal Wave Function

# Ellipsoidal Wave Functions for the HP-41

Overview

0°)  Jacobian Elliptic Function Sn(x)

a) Focal Program
b) M-Code Routine

1°)  Ellipsoidal Wave Differential Equation
2°)  Ellipsoidal Wave Function

0°)  Jacobian Elliptic Function Sn ( x | k^2 ) k < 1

a)  Focal Program

-The ellipsoidal wave function W(x) involves Sn(x)
-So I wrote a specific program to calculate this function.

Data Registers:       R19 to R27: temp
Flags: /
Subroutines: /

 01 LBL "SNX"  02 DEG  03 20.02  04 STO 19   05 STO 20  06 RDN  07 LBL 01           08 1  09 RCL Z  10 -  11 SQRT  12 1 13 X<>Y  14 -  15 X=0?  16 GTO 00  17 1  18 ST+ 20   19 LASTX  20 +  21 /  22 STO IND 20  23 X^2  24 X<>Y 25 LASTX  26 1  27 +  28 /  29 GTO 01  30 LBL 00           31 X<>Y  32 R-D  33 SIN  34 SIGN  35 RCL 19   36 RCL 20 37 X#Y?  38 GTO 00  39 LASTX  40 RTN  41 LBL 00  42 LASTX  43 LBL 02  44 ENTER  45 X^2  46 RCL IND 20  47 *  48 1 49 +  50 X<>Y  51 RCL IND 20  52 LASTX  53 +  54 *  55 X<>Y  56 /  57 DSE 20  58 GTO 02  59 END

( 85 bytes / SIZE 028 )

 STACK INPUTS OUTPUTS Y k^2 / X x Sn ( x | k^2 )

Example:

0.8   ENTER^
0.7   XEQ "SNX"  >>>>  0.612365484                ---Execution time = 7s---

Notes:

-Though it will not work if k = 1, it will with k2 = .9999999999

.9999999999  ENTER^
0.7            R/S        >>>>     0.604367777   which is actually  Tanh 0.7

-With k = 0 and x = 0.7  again,  it yields  0.644217687 = Sin 0.7

-Registers R19 thru R27 are used

b)  M-Code Routine

-This M-Code routine has an advantage: it doesn't use any data register.

098 "X"
00E "N"
013 "S"
04E  C=0 ALL
028  T=C
0B8  C=Y
2FE  ?C#0 MS
0B5  ?CGO                     DATA ERROR
0A3  DATAERROR        if  m < 0
228  P=C
10E  A=C ALL
04E  C
35C  =
050   1
36E  ?A#C ALL
0B5  ?NCGO                   DATA ERROR
0A2  DATAERROR         if m = 1
0F8  C=X
070  N=C ALL
2A0  SETDEC
238  C=P
2BE  C=-C
02E  C
0FA ->
0AE  AB
001  C=
060  AB+1
305  C=
060  sqrt(AB)
268  Q=C
1BE  AB=-AB
001  C=
061  AB+1
2EE  ?C#0 ALL
133  JNC+38d
1E8  O=C
278  C=Q
02E  C
0FA ->
0AE  AB
001  C=
060  AB+1
10E  A=C ALL
1F8  C=O
0AE  A<>C ALL
261  C=
060  A/C
1E8  O=C
001  C=
060  AB+1
10E  A=C ALL
0B0  C=N ALL
0AE  A<>C ALL
261  C=
060  A/C
070  N=C ALL
1F8  C=O
10E  A=C ALL
135  C=
060  A*C
228  P=C
1F8  C=O
10E  A=C ALL
046  C
270  =
038  T
226  C=C+1 S&X
028  T=C
270  RAMSLCT
0AE  A<>C ALL
2F0  WRITDATA
26B  JNC-51d
130  LDI S&X
090  090h
358  ST=C XP
144  CLRF 6
284  CLRF 7
0B0  C=N ALL
229  C=
048  Sin(C)
070  N=C ALL
046  C
270  =
038  T
2E6  ?C#0 S&X
123  JNC+36d
0B0  C=N ALL
10E  A=C ALL
135  C=
060  A*C
046  C
270  =
038  T
270  RAMSLCT
1E8  O=C
13D  C=
060  AB*C
001  C=
060  AB+1
268  Q=C
1F8  C=O
02E  C
0FA  ->
0AE  AB
001  C=
060  AB+1
0B0  C=N ALL
13D  C=
060  AB*C
278  C=Q
269  C=
060  AB/C
070  N=C ALL
046  C
270  =
038  T
266  C=C-1 S&X
028  T=C
2E6  ?C#0 S&X
2F7  JNC-34d
0B0  C=N ALL
0E8  X=C
345  ?NCGO
042  CLA

( 125 words / SIZE 000 )

 STACK INPUTS OUTPUTS Y k^2 < 1 / X x Sn ( x | k^2 )

Example:

0.8   ENTER^
0.7   XEQ "SNX"  >>>>  0.612365484                ---Execution time = 2.5s---

Notes:

-Though it will not work if k = 1, it will with k2 = .9999999999

.9999999999  ENTER^
0.7         XEQ "SNX"  >>>>     0.604367777   which is actually  Tanh 0.7

-With k = 0 and x = 0.7  again,  it yields  0.644217687 = Sin 0.7

-This routine uses registers Z Y X L M N O instead of R01 .............. R07
-On the other hand, m & x are not saved.

-There will be a DATA ERROR message if m < 0 or m = 1 or m > 1
-There is no check for alpha data, overflow or underflow.

-The alpha regiser is cleared and synthetic register Q is used.

1°)  Ellipsoidal Wave Differential Equation ( RK6 )

"EWF" uses a Runge-Kutta-Nystrom method of order 6 to solve the differential equation:

W''(x) = W(x) [ -h + n(n+1) k2 Sn2 (x) - q k4 Sn4 (x) ]     where  Sn (x) is a Jacobian elliptic function.

Data Registers:              R00 = temp          ( Registers R01 thru R09 are to be initialized before executing "EWF" )

•  R01 = h                •  R05 = x0           •  R08 = H = stepsize          R10 to R14: temp         R19 to R27: used by "SNX"
•  R02 = n                •  R06 = W(x0)     •  R09 = N = nb of steps     R15 to R18: unused
•  R03 = k2 < 1        •  R07 = W'(x0)
•  R04 = q
Flags: /
Subroutine:  "SNX"

 01 LBL "EWF"  02 RCL 09  03 STO 13            04 LBL 13  05 RCL 06  06 RCL 05  07 XEQ 10  08 STO 00  09 RCL 08  10 *  11 8  12 /  13 RCL 07  14 +  15 RCL 08   16 *  17 4  18 /  19 RCL 06  20 +  21 RCL 08  22 4  23 /  24 RCL 05  25 + 26 XEQ 10  27 STO 10  28 6  29 /  30 RCL 00            31 24  32 /  33 -  34 RCL 08  35 *  36 RCL 07  37 2  38 /  39 +  40 RCL 08   41 *  42 RCL 06  43 +  44 RCL 08  45 2  46 /  47 RCL 05  48 +  49 XEQ 10  50 STO 11 51 ST+ X  52 RCL 10  53 4  54 *  55 +  56 RCL 00            57 3  58 *  59 +  60 8  61 /  62 RCL 08  63 *  64 RCL 07  65 3  66 *  67 +  68 RCL 08  69 *  70 4  71 /  72 RCL 06  73 +  74 RCL 08  75 .75 76 *  77 RCL 05  78 +  79 XEQ 10  80 STO 12   81 ST+ X  82 RCL 11  83 -  84 RCL 10            85 6  86 *  87 +  88 14  89 /  90 RCL 08  91 *  92 RCL 07  93 +  94 RCL 08  95 ST+ 05  96 *  97 RCL 06  98 +  99 RCL 05 100 XEQ 10 101 RCL 00 102 + 103 7 104 * 105 RCL 10  106 RCL 12 107 + 108 32 109 * 110 + 111 RCL 11          112 12 113 * 114 + 115 90 116 / 117 RCL 08 118 * 119 RCL 12 120 8 121 * 122 RCL 11 123 6 124 * 125 + 126 RCL 10 127 24 128 * 129 + 130 RCL 00  131 7 132 * 133 + 134 90 135 / 136 RCL 08          137 * 138 RCL 07 139 + 140 RCL 08 141 * 142 ST+ 06 143 X<>Y 144 ST+ 07 145 DSE 13 146 GTO 13 147 GTO 12 148 LBL 10 149 X<>Y 150 STO 14 151 RCL 03 152 RCL Z 153 XEQ "SNX" 154 X^2 155 RCL 03  156 * 157 STO Y 158 RCL 04 159 * 160 RCL 02 161 ENTER 162 ST* Y 163 + 164 X<>Y 165 - 166 * 167 RCL 01          168 - 169 RCL 14 170 * 171 RTN 172 LBL 12 173 RCL 07 174 RCL 06 175 RCL 05 176 END

( 219 bytes / SIZE 028 )

 STACK INPUTS OUTPUTS Z / W'(x0 + N.H) Y / W(x0 + N.H) X / x0 + N.H

Example:    h = 1.2    n = 1.7   k2 = 0.8   q = sqrt(2)     W(0) = 1   W'(0) = 0

1.2  STO 01                       0  STO 05
1.7  STO 02                       1  STO 06
0.8  STO 03                       0  STO 07
2   SQRT  STO 04

-If you want to compute W(1) with  H = 0.1 and therefore,  N = 10

0.1  STO 08
10   STO 09

XEQ "EWF"  >>>>      X   =  1
RDN   W(1) =  0.640627307                                    ---Execution time = 7m15s---
RDN  W'(1) = -0.401079727

-Simply press R/S to continue with the same parameters, it yields:

W(2) = 0.542730302
W'(2) = 0.260776593

Note:

-Though "EWF" doesn't work if k2 = 1 , it will work with k2 = .9999999999 and the difference will be negligible.

2°)  Ellipsoidal Wave Function ( Series )

-"EWF2" is less general than "EWF".
-The solution is of the form  W(x) = 1 + A1 t + A2 t2 + ............. + Am tm + ................

where  t = Sn2 (x)

-This leads to the recurrence relation

( 2m+1 ) (2m+2 ) Am+1 = - q k4 Am-2 + Am-1 k2 [ n(n+1) - 2 ( m-1 ) ( 2m - 1 ) ] + Am [ - h + 4 m2 ( k2 + 1 ) ]

Data Registers:              R00 = x               ( Registers R01 thru R04 are to be initialized before executing "EWF" )

•  R01 = h                R05 to R13: temp
•  R02 = n
•  R03 = k2 < 1        R10 = t = Sn2 (x)   R12 = Sum
•  R04 = q
Flags: /
Subroutine:  "SNX"

 01 LBL "EWF2"  02 STO 00  03 1  04 STO 08  05 STO 12  06 RCL 03  07 RCL 00  08 XEQ "SNX"  09 X^2  10 STO 10   11 CLX  12 STO 05           13 STO 06  14 STO 07  15 SIGN  16 STO 09  17 3 18 STO 11  19 LBL 03  20 RCL 11  21 STO 13  22 LBL 04  23 RCL 03  24 X^2  25 RCL 04   26 *  27 RCL 06  28 *  29 RCL 05           30 ENTER  31 ST+ Y  32 1  33 ST- Z  34 - 35 *  36 ST+ X  37 RCL 02  38 ENTER  39 ST* Y  40 +  41 X<>Y  42 -  43 RCL 03   44 *  45 RCL 07  46 STO 06           47 *  48 X<>Y  49 -  50 RCL 03  51 RCL 05 52 ST+ X  53 X^2  54 ST* Y  55 +  56 RCL 01  57 -  58 RCL 08  59 STO 07   60 *  61 +  62 RCL 05           63 1  64 +  65 STO 05  66 ST+ X  67 ENTER  68 DSE X 69 *  70 /  71 STO 08  72 RCL 09  73 RCL 10  74 *  75 STO 09   76 *  77 RCL 12  78 +  79 STO 12           80 LASTX  81 X#Y?  82 GTO 03  83 DSE 13  84 GTO 04  85 END

( 109 bytes / SIZE 028 )

 STACK INPUTS OUTPUTS X x W(x)

Example:             h = 1.2    n = 1.7   k2 = 0.8   q = sqrt(2)

1.2  STO 01
1.7  STO 02
0.8  STO 03
2   SQRT  STO 04    and to calculate  W(1)

1   XEQ "EWF2"  >>>>   W(1) = 0.640627308                                    ---Execution time = 79s---

Notes:

-In this example, "EWF2" is faster than "EWF", but it's not always the case...
-For instance, W(2) = 0.542730301  is returned after a very long time.

-Though "EWF" doesn't work if k2 = 1 , it will work with k2 = .9999999999 and the difference will be negligible.