# hp41programs

Clebsch-Gordan & Wigner 3j- 6j- 9j-

# Clebsch-Gordan Coefficients & Wigner 3j-6j-9j-symbols for the HP-41

Overview

1°)  Clebsch-Gordan Coefficients & Wigner 3j-symbol  ( Program#1 )
2°)  Triangular Coefficients
3°)  Clebsch-Gordan Coefficients & Wigner 3j-symbol  ( Program#2 )
4°)  Wigner 6j-symbol
5°)  Wigner 9j-symbol

-These numbers are useful to describe the angular momentum in quantum physics.
-All the parameters ji , mi must be integers or half-integers.

-The standard notations for Wigner symbols use parentheses and brace brackets beside 2x3 or 3x3 arrays. This is not easy in html.
-So, I've used non-standard notations like W6j ( j1 j2 j3 j4 j5 j6 ) for the Wigner 6j- symbols...

1°)  Clebsch-Gordan Coefficients & Wigner 3j-symbol  ( Program#1 )

-This program returns  <  j1 j2  m1 m2 |   j1 j2 j m  >  in X-register
and      W3J (  j1 j2 j  m1 m2 -m )    in Y-register

Formulae:

•  Clebsch-Gordan Coefficients  = <  j1 j2  m1 m2 |   j1 j2 j m  >
= [ ( 2j +1) D( j1 j2 j  ) ]1/2  [ ( j1 + m1 )! ( j1 - m1 )! ( j2 + m2 )! ( j2 - m2 )! ( j + m )! ( j - m )! ]
x Sumk  (-1)k / [ k! ( j1 + j2 - j - k )! ( j1 - m1 - k )! ( j2 + m2 - k )! ( j - j2 + m1 + k )! ( j - j1 - m2 + k )! ]         if   m1 + m2 = m

( The sum is to be performed for all integers k that produce non-negative integer arguments for the factorials ).

•  <  j1 j2  m1 m2 |   j1 j2 j m  > = 0            if   m1 + m2 # m

•  Wigner 3j-symbol  = W3J (  j1 j2 j  m1 m2 -m ) = (-1)^( j1 - j2 + m )  ( 2j +1) -1/2  <  j1 j2  m1 m2 |   j1 j2 j m  >

•   Triangular coefficients  D( a b c )  are defined in paragraph 2

Data Registers:              R00 & R07 thru R13: temp                      ( Registers R01 thru R06 are to be initialized before executing "CGW3J" )

•  R01 = j1     •  R02 = j2    •  R03 = j
•  R04 = m1   •  R05 = m2   •  R06 = m
Flags: /
Subroutines: /

-If  m1 + m2 # m , the results = 0
-Otherwise, they are less easy to compute...

 01  LBL "CGW3J"   02  RCL 04   03  RCL 05   04  +   05  RCL 06   06  X=Y?   07  GTO 00   08  CLST   09  RTN   10  LBL 00   11  RCL 02   12  RCL 03   13  -   14  STO 00   15  RCL 04   16  -   17  RCL 01   18  RCL 03   19  -   20  STO 09   21  RCL 05   22  +   23  XY   25  X<0?   26  CLX   27  STO 07   28  RCL 01   29  RCL 02   30  +   31  STO 10   32  RCL 03   33  -   34  STO 11 35  RCL 01   36  RCL 04   37  -   38  STO 12            39  X>Y?   40  X<>Y   41  RCL 02   42  RCL 05   43  +   44  STO 13   45  X>Y?   46  X<>Y   47  STO 08   48  CLX   49  LBL 01   50  RCL 07   51  RCL 08   52  XY   91  /   92  -   93  GTO 01   94  LBL 02   95  X<> Z   96  RCL 11   97  FACT   98  SQRT   99  * 100  RCL 01 101  RCL 00 102  - 103  FACT 104  SQRT 105  * 106  RCL 02          107  RCL 09 108  - 109  FACT 110  SQRT 111  * 112  RCL 01 113  RCL 04 114  + 115  FACT 116  SQRT 117  * 118  RCL 12 119  FACT 120  SQRT 121  * 122  RCL 13 123  FACT 124  SQRT 125  * 126  RCL 02 127  RCL 05 128  - 129  FACT 130  SQRT 131  * 132  RCL 03 133  RCL 06 134  + 135  FACT 136  SQRT 137  * 138  RCL 03 139  RCL 06          140  - 141  FACT 142  SQRT 143  * 144  RCL 03 145  RCL 10 146  + 147  1 148  + 149  FACT 150  SQRT 151  / 152  1 153  CHS 154  RCL 01 155  RCL 02 156  - 157  RCL 06 158  + 159  Y^X 160  * 161  X<>Y 162  RCL 03 163  ST+ X 164  1 165  + 166  SQRT 167  * 168  END

( 184 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS Y / W3j X / ClG

where    W3j = W3j (  j1 j2 j  m1 m2 -m )     and   ClG = <  j1 j2  m1 m2 |   j1 j2 j m  >

Examples:

-If   j1 = 12  ,   j2 = 24  ,  j  = 31
m1 = 1   ,  m2 = 16  , m = 17

XEQ "CLGW3J"  returns   ClG =  0.206754081                 ---Execution time = 15s---
X<>Y   W3j = -0.026048566

-With   j1 = 1  ,   j2 = 3/2  ,  j  = 5/2
m1 = 0 ,  m2 = 3/2  , m = 3/2

XEQ "CLGW3J"  yields   ClG =  0.632455533                 ---Execution time = 5s---
X<>Y   W3j = -0.258198890

Notes:

-All the j's & m's must be integers or half-integers
-They must satisfy several inequalities so that  FACT is applied to non-negative integers only.
-Otherwise, DATA ERROR will be displayed.
-If you want to get  W3j (  j1 j2 j  m1 m2 m )  store (-m) into R06

-You can save a few bytes if you perform SQRT after the different products and quotient, but there will be more risks of overflows.
-Furthermore, the HP-41 actually computes n! up to n = 99 internally, so we can write a M-Code function to calculate sqrt(n!) and avoid an "OUT OF RANGE"

-Replace for instance  FACT  SQRT   by the M-Code function  SQRTF  listed below ( it does not check for alpha data ):

086   "F"
014   "T"
012   "R"
011   "Q"
013   "S"
2A0  SETDEC
006   A=0 S&X
166   A=A+1 S&X
306   ?A<C S&X
0B5   ?C GO
0A3   282D            GTO DATA ERROR if x is not an integer or if x > 99
128   WRIT 4(L)
3B1   ?NCXQ         C=
060   18EC            fact(C)
2F9   ?NCQ            C=
060   18BE            sqrt(C)
0E8   WRIT 3(X)
3E0   RTN

-After these modifications,

with   j1 = 14  ,   j2 = 29  ,  j  = 41
m1 = 1   ,  m2 = 16  , m = 17

XEQ "CLGW3J"   gives    ClG =  0.361681799                 ---Execution time = 10s---
X<>Y   W3j =  0.039699735

( sqrt ( 85! ) has been computed )

2°)  Triangular coefficients

-Wigner 3j- 6j- and 9j- symbols use the calculation of

Delta ( a , b , c ) = [ ( a + b - c )! ( a - b + c )! ( -a + b + c )! ] / ( a + b + c + 1 )!

-The routine hereunder computes  SQRT [ Delta ( a , b , c ) ]
-FACT  SQRT  may be replaced by the M-Code routine  SQRTF listed above...
a , b , c  are stored into  R11 , R12 , R13  respectively
we can also use synthetic registers M , N , O

 01  LBL "DABC" 02  STO 13         03  X<>Y 04  STO 12 05  - 06  X<>Y 07  STO 11 08  + 09  FACT 10  SQRT 11  RCL 11 12  RCL 12         13  + 14  RCL 13 15  - 16  FACT 17  SQRT 18  * 19  RCL 12         20  RCL 11 21  - 22  RCL 13 23  + 24  FACT 25  SQRT 26  * 27  RCL 11 28  RCL 12         29  + 30  RCL 13 31  + 32  1 33  + 34  FACT 35  SQRT           36  / 37  END

( 46 bytes )

 STACK INPUTS OUTPUTS T T T Z a T Y b T X c sqrt(delta(a,b,c))

Example:

12  ENTER^
24  ENTER^
31  XEQ "DABC"  returns  5.963240385 E-13

3°)  Clebsch-Gordan Coefficients & Wigner 3j-symbol  ( Program#2 )

-We can now rewrite "CGW3J" to use "DABC"
-This will save bytes if you plan to use also "W6J" and/or "W9J".
-Otherwise, it's wasteful!

Data Registers:              R00 & R07 thru R13: temp                      ( Registers R01 thru R06 are to be initialized before executing "CGW3J" )

•  R01 = j1     •  R02 = j2    •  R03 = j
•  R04 = m1   •  R05 = m2   •  R06 = m
Flags: /
Subroutine:  "DABC"

 01  LBL "CGW3J"   02  RCL 04   03  RCL 05   04  +   05  RCL 06   06  X=Y?   07  GTO 00   08  CLST   09  RTN   10  LBL 00   11  RCL 02   12  RCL 03   13  -   14  RCL 04   15  -   16  STO 00   17  RCL 01   18  RCL 03   19  -   20  RCL 05   21  +   22  STO 08   23  XY   25  X<0?   26  CLX   27  STO 07   28  RCL 01   29  RCL 02 30  +   31  RCL 03   32  -   33  STO 09   34  RCL 01   35  RCL 04            36  -   37  STO 10   38  X>Y?   39  X<>Y   40  RCL 02   41  RCL 05   42  +   43  STO 11   44  X>Y?   45  X<>Y   46  STO 12   47  CLX   48  LBL 01   49  RCL 07   50  RCL 12   51  XY   86  /   87  - 88  GTO 01   89  LBL 02   90  X<> Z   91  RCL 01   92  RCL 04   93  +   94  FACT   95  SQRT   96  *   97  RCL 10            98  FACT   99  SQRT 100  * 101  RCL 11 102  FACT 103  SQRT 104  * 105  RCL 02 106  RCL 05 107  - 108  FACT 109  SQRT 110  * 111  RCL 03 112  RCL 06 113  + 114  FACT 115  SQRT 116  * 117  RCL 03 118  RCL 06 119  - 120  FACT 121  SQRT 122  * 123  RCL 01 124  RCL 02 125  RCL 03 126  XEQ "DABC" 127  * 128  1 129  CHS 130  RCL 01 131  RCL 02 132  - 133  RCL 06 134  + 135  Y^X 136  * 137  X<>Y 138  RCL 03 139  ST+ X 140  1 141  + 142  SQRT 143  * 144  END

( 165 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS Y / W3j X / ClG

where    W3j = W3j (  j1 j2 j  m1 m2 -m )     and   ClG = <  j1 j2  m1 m2 |   j1 j2 j m  >

Example:

-If   j1 = 12  ,   j2 = 24  ,  j  = 31
m1 = 1   ,  m2 = 16  , m = 17

XEQ "CLGW3J"  returns   ClG =  0.206754081                 ---Execution time = 15s---
X<>Y   W3j = -0.026048566

4°)  Wigner 6j-symbol

Formula:

W6J ( j1 j2 j3 j4 j5 j6 )  =  [ D ( j1 j2 j3 ) D ( j1 j5 j6 ) D ( j2 j3 j4 ) D ( j3 j4 j5 ) ]1/2

x Sumk=k' to k"  ( -1 )k ( k + 1 )! / [ ( k - j1 - j2 - j3 )! ( k - j1 - j5 - j6 )! ( k - j2 - j3 - j4 )! ( k - j3 - j4 - j5 )!
( j1 + j2 + j4 + j5 - k )! ( j1 + j3 + j4 + j6 - k )! ( j2 + j3 + j5 + j6 - k )! ]

