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  X<Y?
  24  X<>Y
  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  X<Y?
  53  GTO 02
  54  RDN
  55  FACT
  56  RCL 11
  57  RCL 07
  58  -
  59  FACT
  60  *
  61  RCL 12
  62  RCL 07
  63  -
  64  FACT
  65  *
  66  RCL 13
  67  RCL 07
  68  -
  69  FACT
  70  *
  71  RCL 04
  72  RCL 07         
  73  +
  74  RCL 00
  75  -
  76  FACT
  77  *
  78  RCL 07
  79  RCL 05
  80  -
  81  RCL 09
  82  -
  83  FACT
  84  *
  85  1
  86  ST+ 07
  87  CHS
  88  RCL 07
  89  Y^X
  90  X<>Y
  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
0F8   READ 3(X)
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  X<Y?
  24  X<>Y
  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  X<Y?
  52  GTO 02
  53  RDN
  54  FACT
  55  RCL 09
  56  RCL 07
  57  -
  58  FACT
  59  *
  60  RCL 10
  61  RCL 07
  62  -
  63  FACT
  64  *
  65  RCL 11         
  66  RCL 07
  67  -
  68  FACT
  69  *
  70  RCL 07
  71  RCL 00
  72  -
  73  FACT
  74  *
  75  RCL 07
  76  RCL 08
  77  -
  78  FACT
  79  *
  80  1
  81  ST+ 07
  82  CHS
  83  RCL 07
  84  Y^X
  85  X<>Y
  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  X<Y?
  22  X<>Y
  23  RCL 01
  24  LASTX
  25  RCL 06
  26  ST+ 12
  27  +
  28  ST+ 13
  29  +
  30  STO 09         
  31  X<Y?
  32  X<>Y
  33  RCL 02
  34  RCL 04
  35  +
  36  RCL 06
  37  +
  38  STO 10
  39  X<Y?
  40  X<>Y
  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  X<Y?
  55  GTO 02
  56  CLX
  57  RCL 00
  58  -
  59  FACT
  60  RCL 07
  61  RCL 08
  62  -
  63  FACT
  64  *
  65  RCL 07
  66  RCL 09
  67  -
  68  FACT
  69  *
  70  RCL 07
  71  RCL 10
  72  -
  73  FACT
  74  *
  75  RCL 11
  76  RCL 07
  77  -
  78  FACT
  79  *
  80  RCL 12
  81  RCL 07         
  82  -
  83  FACT
  84  *
  85  RCL 13
  86  RCL 07
  87  -
  88  FACT
  89  *
  90  1
  91  ST+ 07
  92  CHS
  93  RCL 07
  94  Y^X
  95  LASTX
  96  FACT
  97  *
  98  X<>Y
  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  X<Y?
14  X<>Y
15  RCL 02
16  RCL 06
17  -
18  ABS
19  X<Y?
20  X<>Y
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<Y?
42  GTO 02
43  RCL 15
44  STO 01      
45  RCL 23
46  STO 02
47  RCL 22
48  STO 04
49  RCL 18
50  STO 05
51  RCL 21
52  STO 06
53  XEQ "W6J"
54  STO 26
55  RCL 04
56  STO 01
57  RCL 05
58  STO 02
59  RCL 20
60  STO 04
61  RCL 16
62  STO 05
63  RCL 19
64  STO 06
65  XEQ "W6J"
66  ST* 26
67  RCL 04
68  STO 01
69  RCL 05
70  STO 02
71  RCL 15
72  STO 04
73  RCL 23
74  STO 05
75  RCL 17
76  STO 06
77  XEQ "W6J"
78  RCL 26
79  *
80  1
81  CHS
82  RCL 03
83  ST+ X
84  Y^X
85  LASTX
86  1
87  ST+ 03
88  +
89  *
90  *
91  ST+ 24
92  GTO 01
93  LBL 02
94  RCL 27      
95  REGSWAP
96  RCL 24
97  END

 
    ( 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/