hp41programs

Ellipsoidal Wave Function

Ellipsoidal Wave Functions for the HP-41


Overview
 

 0°)  Jacobian Elliptic Function Sn(x)
 1°)  Ellipsoidal Wave Differential Equation
 2°)  Ellipsoidal Wave Function
 
 

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

-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 = Cos 0.7

-Registers R19 thru R27 are 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.