with   k' = max { j1 + j2 + j3 , j1 + j5 + j6 , j2 + j3 + j4 , j3 + j4 + j5 }
and    k" = min { j1 + j2 + j4 + j5 , j1 + j3 + j4 + j6 , j2 + j3 + j5 + j6 }

Data Registers:              R00 & R07 thru R14: temp                      ( Registers R01 thru R06 are to be initialized before executing "W6J" )

•  R01 = j1     •  R02 = j2    •  R03 = j3
•  R04 = j4     •  R05 = j5    •  R06 = j6
Flags: /
Subroutine:  "DABC"

 01  LBL "W6J"   02  RCL 01   03  STO 12   04  RCL 02   05  STO 13   06  +   07  STO 11            08  RCL 03   09  ST+ 13   10  +   11  STO 00   12  LASTX   13  RCL 04   14  ST+ 11   15  +   16  ST+ 12   17  RCL 05   18  ST+ 11   19  +   20  STO 08   21  XY   23  RCL 01   24  LASTX   25  RCL 06 26  ST+ 12   27  +   28  ST+ 13   29  +   30  STO 09            31  XY   33  RCL 02   34  RCL 04   35  +   36  RCL 06   37  +   38  STO 10   39  XY   41  STO 07   42  RCL 11   43  RCL 12   44  X>Y?   45  X<>Y   46  RCL 13   47  X>Y?   48  X<>Y   49  STO 14   50  CLX 51  LBL 01   52  RCL 07   53  RCL 14            54  XY   99  / 100  - 101  GTO 01 102  LBL 02 103  X<> Z 104  RCL 01 105  RCL 02 106  RCL 03 107  XEQ "DABC" 108  * 109  RCL 03 110  RCL 04 111  RCL 05 112  XEQ "DABC" 113  * 114  RCL 01 115  RCL 05 116  RCL 06 117  XEQ "DABC" 118  * 119  RCL 02 120  RCL 04 121  RCL 06 122  XEQ "DABC" 123  * 124  END

