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