Coulomb Wave Functions for the HP-41

Overview
 

 1°) Regular Coulomb Wave Function

   a)  Program#1
   b)  Program#2
   c)  Non-integer values of L

 2°) Irregular Coulomb Wave Function

   a)  Program#1
   b)  Program#2
   c)  Non-integer values of L

 3°) Asymptotic Expansion
 

-The regular Coulomb wave function  F_L(n,r)  and the irregular Coulomb wave function G_L(n,r) are 2 independant solutions  of the differential equation:

           d²y/dr² + [ 1 - 2n/r - L(L+1)/r² ] y = 0
 

1°) Regular Coulomb Wave Function
 

       a)  Program#1
 

Formulae:         F_L(n,r) = C_L(n) r ^L+1 Sum _k>L A_k^L (n) r ^k-L-1

    with   C_L(n) = (1/Gam(2L+2)) 2^L e ^-pi.n/2 | Gam(L+1+i.n) |
    and    A_L+1^L = 1 ;  A_L+2^L = n/(L+1)  ;  (k+L)(k-L-1) A_k^L = 2n A_k-1^L - A_k-2^L   ( k > L+2 )
 

Data Registers:   R00 thru R05: temp

                              R06 = r  ;  R07 = n  ;  R08 = L  ;  R09 = C_L(n)
Flags: /
Subroutine:   "GAMZ"  ( Gamma Function in the complex plane )
 
 

 
   ( 96 bytes / SIZE 010 )
 
 

      L is a non-negative integer,
      n  is real,
      r  is non-negative.

Example:

     2   ENTER^
   0.7  ENTER^
   1.8  XEQ "RCWF"  >>>>   F₂( 0.7 , 1.8 ) = 0.141767744   ( in 34 seconds )
 

       b)  Program#2
 

-This program uses the same basic formulae but  | Gam( L+1+i n ) |  is calculated without  "GAMZ"

-Instead, we have also   | Gam( 1+i y ) |² = (Pi.y) / Sinh (Pi y)   and

    | Gam( 1+L+i y ) |² = [ L² + y² ] [ (L-1)² + y² ] .................. [ 1 + y² ]  (Pi.y) / Sinh (Pi y)
 

Data Registers:   R01 = r  ;  R02 = n  ;  R03 = L  ;  R00 & R04 thru R06: temp
Flag:  F10
Subroutines:  /

-Line 30 is a synthetic  TEXT0  or another  NOP  like  LBL 00
 
 

 
   ( 125 bytes / SIZE 007 )
 
 

       L is a non-negative integer,
       n  is real,
       r  is non-negative.

Example:

     2   ENTER^
   0.7  ENTER^
   1.8  XEQ "RCWF2"  >>>>   F₂( 0.7 , 1.8 ) = 0.141767746                ---Execution time = 17s---

-"RCWF2" is both faster and more accurate than "RCWF"
 

       c)  Non-integer Values of L
 

Formula:
 

   F_L(n,r) = (1/2) [ 1/(2L+1)! ] | Gam(L+1+i.n) | exp [ -(PI) n/2 - i (L+1)(PI)/2 ]  M_i.n,L+1/2 (2i.r)

   where  M_k,µ = Whittaker's function of the 1st kind
 

-Since it involves complex numbers, we cannot use the program "WHIM" listed in "Whittaker's Function for the HP-41"
 

Data Registers:   R00 thru R11: temp

                              R06 = r  ;  R07 = n  ;  R08 = L
Flags: /
Subroutines:  "GAM3"  or another "gamma-function program"  ( cf "Gamma Function for the HP-41" )
                         "GAMZ+" or "GAMZ"    ( gamma function in the complex plane )
 

-The M-Code routine  Z*Z  may be replaced by  XEQ "Z*Z"  ( cf "Complex Functions for the HP-41" )
-Lines 12-13-14 & 17-18  are only useful to compute the irregular Coulomb wave function with  "ICWF3"
-Otherwise, they may be deleted.
 
 

 
   ( 134 bytes / SIZE 012 )
 
 

 
Example:

   1.4   ENTER^
  -2.1   ENTER^
   1.6   XEQ "RCWF3"  >>>>   F_1.4( -2.1 , 1.6 ) = -0.779739876                ---Execution time = 43s---

