I’ve been trying to reincarnate an old FORTRAN Atomic Physics code from many years ago that I had kept around but had misplaced one subroutine called **ANGL**. It calculated Wigner 3J coefficients and I had spent some time figuring out how to do that in the dark ages.

With the modern Internet, I thought it would be a simple search and I’d find someone who had created a standard version of this. I’d then download it, adapt its argument list to match the requirements in my code and be done. Silly me!

It appears that while the definitions are available and there are Wolfram calculators and javascript calculators on line, there isn’t a FORTRAN version that I could find that would compile and work. I found a couple in C++, but they didn’t compile and work on my MacBook Pro, either.

It seems like it is some sort of dark secret how to actually calculate a given 3J symbol for use in a FORTRAN code. When looking back through the literature on line, each author seems to delight is switching notations and giving you definitions in terms of other definitions. Wikipedia isn’t much help here!

So, repeating the quest I did so many years ago, I worked out a version that works for my purposes and might be helpful as a starting point for some other poor slob in a similar boat. I’ll include the definition by Racah of the Wigner 3J symbols, and then show how they are implemented in the FORTRAN subroutine **ANGL**. Of course, it should be a simple task to take what I’m presenting here and port it into C, C++, Objective-C, or Java.

We start with the Racah definition of the Wigner 3J symbol.

In our calculations this is done by calculating tmin and tax as

` t1 = j2 - j3 - m1`

t2 = j1 - j3 + m2

t3 = j1 + j2 - j3

t4 = j1 - m1

t5 = j2 + m2

tmin = min( 0, min( t1, min(t2, min(t3, min(t4, t5)))))

tmax = max( t1, max(t2, max(t3, max(t4, t5))))

The codes below are written in FORTRAN 66 maybe or FORTRAN 77. That’s what I was going back to to start. I’m sure someone can convert these to more modern FORTRAN, but I liked the simplicity and wanted to show every step, not be clever!

The file ANGL.f contains subroutine **ANGL**. It does not check anything but conditions about all arguments being integer or half-integer. The accompaning subroutine **QANG**, which would be in the file QANG.f, does the checking and is the actual call to **ANGL**, determining those which are automatically zero and also has the framework for storing and retrieving previously calculated symbols. This will be updated in a later post. In addition, I’ve included the file Wigner3JMain.f to show how **QANG** is called, along with a shell script, run3J.sh which compiles and runs the software.

File ANGL.f