( 162 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS X / W6j

with  W6j = W6j ( j1 j2 j3 j4 j5 j6 )

Example:

-If   j1 = 10  ,   j2 = 16  ,  j3  = 21    store them     R01   R02   R03
j4 = 24  ,   j5 = 12  ,  j6  = 14         into           R04   R05   R06

XEQ "W6J"  >>>>   W6j  =  0.006251585                 ---Execution time = 21s---

Notes:

-All the j's must be integers or half-integers
-They must satisfy several inequalities so that FACT is applied to non-negative integers only.
-Otherwise, DATA ERROR will be displayed.

5°)  Wigner 9j-symbol

Formula:

W9J ( j1 j2 j3 j4 j5 j6 j7 j8 j9 )  =  Sumk=k' to k"  (-1)2k (2k+1)  W6J ( j1 j9 k j8 j4 j7 )  W6J ( j8 j4 k j6 j2 j5 )  W6J ( j6 j2 k j1 j9 j3 )

where    k' = max { | j1 - j9 | , | j4 - j8 | , | j2 - j6 | }    and    k" =  min {  j1 + j9  ,  j4 + j8  ,  j2 + j6 }

Data Registers:              R00 & R10 thru R27: temp                      ( Registers R01 thru R09 are to be initialized before executing "W9J" )