-If there were no roundoff errors, Y-output = 0
-Actually, the small number in Y-register gives an idea of the accuracy of the results:
-If, for instance, Y = E-4 , there are probably no more than 4 exact decimals.
-Here, Y = -2 E-9 , so the result is probably correct to 8D ( an HP-48 returns -0.779739873 )

-This program also works if L is an integer.
 

2°) Irregular ( Logarithmic ) Coulomb Wave Function
 

       a)  Program#1
 

Formulae:         G_L(n;r) = (1/Pi).(exp(2.Pi.n) - 1) F_L(n;r) [ Ln(2r) + q_L(n)/p_L(n) ] + (1/C_L(n)/r^L/(2L+1))  Sum _k>-L-1  a_k^L (n) r ^k+L

    with      a _-L-1^L = 0 ;  a _-L^L = 1  ;  (k+L)(k-L-1) a_k^L = 2n a_k-1^L - a_k-2^L - (2k-1).p_L(n).A_k^L    and   a _L+1^L = 0

                p_L(n) = 2n ( 1+n² )( 4+n² ) ........... ( L²+n² ) 2^2L / ( 2L+1) / [(2L)!]²

       q_L(n)/p_L(n) = Sum _s=1,...,L s/(s²+n²)  - 1/1 - 1/2 - 1/3 - ...... - 1/(2L+1) + Ré ( Psi(1+i.n) ) + 2.gamma + r_L(n)/p_L(n)   ( gamma = Euler's Constant )

  where    r_L(n) = ( (-1)^L+1 / (2L)! )  Im [ 1/(2L+1) + 2(i.n-L)/1!/(2L) + 2²(i.n-L)(i.n-L+1)/2!/(2L-1)
                                                               + ...................................... + 2^2L(i.n-L)(i.n-L+1) ....... (i.n+L-1)/(2L)! ]
 

Data Registers:   R00 thru R05 & R11 thru R15: temp

                              R06 = r  ;  R07 = n  ;  R08 = L  ;  R09 = C_L(n)  ;  R10 = F_L(n;r)
Flags: /
Subroutines:  "RCWF"  ( regular Coulomb wave function )
                         "PSIZ"    ( Digamma function in the complex plane )

-Line 106 is 2 x Euler's Constant
-Line 178 is  TEXT 0  or another  NOP  instruction like STO X
 
 

 
   ( 259 bytes / SIZE 016 )
 
 

  L is a non-negative integer,  n  is real  and  r  is positive.      F_L(n,r) is stored in R10

Example:

            3   ENTER^
         -0.4  ENTER^
          1.2  XEQ "ICWF"  >>>>  G₃( -0.4 , 1.2 ) = 6.563265438        ( in 102 seconds )
                               and we have  F₃( -0.4 , 1.2 ) = 0.029547268          in register R10.

-Unfortunately, several significant digits are sometimes cancelled when the last operation ( line 199 ) is performed.
 ( Check the value in LASTX or in R11 )
 

       b)  Program#2
 

-Here,  Ré [ Psi( 1+i.x) ]  is computed without calling "PSIZ"

-We use:  Ré [ Psi( 1+i.x) ] = -gamma + x² SUM_{k=1,2,.........} 1/[k(k²+x²)]                              ( gamma = Euler's Constant )

-This infinite sum is actually calculated with the help of Euler-Mac Laurin formula, and after some algebra it yields:

 Ré [ Psi( 1+i.x) ] =  - SUM_{k=1,2,.....15}  k/(k²+x²) + (1/2) Ln ( 16²+x² ) - (16/2)( 16²+x² ) ^-1  + (1/12)( 16²+x² ) ^-2 [ 0.1 + x² - 16² - (24/30) 16² x² ( 16²+x² ) ^-2 ]
                               + eps

  where  eps is a very small number of the order of  -2.9 10 ^-10
 

Data Registers:   R01 = r  ,  R02 = n  ,  R03 = L  ,  R04 = F_L(n,r)    R00 & R05 thru R15: temp
Flag:   F10
Subroutine:  "RCWF2"    ( § 1°) b) )

-You could perhaps add   3 E-10  -   after line 58  though it dosn't seem to change the results significantly.
 
 

 
   ( 313 bytes / SIZE 016 )
 
 

    L is a non-negative integer,
    n  is real,
    r  is positive    &    F_L(n,r) is stored in R04

Example:

            3   ENTER^
         -0.4  ENTER^
          1.2  XEQ "ICWF"  >>>>  G₃( -0.4 , 1.2 ) = 6.563265348                                         ---Execution time = 66s---
                               and we have  F₃( -0.4 , 1.2 ) = 0.029547268          in register R04.