REAL*8 FUNCTION ANGL(j1, j2, j3, m1, m2, m3) ANGL0010
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 wigner, delta
COMMON numFactorials, ifact(10)
INTEGER t1, t2, t3, t4, t5, tmin, tmax
9999 FORMAT(I5)
ANGL = -999.9
if (j1+j2+j3+1 > 10 ) GOTO 100
ANGL = 0.0
if( 2*j1 .NE. floor(2.0*j1) ) GOTO 100
if( 2*j2 .NE. floor(2.0*j2) ) GOTO 100
if( 2*j3 .NE. floor(2.0*j3) ) GOTO 100
if( 2*m1 .NE. floor(2.0*m1) ) GOTO 100
if( 2*m2 .NE. floor(2.0*m2) ) GOTO 100
if( 2*m3 .NE. floor(2.0*m3) ) GOTO 100
10 CONTINUE
t1 = j2 - j3 - m1
t2 = j1 - j3 + m2
t3 = j1 + j2 - j3
t4 = j1 - m1
t5 = j2 + m2
tmin = min( 0, min( t1, min(t2, min(t3, min(t4, t5)))))
tmax = max( t1, max(t2, max(t3, max(t4, t5))))
wigner = 0.0
C Calculate Sum over I of (-1)^I/x
C Where x is I!(j3-j2+I+m1)!(j3-m1+I-m2)!(j1+j2=m3-I)!(j1-I-m1)!(j2-I+m2)!
DO 9 I = tmin, tmax
IX = ifactr(I)
IX = IX * ifactr(j3-j2+I+m1)
IX = IX * ifactr(j3-j1+I-m2)
IX = IX * ifactr(j1+j2-j3-I)
IX = IX * ifactr(j1-I-m1)
IX = IX * ifactr(j2-I+m2)
C CALCULATE (-1)**I
IY = (-1)**I
wigner = wigner + float(IY)/float(IX)
9 CONTINUE
C Calculate the triangle_coeff idelta
delta = float(ifactr(j1+j2-j3)*ifactr(j1-j2+j3)*ifactr(-j1+j2+j3))
delta = delta/float(ifactr(j1+j2+j3+1))
C Calculate sign
I1 = (-1)**(j1-j2-m3)
C Calculate the factorials in the square root with delta
IY = ifactr(j1+m1)
IY = IY * ifactr(j1-m1)
IY = IY * ifactr(j2+m2)
IY = IY * ifactr(j2-m2)
IY = IY * ifactr(j3+m3)
IY = IY * ifactr(j3-m3)
ANGL = I1 * wigner*sqrt(FLOAT(IY)*delta)
1000 FORMAT(1H1,46HAll arguments must be integer or half-integer.6(I5))
100 RETURN
END
INTEGER FUNCTION ifactr(K)
INTEGER numFactorials
COMMON numFactorials, ifact(10)
ifactr = 1;
if( K == 0 ) GOTO 10
ifactr = ifact(K)
10 CONTINUE
RETURN
END
INTEGER FUNCTION ifactorial( K )
J = 1;
IF ( K .EQ. 0 ) GO TO 11
DO 10 I=1,K
J = J*I
10 CONTINUE
11 ifactorial = J
return
END

File QANG.f

FUNCTION QANG(J1,J2,J3,M1,M2,M3) QANG0020
IMPLICIT REAL*8(A-H,O-Z) QANG0040
REAL*8 ANGL
INTEGER*2 JTABLE(700,3) QANG0060
DIMENSION TABLE(74) QANG0080
DATA NUM/503/ QANG0100
DATA IFLAG/1/ QANG0120
DATA ITEST/1/
IF( ITEST .NE. 1 ) GO TO 11
DO 9 I = 1, 700
DO 9 J = 1,3
JTABLE(I,J) = -1
9 CONTINUE
11 CONTINUE
IF(IFLAG.NE.0)GO TO 989 QANG0480
C... INITIALATION AND SORTED TABLE AND ITABLE GO HERE QANG0500
READ(10)((JTABLE(I,J),J=1,3),I=1,700) QANG0520
READ(10)(TABLE(I),I=1,74) QANG0540
989 IFLAG=1 QANG0560
QANG = 0.D0
IF( IABS(M1) .GT. J1 ) GOTO 973
IF( IABS(M2) .GT. J2 ) GOTO 973
IF( IABS(M3) .GT. J3 ) GOTO 973
IF( IABS(J1-J2) .GT. J3 ) GOTO 973
IF( J3 .GT. J1+J2 ) GOTO 973
IF( M1+M2+M3 .NE. 0) GOTO 973
GOTO 1073
INDEX=J1*100000+J2*10000+J3*1000+M1*100+M2*10+M3 QANG0700
IQUOT=INDEX/NUM QANG0720
IPROB=INDEX-IQUOT*NUM+1 QANG0740
C WRITE(6,9999) INDEX
C WRITE(6,9999) IPROB
IF(JTABLE(IPROB,1).EQ.-1)GO TO 999 QANG0760
100 CONTINUE QANG0780
IF(JTABLE(IPROB,1).NE.IQUOT)GO TO 200 QANG0800
QANG=TABLE(JTABLE(IPROB,2)) QANG0820
GO TO 973 QANG0860
200 IF(JTABLE(IPROB,3).EQ.-1)GO TO 999 QANG0880
IPROB=JTABLE(IPROB,3) QANG0900
GO TO 100 QANG0920
999 QANG=0.D0 QANG0940
1073 CONTINUE QANG1000
QANG=ANGL(J1,J2,J3,M1,M2,M3) QANG1020
C TABLE(IPROB,1) = IQUOT
C TABLE(IPROB,2) = QANG
973 CONTINUE QANG1040
9999 FORMAT(1H1,I10)
RETURN QANG1060
END QANG1080