•  R01 = j1     •  R02 = j2    •  R03 = j3
•  R04 = j4     •  R05 = j5    •  R06 = j6              registers R01 to R09 are saved in registers R15 to R23
•  R07 = j7     •  R08 = j8    •  R09 = j9              and they are restored at the end
Flags: /
Subroutine:  "W6J"

-If you don't have an X-Functions Module, lines 02-03-04  may be replaced by  1.009   15   XEQ "LCO"
-And lines 94-95  may be replaced by    15.023   1   XEQ "LCO"
( cf "Miscellaneous Short Routines for the HP-41" )

 01  LBL "W9J" 02  1.015009 03  STO 27 04  REGMOVE 05  RCL 01       06  RCL 09 07  - 08  ABS 09  RCL 04 10  RCL 08 11  - 12  ABS 13  XY 15  RCL 02 16  RCL 06 17  - 18  ABS 19  XY 21  STO 03 22  RCL 01 23  RCL 09 24  + 25  RCL 04       26  RCL 08 27  + 28  X>Y? 29  X<>Y 30  RCL 02 31  RCL 06 32  + 33  X>Y? 34  X<>Y 35  STO 25 36  CLX 37  STO 24 38  LBL 01 39  RCL 03 40  RCL 25 41  X

( 149 bytes / SIZE 028 )

 STACK INPUTS OUTPUTS X / W9j

where  W9j = W9J ( j1 j2 j3 j4 j5 j6 j7 j8 j9 )

Examples:

j1 = 3  ,   j2 = 7  ,  j3  = 5              R01   R02   R03
-If     j4 = 6  ,   j5 = 8  ,  j6  = 9      =      R04   R05   R06
j7 = 4  ,   j8 = 5  ,  j9  = 7              R07   R08   R09

XEQ "W9J"  >>>>   W9j  =  0.001032402067               ---Execution time = 4m32s---

j1 = 1/2  ,   j2 =  1    ,  j3  = 3/2
-With      j4 =  1    ,   j5 = 1/2  ,  j6  = 1/2
j7 = 3/2  ,   j8 = 1/2  ,  j9  =  1

XEQ "W9J"  >>>>   W9j  =  -0.02777777771   ( in fact -1/36 )         ---Execution time = 58s---

Notes:

-All the j's must be integers or half-integers
-They must satisfy several inequalities so that  FACT is applied to non-negative integers only.
-Otherwise, DATA ERROR will be displayed.

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   http://functions.wolfram.com/