%%HP: T(3)A(D)F(.);
@ QUANTUM by Paul Edgar Raines
DIR

W3J
@ This program calculates Wigner 3-j symbols.
@       It uses the following formula.
@
@       ³ J1 J2 J3 ³         J1-J2-M3
@       ³ M1 M2 M3 ³  M  (-1)         delta(J1,J2,J3) w3jf(J1,J2,J3,M1,M2,M3)
@
@                           ³                                       ³1/2
@   where                   ³ (J1+J2-J3)! (J1-J2+J3)! (-J1+J2+J3)!  ³
@         delta(J1,J2,J3) = ³ ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ ³
@                           ³            (J1+J2+J3+1)!              ³
@                           ³                                       ³
@                ³                                                       ³1/2
@   and   w3jf = ³ (J1+M1)! (J1-M1)! (J1+M1)! (J1-M1)! (J1+M1)! (J1-M1)! ³
@                ³                                                       ³
@      ÄÄÄÄÄ                               V
@       \                              (-1)
@    *   >  ÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ
@       /   V! (J1+J2-J3-V)! (J1-M1-V)! (J2+M2-V)! (J3-J2+M1+V)! (J3-J1-M2+V)!
@      ÄÄÄÄÄ
@        V
@
@   where V is over all values such that there are no negative factorials.
@
@   Enter Wigner 3-j symbols as a 2x3 matrix as shown below.
@                1: [[ 1.5 1 .5 ]
@                    [ .5 -1 .5 ]]
@   which produces
@                1:  .288675134595   or '{(1/12)'
@
\<< OBJ\-> DROP                      @ decompose w3j symbol
  \-> J1 J2 J3 M1 M2 M3              @ read in w3j elements
  \<< -1 J1 J2 M3 + - ^              @ calculate over all sign
      J1 J2 + J3 - !                 @ calculate delta part
      J1 J2 - J3 + ! *
      J2 J1 - J3 + ! *
      J1 J2 + J3 + 1 + ! /
      \v/ *

      J1 M1 + !                      @ calculate first part of w3jf
      J1 M1 - ! *
      J2 M2 + ! *
      J2 M2 - ! *
      J3 M3 + ! *
      J3 M3 - ! * \v/ *

      J1 J2 + J3 -                   @ calculate elements for the
      J1 M1 - J2 M2 +                @    second part of w3jf
      J3 J2 - M1 +                   @    to determine iteration
      J3 J1 - M2 -                   @    bounds for 'V'
      \-> W1 W2 W3 W4 W5             @ read in elements
      \<< W1 W2 MIN W3 MIN           @ determine lower bound for 'V'
          0 -1 W4 * MAX              @ determine upper bound for 'V'
          -1 W5 * MAX
          \-> XH XL                  @ read in bounds
          \<<
          IF 'XH<XL'                 @ if bounds improper
          THEN 1                     @      return 1
          ELSE 0 XL XH               @ else calculate second part
            FOR V V ! W1 V - ! *     @      of w3jf
              W2 V - ! * W3 V - ! *
              W4 V + ! * W5 V + ! *
              INV -1 V ^ * +
            NEXT
          END
          \>>
      \>>
      *                              @ multiply w3jf 2nd part to rest
   \>>
\>>

C&G
@ This program calculates Clebsch-Gordan coefficients.
@       Put the C-G coeff. on level one of the stack in
@       < J1 J2 M1 M2 3 J3 M3> form as a list.
@
@                        1:  { J1 J2 M1 M2 J3 M3 }
@       which gives
@                        1:
@
@       If you prefer < J1 M1 J2 M2 3 J3 M3> form or any other,
@       just switch the order on line 2 of the program to your preference.
@
\<< OBJ\-> DROP NEG                  @ decompose Clebsch-Gordan coefficient
  \-> J1 J2 M1 M2 J3 M3              @ read in C-G elements
  \<< -1 J1 J2 M3 + - ^              @ calculate over all sign
      2 J3 * 1 + \v/ *               @ calculate conversion factor
      J1 J2 J3 M1 M2 M3              @ calculate Wigner 3-j symbol
      { 2 3 } \->ARRY W3J *
  \>>
\>>

JGEN
@ J-matrices generator
@        This program generates the J^2, Jz, J-, and J+ matrices
@        for an arbitrary spin/ang.mom/j given on stack.
@        For a spin of 3/2 but 1.5 in level one.
\<< \-> S                       @ read spin/ang.mom./j into 'S'
  \<< S -1 * S FOR R            @ iterate 'R' from '-S' to 'S-1' {row}
      S -1 * S FOR C            @ iterate 'C' from '-S' to 'S-1' {col}
      R C - 1
      IF SAME THEN              @ if on bottom off diagonal calc
      S C - S C + 1 + * \v/     @    proper element for matrix
      ELSE 0 END                @ else put a zero
    NEXT NEXT
    S 2 * 1 +                   @ put together Jminus
    DUP 2 \->LIST \->ARRY
    DUP 'JM' STO                @ duplicate matrix and store as 'JM'
    TRN 'JP' STO                @ transpose Jminus and store as 'JP'
    0 S 2 * FOR R               @ iterate 'R' from '0' to '2*S' {row}
    0 S 2 * FOR C               @ iterate 'C' from '0' to '2*S' {col}
      R C SAME S R - *          @ if on diagonal, calc proper element
    NEXT NEXT
    S 2 * 1 +                   @ put together Jz
    DUP 2 \->LIST \->ARRY
    'Jz' STO                    @ store as 'Jz'
    S 2 * 1 + IDN S S 1 +       @ build Jsquared matrix
    * * 'J2' STO                @ store as 'J2'
  \>>
\>>

END