File Wigner3JMain.f

C PROGRAM WIGNER
C..... MAIN PROGRAM FOR TESTING AND RUNNING QANG
REAL*8 X, QANG, ANGL
INTEGER numFactorials
COMMON numFactorials, ifact(10)
C
numFactorials = 10
C
DO 8, I=1,numFactorials
ifact(I) = ifactorial(I)
8 CONTINUE
WRITE(6,9999) ifactr(0), ifactr(1), ifactr(2),ifactr(10)
C GOTO 101
9999 FORMAT(3I5I30)
I=2
J=2
K=0
II=0
JJ=0
KK=0
X = ANGL(I,J,K,II,JJ,KK)
C FOR (2,2,0 / 0,0,0 ) = sqrt(1/5)
WRITE(6,1001)I,J,K,II,JJ,KK,X
I=1
J=1
K=0
X = ANGL(I,J,K,II,JJ,KK)
WRITE(6,1001)I,J,K,II,JJ,KK,X
X = QANG(I,J,K,II,JJ,KK)
WRITE(6,1001)I,J,K,II,JJ,KK,X
C GO TO 101
Write(6,1000)
1000 FORMAT(1H1,"Starting to calculate some 3J symbols\n")
DO 100 I = 0,2
DO 100 J = 0,2
DO 100 K = 0,2
DO 100 II = 0,2
DO 100 JJ = 0,2
DO 100 KK = 0,2
X = QANG(I,J,K,II,JJ,KK)
IF( X .EQ. 0.0 ) GOTO 100
WRITE(6,1001)I,J,K,II,JJ,KK,X
1001 FORMAT(1H1,I3,1X,I3,1X,I3,1X,I3,1X,I3,1X,I3,1X,F10.3)
100 CONTINUE
101 CONTINUE
WRITE(6,300)
300 FORMAT(1H1,"ALL DONE!\n")
STOP
END

The shell script to run this on Mac OS X 10.10.3 using the gcc-4.9.bin after installing XCode 6.3 and the command line tools ( Mac OS X has made it difficult to use the underlying Unix functionality by moving things to unconventional places. Unless you install XCode and the command line tools, you are missing most of the Unix command line features. )

File run3j.sh

gfortran -v
gfortran -o Wigner -lcrt1.o -ffpe-trap='underflow' Wigner3JMain.f QANG.f ANGL.f
Wigner

The output for running the above shell, if everything works right and all the fortran guts are installed correctly, etc. is:

19:48:40 [woo:~/Development … e/fortran/3J] > run3J.sh
Using built-in specs.
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=/usr/local/libexec/gcc/x86_64-apple-darwin13.0.0/4.9.0/lto-wrapper
Target: x86_64-apple-darwin13.0.0
Configured with: ../gcc-4.9-20131215/configure --enable-languages=c++,fortran
Thread model: posix
gcc version 4.9.0 20131215 (experimental) (GCC)
1 1 2 3628800
1 2 2 0 0 0 0 0.447
1 1 1 0 0 0 0 -0.577
1 1 1 0 0 0 0 -0.577
1Starting to calculate some 3J symbols\n
1 0 0 0 0 0 0 1.000
1 0 1 1 0 0 0 -0.577
1 0 2 2 0 0 0 0.447
1 1 0 1 0 0 0 -0.577
1 1 1 0 0 0 0 -0.577
1 1 1 2 0 0 0 0.365
1 1 2 1 0 0 0 0.365
1 2 0 2 0 0 0 0.447
1 2 1 1 0 0 0 0.365
1 2 2 0 0 0 0 0.447

I’ll put these files into a download directory soon and arrange a link to them. Good luck!

### Like this:

Like Loading...