-Like with "ICWF", several significant digits may be cancelled when the last operation is performed.
 ( Check the value in LASTX or in R06 )
 

       c)  Non-integer Values of L
 

Formula:
 

   G_L(n,r) = i F_L(n,r) + | Gam(L+1+i.n) |  [ 1/Gam(L+1+i.n) ]  exp [ (PI) n/2 + i L(PI)/2 ]  W_i.n,L+1/2 (2i.r)

   where  W_k,µ = Whittaker's function of the 2nd kind
 

Data Registers:   R00 thru R15: temp
                              R06 = r  ,  R07 = n  ,  R08 = L
Flags: /
Subroutines:  "GAM3"  or another "gamma-function program"  ( cf "Gamma Function for the HP-41" )
                         "GAMZ+" or "GAMZ"    ( gamma function in the complex plane )
                         "RCWF3" cf § 1-c)

-The M-Code routines  Z*Z & Z/Z  may be replaced by  XEQ "Z*Z" & XEQ "Z/Z"  ( cf "Complex Functions for the HP-41" )
 
 

 
   ( 191 bytes / SIZE 016 )
 
 

 
Example:

   1.4   ENTER^
  -2.1   ENTER^
   1.6   XEQ "ICWF3"  >>>>   G_1.4( -2.1 , 1.6 ) = -0.206368268                         ---Execution time = 107s---

-The first term:  i F_L(n,r)  is not taken into account since this complex number is purely imaginary.
-We can, however, add this term to the final result to get an idea of the accuracy:
-Add  X<>Y  RCL 16  +  X<>Y  after line 141  and add  STO 16  after line 02
-Here, it yields  Y-output = 2 E-9 suggesting 8 or 9 correct decimals ( it should be 0 if there were no roundoff errors ).
-The HP-48 returns  -0.206368268  so the result is exact to 9D.

-This program does not work if L is an integer because we cannot use the same formulae to compute Whittaker's function of the 2nd kind.
 

3°)  Asymptotic Expansion
 

-This routine is useful to compute the Coulomb Wave Functions for large value of r.
-It also works if L is not an integer.
 

Formulae:

    F_L(n,r)  =  g cos µ + f sin µ       where       µ     =   r - n Ln (2r) - L.(Pi)/2 + Arg Gam ( 1 + L + i n )
    G_L(n,r)  =  f cos µ - g sin µ         and      f + i.g ~  u₀ + u₁ + u₂ + ............ + u_k + ......................

    with  u₀ = 1  and   u_k / u_k-1 = ( i n - L + k - 1 ) ( i n + L + k ) / ( 2k i r )        ( k > 0 )
 

Data Registers:  R00 thru R08: temp
Flags: /
Subroutines:  "GAMZ+"  ( Gamma Function in the Complex Plane - cf "Gamma Function for the HP-41" )

-The M-Code routine  Z*Z  may be replaced by  XEQ "Z*Z"  ( cf "Complex Functions for the HP-41" )
 
 

 
   ( 110 bytes / SIZE 009 )
 
 

 
Example:      L = 3 , n = 4 , r = 16

        FIX 8

   3   ENTER^
   4   ENTER^
  16  XEQ "AECWF"  >>>>    F₃(4,16) = -0.997776721                              ---Execution time = 38s---
                                   RDN     G₃(4,16) = -0.694942604

Notes:

-In fact, these series are divergent and line 72 ( RND ) is used to stop the loop when 2 rounded consecutive sums are equal
-Check registers R07 & R08 to see the last u_k =  -2 E-11 - 4 E-9 i  in the example above.

-If r is even larger, the routine works in FIX 9 or even in SCI 9
-For instance, "AECWF" returns  F₀(1,20) = -0.3292255394  and  G₀(1,20) = -0.9724283978

-This program can be modified so that the loop stops when 2 successive u_k start to increase:

  Replace line 83-84 by  RCL 09  X<>Y  RDN  RDN  +
  Replace line 72-73-74 by  X=0?  GTO 00  ENTER^  X<> 09  X>Y?  GTO 01  LBL 00
  and add  STO 09  after line 14

-Register Z and R09 will give an idea of the accuracy.

-It works often ... but not always! A premature stop is possible.
-But you can press XEQ 01 to check if the loop is still going on.
 

Reference:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4