From b262e89b717a90eb5130ef48e0acc47716fdd433 Mon Sep 17 00:00:00 2001
From: Joey Dumont <joey.dumont@gmail.com>
Date: Tue, 6 May 2014 10:42:43 -0400
Subject: [PATCH] First commit of working code. Needs cleanup.

---
 fortranLinkage.h  |  10 +
 testWigner.cpp    | 220 ++++++++++++
 wignerSymbols.cpp | 463 +++++++++++++++++++++++++
 wignerSymbols.f   | 867 ++++++++++++++++++++++++++++++++++++++++++++++
 wignerSymbols.h   | 166 +++++++++
 wignerTest.f90    |  11 +
 wigner_wrap.f90   |  51 +++
 7 files changed, 1788 insertions(+)
 create mode 100644 fortranLinkage.h
 create mode 100644 testWigner.cpp
 create mode 100644 wignerSymbols.cpp
 create mode 100644 wignerSymbols.f
 create mode 100644 wignerSymbols.h
 create mode 100644 wignerTest.f90
 create mode 100644 wigner_wrap.f90

diff --git a/fortranLinkage.h b/fortranLinkage.h
new file mode 100644
index 0000000..be4973e
--- /dev/null
+++ b/fortranLinkage.h
@@ -0,0 +1,10 @@
+#ifndef FORTRANLINKAGE_H
+#define FORTRANLINKAGE_H
+
+extern "C"
+{
+	extern void drc3jj_wrap(double,double,double,double,double*,double*,double*,int,int*);
+	extern void drc6j_wrap(double,double,double,double,double*,double*,double*,int,int*);
+}
+
+#endif // FORTRANLINKAGE_H
\ No newline at end of file
diff --git a/testWigner.cpp b/testWigner.cpp
new file mode 100644
index 0000000..f3e8d27
--- /dev/null
+++ b/testWigner.cpp
@@ -0,0 +1,220 @@
+#include "wignerSymbols.h"
+
+#include <iostream>
+#include <iomanip>
+#include <armadillo>
+#include <time.h>
+
+#define EPS 1.0e-10
+
+using namespace std;
+using namespace arma;
+
+// Orthogonality relations and sum rules for the 3j symbols [Brink & Satchler, pp. 136+139].
+double firstOrthoRelation3j(double l1, double l2, double l3, double m3);
+double sumOverM23j(double l1, double l2, double l3, double m3);
+
+// Orthogonality relations for the 6j symbols [Brink & Satchler pp. 142-3].
+double firstSumOverL3(double l1, double l2, double l6);
+double seconSumOverL3(double l1, double l2, double l6);
+
+// Timing computation of Wigner symbols. 
+double timingWignerSymbols(double l2, double l3, double m1, double m2, double m3);
+double timingWignerSymbolsF(double l2, double l3, double m1, double m2, double m3);
+
+int main(int argc, char* argv[])
+{
+	// We time the computation for different sizes.
+	int N = 500;
+	double l2(500.0),m1(-10.0),m2(60.0),m3(-50.0);
+	mat times(N,2);
+	mat timesF(N,2);
+
+	for (int i=0;i<N;i++)
+	{
+		int l3 = i+std::fabs(m3);
+		times(i,0) = (int)std::ceil(l2+l3-std::max(std::fabs(l2-l3),std::fabs(m1)))+1;
+		timesF(i,0) = times(i,0);
+		times(i,1) = timingWignerSymbols(l2,l3,m1,m2,m3);
+		timesF(i,1) = timingWignerSymbolsF(l2,l3,m1,m2,m3);
+	}
+
+	times.save("times.dat", raw_ascii);
+	timesF.save("timesF.dat",raw_ascii);
+
+	// Generate values for l1, l2 and derive the rest of the allowed values.
+	double lMax = atof(argv[1]);
+	colvec lVec = linspace(0,lMax,lMax+1);
+
+	// Run over l1
+	for (int i=0;i<lMax+1;i++)
+	{
+		double l1 = lVec(i);
+
+		// Run over l2
+		for (int j=0;j<lMax+1;j++)
+		{
+			
+			double l2 = lVec(j);
+			double l3min = std::fabs(l1-l2);
+			double l3max = l1+l2;
+
+			// Run over l3
+			for (int k=0;k<(l3max-l3min+1);k++)
+			{
+				double l3 = l3min+k;
+
+				double test = firstSumOverL3(l1,l2,l3); 
+
+				if (!test)
+					cout << "The first sum over L3 test failed for \n\tl1 = " << l1 << "\n\tl2 = " << l2 << "\n\tl6 = " << l3 << endl;
+
+				// Run over m3
+				mat precisionWigner(2*l3+1,2);
+				for (int l=0;l<(2*l3+1);l++)
+				{
+					double m3 = -l3+l;
+					test = firstOrthoRelation3j(l1,l2,l3,m3);
+
+					precisionWigner(l,0) = m3;
+					precisionWigner(l,1) = test;
+
+					if (!test)
+						cout << "The first orthogonality test failed for \n\tl1 = " << l1 << "\n\tl2 = " << l2 << "\n\tl3 = " << l3 << "\n\tm3 = " << m3 << endl;
+
+					//if (l3!=0.0) test = sumOverM23j(l1,l2,l3,m3);
+
+					if (!test)
+						cout << "The sum over M2 test failed for \n\tl1 = " << l1 << "\n\tl2 = " << l2 << "\n\tl3 = " << l3 << "\n\tm3 = " << m3 << endl;
+				}
+				precisionWigner.save("precision.dat", raw_ascii);
+			}
+
+			double l6max = std::min(2.*l1,2.*l2);
+			
+			// Run over l6
+			for (int k=0;k<l6max+1;k++)
+			{
+				bool test = seconSumOverL3(l1,l2,k);
+
+				if (!test)
+					cout << "The second sum over L3 test has failed for \n\tl1 = " << l1 << "\n\tl2 = " << l2 << "\n\tl6 = " << k << endl;
+			}
+		}
+	}
+
+	return 0;
+}
+
+// Sum_{m1,m2} (2*l3+1)*Wigner3j(l1,l2,l3,m1,m2,m3)*3j(l1,l2,l3',m1,m2,m3') = delta(l3,l3')*delta(m3,m3')
+double firstOrthoRelation3j(double l1, double l2, double l3, double m3)
+{
+	// We determine all possible values of m1 and m2.
+	int m1Size = ((int)(2.0*l1+1.0));
+	int m2Size = ((int)(2.0*l2+1.0)); 
+
+	double sum = 0.0;
+	for (int i=0;i<m1Size;i++)
+	{
+		for (int j=0;j<m2Size;j++)
+		{
+			sum += (2.0*l3+1.0)*varPhase::c_wigner3j(l1,l2,l3,((double)(-l1+i)),((double)(-l2+j)),m3)
+				*varPhase::c_wigner3j(l1,l2,l3,((double)(-l1+i)),((double)(-l2+j)),m3);
+		}
+	}
+
+	double diff = std::fabs(1.0-sum);
+
+	return diff;
+}
+
+// sum_{m1,m2} m2*(2*l3+1)c_wigner3j(l1,l2,l3,m1,m2,-m3)^2 = m3*(l3*(l3+1)+l2*(l2+1)-l1*(l1+1))/(2*l3*(l3+1))
+double sumOverM23j(double l1, double l2, double l3, double m3)
+{
+	double sum = 0.0;
+	// We sum over all allowable values of m1+m2+m3=0.
+	for (int i=0;i<(2*l1+1);i++)
+	{
+		double m1 = -l1+i;
+
+		for (int j=0;j<(2*l2+1);j++)
+		{
+			double m2 = -l2+j;
+
+			if (std::fabs(m1+m2-m3)<1.0e-6)
+				sum += m2*(2.0*l3+1.0)*pow(varPhase::c_wigner3j(l1,l2,l3,m1,m2,-m3),2.0);
+		}
+	}
+
+	double value = m3*(l3*(l3+1.)+l2*(l2+1.)-l1*(l1+1.))/(2.*l3*(l3+1.));
+
+	return std::fabs(sum-value);
+}
+
+// sum_{l3} = (-1)^(2*l3)*(2*l3+1)*Wigner6j(l1,l2,l3,l1,l2,l6) = 1
+double firstSumOverL3(double l1, double l2, double l6)
+{
+	// We sum over the value of l3. 
+	double l3min = std::fabs(l2-l1);
+	double l3max = l1+l2;
+
+	double sum = 0.0;
+	for (int i=0;i<l3max-l3min+1;i++)
+	{
+		double l3 = l3min+i;
+		sum += pow(-1.,2.0*l3)*(2.0*l3+1.0)*varPhase::c_wigner6j(l1,l2,l3,l1,l2,l6);
+	}
+
+	double diff = std::fabs(1.0-sum);
+
+	return diff;
+}
+
+// sum_{l3} = (-1)^(l1+l2+l3)*(2*l3+1)*Wigner6j(l1,l2,l3,l2,l1,l6) = delta(l6,0)*sqrt((2*l1+1)*(2*l2+1))
+double seconSumOverL3(double l1, double l2, double l6)
+{
+		// We sum over the value of l3. 
+	double l3min = std::fabs(l1-l2);
+	double l3max = l1+l2;
+
+	double sum = 0.0;
+	for (int i=0;i<l3max-l3min+1;i++)
+	{
+		double l3 = l3min+i;
+		sum += pow(-1.,l1+l2+l3)*(2.0*l3+1.0)*varPhase::c_wigner6j(l1,l2,l3,l2,l1,l6);
+	}
+
+	double value = (l6==0.0 ? sqrt((2.*l1+1.)*(2.*l2+1.)) : 0.0) ;
+
+	double diff = std::fabs(value-sum);
+
+	return std::fabs(value-sum);
+	
+}
+
+double timingWignerSymbols(double l2, double l3, double m1, double m2, double m3)
+{
+	// We start the timer. 
+	clock_t start = clock();
+
+	// We call the subroutine.
+	colvec test = varPhase::c_wigner3j(l2,l3,m1,m2,m3);
+
+	clock_t end = clock();
+
+	return ((double)end-start)/CLOCKS_PER_SEC;
+}
+
+double timingWignerSymbolsF(double l2, double l3, double m1, double m2, double m3)
+{
+	// We start the timer. 
+	clock_t start = clock();
+
+	// We call the subroutine.
+	double test = varPhase::wigner3j(l2+l3,l2,l3,m1,m2,m3);
+
+	clock_t end = clock();
+
+	return ((double)end-start)/CLOCKS_PER_SEC;
+}
+
diff --git a/wignerSymbols.cpp b/wignerSymbols.cpp
new file mode 100644
index 0000000..da23b45
--- /dev/null
+++ b/wignerSymbols.cpp
@@ -0,0 +1,463 @@
+#include "wignerSymbols.h"
+
+namespace varPhase {
+arma::colvec c_wigner3j(double l2, double l3,
+					double m1, double m2, double m3)
+{
+	// We compute the numeric limits of double precision. 
+	double huge = std::numeric_limits<double>::max();
+	double srhuge = sqrt(huge);
+	double tiny = std::numeric_limits<double>::min();
+	double srtiny = sqrt(tiny);
+	double eps = std::numeric_limits<double>::epsilon();
+
+	// We enforce the selection rules.
+	bool select(true);
+	select = (
+		   std::fabs(m1+m2+m3)<eps
+		&& std::fabs(m2) <= l2+eps
+		&& std::fabs(m3) <= l3+eps
+		);
+
+	if (!select) return arma::zeros<arma::colvec>(1);
+
+	// We compute the limits of l1. 
+	double l1min = std::max(std::fabs(l2-l3),std::fabs(m1));
+	double l1max = l2+l3;
+
+	// We compute the size of the resulting array.
+	int size = (int)std::floor(l1max-l1min+1.0+eps);
+	arma::colvec thrcof(size);
+
+	// If l1min=l1max, we have an analytical formula.
+	if (size==1)
+	{
+		thrcof(0) = pow(-1.0,std::floor(std::fabs(l2+m2-l3+m3)))/sqrt(l1min+l2+l3+1.0);
+	}
+
+	// Another special case where the recursion relation fails.
+	else
+	{
+		// We start with an arbitrary value. 
+		thrcof(0) = srtiny;
+
+		// From now on, we check the variation of |alpha(l1)|. 
+		double alphaOld, alphaNew, beta, l1(l1min);
+		if (l1min==0.0)
+			alphaNew = -(m3-m2+2.0*c_wigner3j_auxB(l1,l2,l3,m1,m2,m3))/c_wigner3j_auxA(1.0,l2,l3,m1,m2,m3);
+		else
+			alphaNew = -c_wigner3j_auxB(l1min,l2,l3,m1,m2,m3)
+							/(l1min*c_wigner3j_auxA(l1min+1.0,l2,l3,m1,m2,m3));
+
+		// We compute the two-term recursion.
+		thrcof(1) = alphaNew*thrcof(0);
+
+		// If size > 2, we start the forward recursion.
+		if (size>2)
+		{
+			// We start with an arbitrary value. 
+			thrcof(0) = srtiny;
+	
+			// From now on, we check the variation of |alpha(l1)|. 
+			double alphaOld, alphaNew, beta, l1(l1min);
+			if (l1min==0.0)
+				alphaNew = -(m3-m2+2.0*c_wigner3j_auxB(l1,l2,l3,m1,m2,m3))/c_wigner3j_auxA(1.0,l2,l3,m1,m2,m3);
+			else
+				alphaNew = -c_wigner3j_auxB(l1min,l2,l3,m1,m2,m3)
+							/(l1min*c_wigner3j_auxA(l1min+1.0,l2,l3,m1,m2,m3));
+	
+			// We compute the two-term recursion.
+			thrcof(1) = alphaNew*thrcof(0);
+	
+			// We compute the rest of the recursion.
+			int i = 1;
+			bool alphaVar = false;
+			do
+			{
+				// Bookkeeping:
+				i++;					// Next term in recursion
+				alphaOld = alphaNew;	// Monitoring of |alpha(l1)|.
+				l1 += 1.0;				// l1 = l1+1
+	
+				// New coefficients in recursion. 
+				alphaNew = -c_wigner3j_auxB(l1,l2,l3,m1,m2,m3)
+							/(l1*c_wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3));
+	
+				beta = -(l1+1.0)*c_wigner3j_auxA(l1,l2,l3,m1,m2,m3)
+						/(l1*c_wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3));
+	
+				// Application of the recursion. 
+				thrcof(i) = alphaNew*thrcof(i-1)+beta*thrcof(i-2);
+	
+				// We check if we are overflowing.
+				if (std::fabs(thrcof(i))>srhuge)
+				{
+					std::cout << "We renormalized the forward recursion." << std::endl;
+					thrcof(arma::span(0,i)) /= srhuge;
+				}
+	
+				// This piece of code checks whether we have reached
+				// the classical region. If we have, the second if 
+				// sets alphaVar to true and we break this loop at the
+				// next iteration because we need thrcof(l1mid+1) to 
+				// compute the scalar.
+				if (alphaVar) break;
+
+				if (std::fabs(alphaNew)-std::fabs(alphaOld)>0.0)
+					alphaVar=true;
+	
+			}	while (i<(size-1));	// Loop stops when we have computed all values.
+	
+			// If this is the case, we have stumbled upon a classical region.
+			// We start the backwards recursion.
+			if (i!=size-1)
+			{
+				// We keep the two terms around l1mid to compute the factor later.
+				double l1midm1(thrcof(i-2)),l1mid(thrcof(i-1)),l1midp1(thrcof(i));
+	
+				// We compute the backward recursion by providing an arbitrary 
+				// startint value.
+				thrcof(size-1) = srtiny;
+	
+				// We compute the two-term recursion.
+				l1 = l1max;
+				alphaNew = -c_wigner3j_auxB(l1,l2,l3,m1,m2,m3)
+								/((l1+1.0)*c_wigner3j_auxA(l1,l2,l3,m1,m2,m3));
+				thrcof(size-2) = alphaNew*thrcof(size-1);
+	
+				// We compute the rest of the backward recursion.
+				int j = size-2;
+				do
+				{
+					// Bookkeeping
+					j--;			// Previous term in recursion. 
+					l1 -= 1.0;		// l1 = l1-1
+	
+					// New coefficients in recursion. 
+					alphaNew = -c_wigner3j_auxB(l1,l2,l3,m1,m2,m3)
+								/((l1+1.0)*c_wigner3j_auxA(l1,l2,l3,m1,m2,m3));
+					beta = -l1*c_wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3)
+									/((l1+1.0)*c_wigner3j_auxA(l1,l2,l3,m1,m2,m3));
+	
+					// Application of the recursion. 
+					thrcof(j) = alphaNew*thrcof(j+1)+beta*thrcof(j+2);
+	
+					// We check if we are overflowing.
+					if (std::fabs(thrcof(j)>srhuge))
+					{
+						std::cout << "We renormalized the backward recursion." << std::endl;
+						thrcof(arma::span(j,size-1)) /= srhuge;
+					}
+			
+				} while (j>(i-2)); // Loop stops when we are at l1=l1mid-1.
+	
+				// We now compute the scaling factor for the forward recursion.
+				double lambda = (l1midp1*thrcof(j+2)+l1mid*thrcof(j+1)+l1midm1*thrcof(j))
+										/(l1midp1*l1midp1+l1mid*l1mid+l1midm1*l1midm1);
+	
+				// We scale the forward recursion.
+				thrcof(arma::span(0,j-1)) *= lambda;
+			}
+		}
+	}
+
+	// We compute the overall factor.
+	double sum = 0.0;
+	for (int k=0;k<size;k++)
+	{
+		sum += (2.0*(l1min+k)+1.0)*thrcof(k)*thrcof(k);
+	}
+	double c1 = pow(-1.0,l2-l3-m1)*sgn(thrcof(size-1))/sqrt(sum);
+
+	return thrcof*c1;
+}
+
+double c_wigner3j(double l1, double l2, double l3, 
+					double m1, double m2, double m3)
+{
+	// We enforce the selection rules.
+	bool select(true);
+	select = ( 
+		   std::fabs(m1+m2+m3)<1.0e-10 
+		&& std::floor(l1+l2+l3)==(l1+l2+l3)
+		&& l3 >= std::fabs(l1-l2) 
+		&& l3 <= l1+l2
+		&& std::fabs(m1) <= l1
+		&& std::fabs(m2) <= l2
+		&& std::fabs(m3) <= l3
+		);
+
+	if (!select) return 0.0;
+
+	// We compute l1min and the position of the array we will want.
+	double l1min = std::max(std::fabs(l2-l3),std::fabs(m1));
+
+	// We fetch the proper value in the array.
+	int index = (int)l1-l1min;
+
+	return c_wigner3j(l2,l3,m1,m2,m3)(index);
+}
+
+arma::colvec c_wigner6j(double l2, double l3,
+					double l4, double l5, double l6)
+{
+	// We compute the numeric limits of double precision. 
+	double huge = std::numeric_limits<double>::max();
+	double srhuge = sqrt(huge);
+	double tiny = std::numeric_limits<double>::min();
+	double srtiny = sqrt(tiny);
+	double eps = std::numeric_limits<double>::epsilon();
+
+	// We enforce the selection rules. 
+	bool select(true);
+
+	// Triangle relations for the four tryads
+	select = ( 
+		std::fabs(l4-l2) <= l6 && l6 <= l4+l2
+		&& std::fabs(l4-l5) <= l3 && l3 <= l4+l5
+		);
+
+	// Sum rule of the tryads
+	select = (
+		std::floor(l4+l2+l6)==(l4+l2+l6)
+		&& std::floor(l4+l5+l3)==(l4+l5+l3)
+		);
+
+	if (!select) return arma::zeros<arma::colvec>(1);
+
+	// We compute the limits of l1.
+	double l1min = std::max(std::fabs(l2-l3),std::fabs(l5-l6));
+	double l1max = std::min(l2+l3,l5+l6);
+
+	// We compute the size of the resulting array.
+	unsigned int size = (int)std::floor(l1max-l1min+1.0+eps);
+	arma::colvec sixcof(size);
+
+	// If l1min=l1max, we have an analytical formula.
+	if (size==1)
+	{
+		sixcof(0) = 1.0/sqrt((l1min+l1min+1.0)*(l4+l4+1.0));
+		sixcof(0) *= ((int)std::floor(l2+l3+l5+l6+eps) & 1 ? -1.0 : 1.0);
+	}
+
+	// Otherwise, we start the forward recursion.
+	else
+	{
+		// We start with an arbitrary value. 
+		sixcof(0) = srtiny;
+
+		// From now on, we check the variation of |alpha(l1)|.
+		double alphaOld, alphaNew, beta, l1(l1min);
+
+		if (l1min==0)
+			alphaNew = -(l2*(l2+1.0)+l3*(l3+1.0)+l5*(l5+1.0)+l6*(l6+1.0)-2.0*l4*(l4+1.0))/c_wigner6j_auxA(1.0,l2,l3,l4,l5,l6);
+
+		else
+			alphaNew = -c_wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+							/(l1min*c_wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6));
+
+		// We compute the two-term recursion.
+		sixcof(1) = alphaNew*sixcof(0);
+
+		if (size>2)
+		{
+			// We start with an arbitrary value. 
+			sixcof(0) = srtiny;
+	
+			// From now on, we check the variation of |alpha(l1)|. 
+			double alphaOld, alphaNew, beta, l1(l1min);
+			if (l1min==0)
+				alphaNew = -(l2*(l2+1.0)+l3*(l3+1.0)+l5*(l5+1.0)+l6*(l6+1.0)-2.0*l4*(l4+1.0))/c_wigner6j_auxA(1.0,l2,l3,l4,l5,l6);
+
+			else
+				alphaNew = -c_wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+							/(l1min*c_wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6));
+	
+			// We compute the two-term recursion.
+			sixcof(1) = alphaNew*sixcof(0);
+	
+			// We compute the rest of the recursion.
+			int i = 1;
+			bool alphaVar = false;
+			do
+			{
+				// Bookkeeping:
+				i++;					// Next term in recursion
+				alphaOld = alphaNew;	// Monitoring of |alpha(l1)|.
+				l1 += 1.0;				// l1 = l1+1
+	
+				// New coefficients in recursion. 
+				alphaNew = -c_wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+							/(l1*c_wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6));
+	
+				beta = -(l1+1.0)*c_wigner6j_auxA(l1,l2,l3,l4,l5,l6)
+						/(l1*c_wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6));
+	
+				// Application of the recursion. 
+				sixcof(i) = alphaNew*sixcof(i-1)+beta*sixcof(i-2);
+	
+				// We check if we are overflowing.
+				if (std::fabs(sixcof(i)>srhuge))
+				{
+					std::cout << "We renormalized the forward recursion." << std::endl;
+					sixcof(arma::span(0,i)) /= srhuge;
+				}
+	
+				// This piece of code checks whether we have reached
+				// the classical region. If we have, the second if 
+				// sets alphaVar to true and we break this loop at the
+				// next iteration because we need sixcof(l1mid+1) to 
+				// compute the scalar.
+				if (alphaVar) break;
+	
+				if (std::fabs(alphaNew)-std::fabs(alphaOld)>0.0)
+					alphaVar=true;
+	
+			}	while (i<(size-1));	// Loop stops when we have computed all values.
+	
+			// If this is the case, we have stumbled upon a classical region.
+			// We start the backwards recursion.
+			if (i!=size-1)
+			{
+				// We keep the two terms around l1mid to compute the factor later.
+				double l1midm1(sixcof(i-2)),l1mid(sixcof(i-1)),l1midp1(sixcof(i));
+	
+				// We compute the backward recursion by providing an arbitrary 
+				// startint value.
+				sixcof(size-1) = srtiny;
+	
+				// We compute the two-term recursion.
+				l1 = l1max;
+				alphaNew = -c_wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+								/((l1+1.0)*c_wigner6j_auxA(l1,l2,l3,l4,l5,l6));
+				sixcof(size-2) = alphaNew*sixcof(size-1);
+	
+				// We compute the rest of the backward recursion.
+				int j = size-2;
+				do
+				{
+					// Bookkeeping
+					j--;			// Previous term in recursion. 
+					l1 -= 1.0;		// l1 = l1-1
+	
+					// New coefficients in recursion. 
+					alphaNew = -c_wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+								/((l1+1.0)*c_wigner6j_auxA(l1,l2,l3,l4,l5,l6));
+					beta = -l1*c_wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6)
+								/((l1+1.0)*c_wigner6j_auxA(l1,l2,l3,l4,l5,l6));
+	
+					// Application of the recursion. 
+					sixcof(j) = alphaNew*sixcof(j+1)+beta*sixcof(j+2);
+	
+					// We check if we are overflowing.
+					if (std::fabs(sixcof(j)>srhuge))
+					{
+						std::cout << "We renormalized the backward recursion." << std::endl;
+						sixcof(arma::span(j,size-1)) /= srhuge;
+					}
+			
+				} while (j>(i-2)); // Loop stops when we are at l1=l1mid-1.
+	
+				// We now compute the scaling factor for the forward recursion.
+				double lambda = (l1midp1*sixcof(j+2)+l1mid*sixcof(j+1)+l1midm1*sixcof(j))
+									/(l1midp1*l1midp1+l1mid*l1mid+l1midm1*l1midm1);
+	
+				// We scale the forward recursion.
+				sixcof(arma::span(0,j-1)) *= lambda;
+			}
+		}
+	}
+
+	// We compute the overall factor.
+	double sum = 0.0;
+	for (int k=0;k<size;k++)
+	{
+		sum += (2.0*(l1min+k)+1.0)*(2.0*l4+1.0)*sixcof(k)*sixcof(k);
+	}
+	double c1 = pow(-1.0,std::floor(l2+l3+l5+l6+eps))*sgn(sixcof(size-1))/sqrt(sum);
+
+	return sixcof*c1;
+}
+
+double c_wigner6j(double l1, double l2, double l3, 
+					double l4, double l5, double l6)
+{
+	// We enforce the selection rules. 
+	bool select(true);
+
+	// Triangle relations for the four tryads
+	select = ( 
+		   std::fabs(l1-l2) <= l3 && l3 <= l1+l2
+		&& std::fabs(l1-l5) <= l6 && l6 <= l1+l5
+		&& std::fabs(l4-l2) <= l6 && l6 <= l4+l2
+		&& std::fabs(l4-l5) <= l3 && l3 <= l4+l5
+		);
+
+	// Sum rule of the tryads
+	select = (
+		   std::floor(l1+l2+l3)==(l1+l2+l3)
+		&& std::floor(l1+l5+l6)==(l1+l5+l6)
+		&& std::floor(l4+l2+l6)==(l4+l2+l6)
+		&& std::floor(l4+l5+l3)==(l4+l5+l3)
+		);
+
+	if (!select) return 0.0;
+
+	// We compute l1min and the position of the array we will want.
+	double l1min = std::max(std::fabs(l2-l3),std::fabs(l5-l6));
+	int index = (int)l1-l1min;
+
+	return c_wigner6j(l2,l3,l4,l5,l6)(index);
+}
+
+double c_wigner3j_auxA(double l1, double l2, double l3,
+						double m1, double m2, double m3)
+{
+	double T1 = l1*l1-pow(l2-l3,2.0);
+	double T2 = pow(l2+l3+1.0,2.0)-l1*l1;
+	double T3 = l1*l1-m1*m1;
+
+	return sqrt(T1*T2*T3);
+}
+
+double c_wigner3j_auxB(double l1, double l2, double l3, 
+						double m1, double m2, double m3)
+{
+	double T1 = -(2.0*l1+1.0);
+	double T2 = l2*(l2+1.0)*m1;
+	double T3 = l3*(l3+1.0)*m1;
+	double T4 = l1*(l1+1.0)*(m3-m2);
+
+	return T1*(T2-T3-T4);
+}
+
+double c_wigner6j_auxA(double l1, double l2, double l3,
+						double l4, double l5, double l6)
+{
+	double T1 = l1*l1-pow(l2-l3,2.0);
+	double T2 = pow(l2+l3+1.0,2.0)-l1*l1;
+	double T3 = l1*l1-pow(l5-l6,2.0);
+	double T4 = pow(l5+l6+1.0,2.0)-l1*l1;
+
+	return sqrt(T1*T2*T3*T4);
+}
+
+double c_wigner6j_auxB(double l1, double l2, double l3,
+						double l4, double l5, double l6)
+{
+	double T0 = 2.0*l1+1.0;
+
+	double T1 = l1*(l1+1.0);
+	double T2 = -l1*(l1+1.0)+l2*(l2+1.0)+l3*(l3+1.0);
+
+	double T3 = l5*(l5+1.0);
+	double T4 = l1*(l1+1.0)+l2*(l2+1.0)-l3*(l3+1.0);
+
+	double T5 = l6*(l6+1.0);
+	double T6 = l1*(l1+1.0)-l2*(l2+1.0)+l3*(l3+1.0);
+
+	double T7 = 2.0*l1*(l1+1.0)*l4*(l4+1.0);
+
+	return (T0*(T1*T2+T3*T4+T5*T6-T7));
+}
+}
\ No newline at end of file
diff --git a/wignerSymbols.f b/wignerSymbols.f
new file mode 100644
index 0000000..80add9b
--- /dev/null
+++ b/wignerSymbols.f
@@ -0,0 +1,867 @@
+*DECK DRC3JJ
+      SUBROUTINE DRC3JJ (L2, L3, M2, M3, L1MIN, L1MAX, THRCOF, NDIM,
+     +   IER)
+C***BEGIN PROLOGUE  DRC3JJ
+C***PURPOSE  Evaluate the 3j symbol f(L1) = (  L1   L2 L3)
+C                                           (-M2-M3 M2 M3)
+C            for all allowed values of L1, the other parameters
+C            being held fixed.
+C***LIBRARY   SLATEC
+C***CATEGORY  C19
+C***TYPE      DOUBLE PRECISION (RC3JJ-S, DRC3JJ-D)
+C***KEYWORDS  3J COEFFICIENTS, 3J SYMBOLS, CLEBSCH-GORDAN COEFFICIENTS,
+C             RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS,
+C             WIGNER COEFFICIENTS
+C***AUTHOR  Gordon, R. G., Harvard University
+C           Schulten, K., Max Planck Institute
+C***DESCRIPTION
+C
+C *Usage:
+C
+C        DOUBLE PRECISION L2, L3, M2, M3, L1MIN, L1MAX, THRCOF(NDIM)
+C        INTEGER NDIM, IER
+C
+C        CALL DRC3JJ (L2, L3, M2, M3, L1MIN, L1MAX, THRCOF, NDIM, IER)
+C
+C *Arguments:
+C
+C     L2 :IN      Parameter in 3j symbol.
+C
+C     L3 :IN      Parameter in 3j symbol.
+C
+C     M2 :IN      Parameter in 3j symbol.
+C
+C     M3 :IN      Parameter in 3j symbol.
+C
+C     L1MIN :OUT  Smallest allowable L1 in 3j symbol.
+C
+C     L1MAX :OUT  Largest allowable L1 in 3j symbol.
+C
+C     THRCOF :OUT Set of 3j coefficients generated by evaluating the
+C                 3j symbol for all allowed values of L1.  THRCOF(I)
+C                 will contain f(L1MIN+I-1), I=1,2,...,L1MAX+L1MIN+1.
+C
+C     NDIM :IN    Declared length of THRCOF in calling program.
+C
+C     IER :OUT    Error flag.
+C                 IER=0 No errors.
+C                 IER=1 Either L2.LT.ABS(M2) or L3.LT.ABS(M3).
+C                 IER=2 Either L2+ABS(M2) or L3+ABS(M3) non-integer.
+C                 IER=3 L1MAX-L1MIN not an integer.
+C                 IER=4 L1MAX less than L1MIN.
+C                 IER=5 NDIM less than L1MAX-L1MIN+1.
+C
+C *Description:
+C
+C     Although conventionally the parameters of the vector addition
+C  coefficients satisfy certain restrictions, such as being integers
+C  or integers plus 1/2, the restrictions imposed on input to this
+C  subroutine are somewhat weaker. See, for example, Section 27.9 of
+C  Abramowitz and Stegun or Appendix C of Volume II of A. Messiah.
+C  The restrictions imposed by this subroutine are
+C       1. L2 .GE. ABS(M2) and L3 .GE. ABS(M3);
+C       2. L2+ABS(M2) and L3+ABS(M3) must be integers;
+C       3. L1MAX-L1MIN must be a non-negative integer, where
+C          L1MAX=L2+L3 and L1MIN=MAX(ABS(L2-L3),ABS(M2+M3)).
+C  If the conventional restrictions are satisfied, then these
+C  restrictions are met.
+C
+C     The user should be cautious in using input parameters that do
+C  not satisfy the conventional restrictions. For example, the
+C  the subroutine produces values of
+C       f(L1) = ( L1  2.5  5.8)
+C               (-0.3 1.5 -1.2)
+C  for L1=3.3,4.3,...,8.3 but none of the symmetry properties of the 3j
+C  symbol, set forth on page 1056 of Messiah, is satisfied.
+C
+C     The subroutine generates f(L1MIN), f(L1MIN+1), ..., f(L1MAX)
+C  where L1MIN and L1MAX are defined above. The sequence f(L1) is
+C  generated by a three-term recurrence algorithm with scaling to
+C  control overflow. Both backward and forward recurrence are used to
+C  maintain numerical stability. The two recurrence sequences are
+C  matched at an interior point and are normalized from the unitary
+C  property of 3j coefficients and Wigner's phase convention.
+C
+C    The algorithm is suited to applications in which large quantum
+C  numbers arise, such as in molecular dynamics.
+C
+C***REFERENCES  1. Abramowitz, M., and Stegun, I. A., Eds., Handbook
+C                  of Mathematical Functions with Formulas, Graphs
+C                  and Mathematical Tables, NBS Applied Mathematics
+C                  Series 55, June 1964 and subsequent printings.
+C               2. Messiah, Albert., Quantum Mechanics, Volume II,
+C                  North-Holland Publishing Company, 1963.
+C               3. Schulten, Klaus and Gordon, Roy G., Exact recursive
+C                  evaluation of 3j and 6j coefficients for quantum-
+C                  mechanical coupling of angular momenta, J Math
+C                  Phys, v 16, no. 10, October 1975, pp. 1961-1970.
+C               4. Schulten, Klaus and Gordon, Roy G., Semiclassical
+C                  approximations to 3j  and 6j coefficients for
+C                  quantum-mechanical coupling of angular momenta,
+C                  J Math Phys, v 16, no. 10, October 1975,
+C                  pp. 1971-1988.
+C               5. Schulten, Klaus and Gordon, Roy G., Recursive
+C                  evaluation of 3j and 6j coefficients, Computer
+C                  Phys Comm, v 11, 1976, pp. 269-278.
+C***ROUTINES CALLED  XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   750101  DATE WRITTEN
+C   880515  SLATEC prologue added by G. C. Nielson, NBS; parameters
+C           HUGE and TINY revised to depend on D1MACH.
+C   891229  Prologue description rewritten; other prologue sections
+C           revised; LMATCH (location of match point for recurrences)
+C           removed from argument list; argument IER changed to serve
+C           only as an error flag (previously, in cases without error,
+C           it returned the number of scalings); number of error codes
+C           increased to provide more precise error information;
+C           program comments revised; SLATEC error handler calls
+C           introduced to enable printing of error messages to meet
+C           SLATEC standards. These changes were done by D. W. Lozier,
+C           M. A. McClain and J. M. Smith of the National Institute
+C           of Standards and Technology, formerly NBS.
+C   910415  Mixed type expressions eliminated; variable C1 initialized;
+C           description of THRCOF expanded. These changes were done by
+C           D. W. Lozier.
+C***END PROLOGUE  DRC3JJ
+C
+      INTEGER NDIM, IER
+      DOUBLE PRECISION L2, L3, M2, M3, L1MIN, L1MAX, THRCOF(NDIM)
+C
+      INTEGER I, INDEX, LSTEP, N, NFIN, NFINP1, NFINP2, NFINP3, NLIM,
+     +        NSTEP2
+      DOUBLE PRECISION A1, A1S, A2, A2S, C1, C1OLD, C2, CNORM, HVAL, D1MACH,
+     +                 DENOM, DV, EPS, HUGE, L1, M1, NEWFAC, OLDFAC,
+     +                 ONE, RATIO, SIGN1, SIGN2, SRHUGE, SRTINY, SUM1,
+     +                 SUM2, SUMBAC, SUMFOR, SUMUNI, THREE, THRESH,
+     +                 TINY, TWO, X, X1, X2, X3, Y, Y1, Y2, Y3, ZERO
+C
+      DATA  ZERO,EPS,ONE,TWO,THREE /0.0D0,0.01D0,1.0D0,2.0D0,3.0D0/
+C
+C***FIRST EXECUTABLE STATEMENT  DRC3JJ
+      IER=0
+C  HUGE is the square root of one twentieth of the largest floating
+C  point number, approximately.
+      HVAL = HUGE(1.D0)
+      SRHUGE = SQRT(HVAL)
+      TINY = 1.0D0/HVAL
+      SRTINY = 1.0D0/SRHUGE
+C
+C     LMATCH = ZERO
+      M1 = - M2 - M3
+C
+C  Check error conditions 1 and 2.
+C      IF((L2-ABS(M2)+EPS.LT.ZERO).OR.
+C     +   (L3-ABS(M3)+EPS.LT.ZERO))THEN
+C         IER=1
+C         CALL XERMSG('SLATEC','DRC3JJ','L2-ABS(M2) or L3-ABS(M3) '//
+C     +      'less than zero.',IER,1)
+C         RETURN
+C      ELSEIF((MOD(L2+ABS(M2)+EPS,ONE).GE.EPS+EPS).OR.
+C     +   (MOD(L3+ABS(M3)+EPS,ONE).GE.EPS+EPS))THEN
+C         IER=2
+C         CALL XERMSG('SLATEC','DRC3JJ','L2+ABS(M2) or L3+ABS(M3) '//
+C     +      'not integer.',IER,1)
+C         RETURN
+C      ENDIF
+C
+C
+C
+C  Limits for L1
+C
+      L1MIN = MAX(ABS(L2-L3),ABS(M1))
+      L1MAX = L2 + L3
+C
+C  Check error condition 3.
+C      IF(MOD(L1MAX-L1MIN+EPS,ONE).GE.EPS+EPS)THEN
+C         IER=3
+C         CALL XERMSG('SLATEC','DRC3JJ','L1MAX-L1MIN not integer.',IER,1)
+C         RETURN
+C      ENDIF
+      IF(L1MIN.LT.L1MAX-EPS)   GO TO 20
+      IF(L1MIN.LT.L1MAX+EPS)   GO TO 10
+C
+C  Check error condition 4.
+C      IER=4
+C      CALL XERMSG('SLATEC','DRC3JJ','L1MIN greater than L1MAX.',IER,1)
+C      RETURN
+C
+C  This is reached in case that L1 can take only one value,
+C  i.e. L1MIN = L1MAX
+C
+   10 CONTINUE
+C     LSCALE = 0
+      THRCOF(1) = (-ONE) ** INT(ABS(L2+M2-L3+M3)+EPS) /
+     1 SQRT(L1MIN + L2 + L3 + ONE)
+      RETURN
+C
+C  This is reached in case that L1 takes more than one value,
+C  i.e. L1MIN < L1MAX.
+C
+   20 CONTINUE
+C     LSCALE = 0
+      NFIN = INT(L1MAX-L1MIN+ONE+EPS)
+C      IF(NDIM-NFIN)  21, 23, 23
+C
+C  Check error condition 5.
+C   21 IER = 5
+C      CALL XERMSG('SLATEC','DRC3JJ','Dimension of result array for '//
+C     +            '3j coefficients too small.',IER,1)
+C      RETURN
+C
+C
+C  Starting forward recursion from L1MIN taking NSTEP1 steps
+C
+   23 L1 = L1MIN
+      NEWFAC = 0.0D0
+      C1 = 0.0D0
+      THRCOF(1) = SRTINY
+      SUM1 = (L1+L1+ONE) * TINY
+C
+C
+      LSTEP = 1
+   30 LSTEP = LSTEP + 1
+      L1 = L1 + ONE
+C
+C
+      OLDFAC = NEWFAC
+      A1 = (L1+L2+L3+ONE) * (L1-L2+L3) * (L1+L2-L3) * (-L1+L2+L3+ONE)
+      A2 = (L1+M1) * (L1-M1)
+      NEWFAC = SQRT(A1*A2)
+      IF(L1.LT.ONE+EPS)   GO TO 40
+C
+C
+      DV = - L2*(L2+ONE) * M1 + L3*(L3+ONE) * M1 + L1*(L1-ONE) * (M3-M2)
+      DENOM = (L1-ONE) * NEWFAC
+C
+      IF(LSTEP-2)  32, 32, 31
+C
+   31 C1OLD = ABS(C1)
+   32 C1 = - (L1+L1-ONE) * DV / DENOM
+      GO TO 50
+C
+C  If L1 = 1, (L1-1) has to be factored out of DV, hence
+C
+   40 C1 = - (L1+L1-ONE) * L1 * (M3-M2) / NEWFAC
+C
+   50 IF(LSTEP.GT.2)   GO TO 60
+C
+C
+C  If L1 = L1MIN + 1, the third term in the recursion equation vanishes,
+C  hence
+      X = SRTINY * C1
+      THRCOF(2) = X
+      SUM1 = SUM1 + TINY * (L1+L1+ONE) * C1*C1
+      IF(LSTEP.EQ.NFIN)   GO TO 220
+      GO TO 30
+C
+C
+   60 C2 = - L1 * OLDFAC / DENOM
+C
+C  Recursion to the next 3j coefficient X
+C
+      X = C1 * THRCOF(LSTEP-1) + C2 * THRCOF(LSTEP-2)
+      THRCOF(LSTEP) = X
+      SUMFOR = SUM1
+      SUM1 = SUM1 + (L1+L1+ONE) * X*X
+      IF(LSTEP.EQ.NFIN)   GO TO 100
+C
+C  See if last unnormalized 3j coefficient exceeds SRHUGE
+C
+      IF(ABS(X).LT.SRHUGE)   GO TO 80
+C
+C  This is reached if last 3j coefficient larger than SRHUGE,
+C  so that the recursion series THRCOF(1), ... , THRCOF(LSTEP)
+C  has to be rescaled to prevent overflow
+C
+C     LSCALE = LSCALE + 1
+      DO 70 I=1,LSTEP
+      IF(ABS(THRCOF(I)).LT.SRTINY)   THRCOF(I) = ZERO
+   70 THRCOF(I) = THRCOF(I) / SRHUGE
+      SUM1 = SUM1 / HVAL
+      SUMFOR = SUMFOR / HVAL
+      X = X / SRHUGE
+C
+C  As long as ABS(C1) is decreasing, the recursion proceeds towards
+C  increasing 3j values and, hence, is numerically stable.  Once
+C  an increase of ABS(C1) is detected, the recursion direction is
+C  reversed.
+C
+   80 IF(C1OLD-ABS(C1))   100, 100, 30
+C
+C
+C  Keep three 3j coefficients around LMATCH for comparison with
+C  backward recursion.
+C
+  100 CONTINUE
+C     LMATCH = L1 - 1
+      X1 = X
+      X2 = THRCOF(LSTEP-1)
+      X3 = THRCOF(LSTEP-2)
+      NSTEP2 = NFIN - LSTEP + 3
+C
+C
+C
+C
+C  Starting backward recursion from L1MAX taking NSTEP2 steps, so
+C  that forward and backward recursion overlap at three points
+C  L1 = LMATCH+1, LMATCH, LMATCH-1.
+C
+      NFINP1 = NFIN + 1
+      NFINP2 = NFIN + 2
+      NFINP3 = NFIN + 3
+      L1 = L1MAX
+      THRCOF(NFIN) = SRTINY
+      SUM2 = TINY * (L1+L1+ONE)
+C374302493
+      L1 = L1 + TWO
+      LSTEP = 1
+  110 LSTEP = LSTEP + 1
+      L1 = L1 - ONE
+C
+      OLDFAC = NEWFAC
+      A1S = (L1+L2+L3)*(L1-L2+L3-ONE)*(L1+L2-L3-ONE)*(-L1+L2+L3+TWO)
+      A2S = (L1+M1-ONE) * (L1-M1-ONE)
+      NEWFAC = SQRT(A1S*A2S)
+C
+      DV = - L2*(L2+ONE) * M1 + L3*(L3+ONE) * M1 + L1*(L1-ONE) * (M3-M2)
+C
+      DENOM = L1 * NEWFAC
+      C1 = - (L1+L1-ONE) * DV / DENOM
+      IF(LSTEP.GT.2)   GO TO 120
+C
+C  If L1 = L1MAX + 1, the third term in the recursion formula vanishes
+C
+      Y = SRTINY * C1
+      THRCOF(NFIN-1) = Y
+      SUMBAC = SUM2
+      SUM2 = SUM2 + TINY * (L1+L1-THREE) * C1*C1
+C
+      GO TO 110
+C
+C
+  120 C2 = - (L1 - ONE) * OLDFAC / DENOM
+C
+C  Recursion to the next 3j coefficient Y
+C
+      Y = C1 * THRCOF(NFINP2-LSTEP) + C2 * THRCOF(NFINP3-LSTEP)
+C
+      IF(LSTEP.EQ.NSTEP2)   GO TO 200
+C
+      THRCOF(NFINP1-LSTEP) = Y
+      SUMBAC = SUM2
+      SUM2 = SUM2 + (L1+L1-THREE) * Y*Y
+C
+C  See if last unnormalized 3j coefficient exceeds SRHUGE
+C
+      IF(ABS(Y).LT.SRHUGE)   GO TO 110
+C
+C  This is reached if last 3j coefficient larger than SRHUGE,
+C  so that the recursion series THRCOF(NFIN), ... ,THRCOF(NFIN-LSTEP+1)
+C  has to be rescaled to prevent overflow
+C
+C     LSCALE = LSCALE + 1
+      DO 130 I=1,LSTEP
+      INDEX = NFIN - I + 1
+      IF(ABS(THRCOF(INDEX)).LT.SRTINY)   THRCOF(INDEX) = ZERO
+  130 THRCOF(INDEX) = THRCOF(INDEX) / SRHUGE
+      SUM2 = SUM2 / HVAL
+      SUMBAC = SUMBAC / HVAL
+C
+C
+      GO TO 110
+C
+C
+C  The forward recursion 3j coefficients X1, X2, X3 are to be matched
+C  with the corresponding backward recursion values Y1, Y2, Y3.
+C
+  200 Y3 = Y
+      Y2 = THRCOF(NFINP2-LSTEP)
+      Y1 = THRCOF(NFINP3-LSTEP)
+C
+C
+C  Determine now RATIO such that YI = RATIO * XI  (I=1,2,3) holds
+C  with minimal error.
+C
+      RATIO = ( X1*Y1 + X2*Y2 + X3*Y3 ) / ( X1*X1 + X2*X2 + X3*X3 )
+      NLIM = NFIN - NSTEP2 + 1
+C
+      IF(ABS(RATIO).LT.ONE)   GO TO 211
+C
+      DO 210 N=1,NLIM
+  210 THRCOF(N) = RATIO * THRCOF(N)
+      SUMUNI = RATIO * RATIO * SUMFOR + SUMBAC
+      GO TO 230
+C
+  211 NLIM = NLIM + 1
+      RATIO = ONE / RATIO
+      DO 212 N=NLIM,NFIN
+  212 THRCOF(N) = RATIO * THRCOF(N)
+      SUMUNI = SUMFOR + RATIO*RATIO*SUMBAC
+      GO TO 230
+C
+  220 SUMUNI = SUM1
+C
+C
+C  Normalize 3j coefficients
+C
+  230 CNORM = ONE / SQRT(SUMUNI)
+C
+C  Sign convention for last 3j coefficient determines overall phase
+C
+      SIGN1 = SIGN(ONE,THRCOF(NFIN))
+      SIGN2 = (-ONE) ** INT(ABS(L2+M2-L3+M3)+EPS)
+      IF(SIGN1*SIGN2) 235,235,236
+  235 CNORM = - CNORM
+C
+  236 IF(ABS(CNORM).LT.ONE)   GO TO 250
+C
+      DO 240 N=1,NFIN
+  240 THRCOF(N) = CNORM * THRCOF(N)
+      RETURN
+C
+  250 THRESH = TINY / ABS(CNORM)
+      DO 251 N=1,NFIN
+      IF(ABS(THRCOF(N)).LT.THRESH)   THRCOF(N) = ZERO
+  251 THRCOF(N) = CNORM * THRCOF(N)
+C
+      RETURN
+      END
+*DECK DRC6J
+      SUBROUTINE DRC6J (L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF, NDIM,
+     +   IER)
+C***BEGIN PROLOGUE  DRC6J
+C***PURPOSE  Evaluate the 6j symbol h(L1) = {L1 L2 L3}
+C                                           {L4 L5 L6}
+C            for all allowed values of L1, the other parameters
+C            being held fixed.
+C***LIBRARY   SLATEC
+C***CATEGORY  C19
+C***TYPE      DOUBLE PRECISION (RC6J-S, DRC6J-D)
+C***KEYWORDS  6J COEFFICIENTS, 6J SYMBOLS, CLEBSCH-GORDAN COEFFICIENTS,
+C             RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS,
+C             WIGNER COEFFICIENTS
+C***AUTHOR  Gordon, R. G., Harvard University
+C           Schulten, K., Max Planck Institute
+C***DESCRIPTION
+C
+C *Usage:
+C
+C        DOUBLE PRECISION L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF(NDIM)
+C        INTEGER NDIM, IER
+C
+C        CALL DRC6J(L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF, NDIM, IER)
+C
+C *Arguments:
+C
+C     L2 :IN      Parameter in 6j symbol.
+C
+C     L3 :IN      Parameter in 6j symbol.
+C
+C     L4 :IN      Parameter in 6j symbol.
+C
+C     L5 :IN      Parameter in 6j symbol.
+C
+C     L6 :IN      Parameter in 6j symbol.
+C
+C     L1MIN :OUT  Smallest allowable L1 in 6j symbol.
+C
+C     L1MAX :OUT  Largest allowable L1 in 6j symbol.
+C
+C     SIXCOF :OUT Set of 6j coefficients generated by evaluating the
+C                 6j symbol for all allowed values of L1.  SIXCOF(I)
+C                 will contain h(L1MIN+I-1), I=1,2,...,L1MAX-L1MIN+1.
+C
+C     NDIM :IN    Declared length of SIXCOF in calling program.
+C
+C     IER :OUT    Error flag.
+C                 IER=0 No errors.
+C                 IER=1 L2+L3+L5+L6 or L4+L2+L6 not an integer.
+C                 IER=2 L4, L2, L6 triangular condition not satisfied.
+C                 IER=3 L4, L5, L3 triangular condition not satisfied.
+C                 IER=4 L1MAX-L1MIN not an integer.
+C                 IER=5 L1MAX less than L1MIN.
+C                 IER=6 NDIM less than L1MAX-L1MIN+1.
+C
+C *Description:
+C
+C     The definition and properties of 6j symbols can be found, for
+C  example, in Appendix C of Volume II of A. Messiah. Although the
+C  parameters of the vector addition coefficients satisfy certain
+C  conventional restrictions, the restriction that they be non-negative
+C  integers or non-negative integers plus 1/2 is not imposed on input
+C  to this subroutine. The restrictions imposed are
+C       1. L2+L3+L5+L6 and L2+L4+L6 must be integers;
+C       2. ABS(L2-L4).LE.L6.LE.L2+L4 must be satisfied;
+C       3. ABS(L4-L5).LE.L3.LE.L4+L5 must be satisfied;
+C       4. L1MAX-L1MIN must be a non-negative integer, where
+C          L1MAX=MIN(L2+L3,L5+L6) and L1MIN=MAX(ABS(L2-L3),ABS(L5-L6)).
+C  If all the conventional restrictions are satisfied, then these
+C  restrictions are met. Conversely, if input to this subroutine meets
+C  all of these restrictions and the conventional restriction stated
+C  above, then all the conventional restrictions are satisfied.
+C
+C     The user should be cautious in using input parameters that do
+C  not satisfy the conventional restrictions. For example, the
+C  the subroutine produces values of
+C       h(L1) = { L1 2/3  1 }
+C               {2/3 2/3 2/3}
+C  for L1=1/3 and 4/3 but none of the symmetry properties of the 6j
+C  symbol, set forth on pages 1063 and 1064 of Messiah, is satisfied.
+C
+C     The subroutine generates h(L1MIN), h(L1MIN+1), ..., h(L1MAX)
+C  where L1MIN and L1MAX are defined above. The sequence h(L1) is
+C  generated by a three-term recurrence algorithm with scaling to
+C  control overflow. Both backward and forward recurrence are used to
+C  maintain numerical stability. The two recurrence sequences are
+C  matched at an interior point and are normalized from the unitary
+C  property of 6j coefficients and Wigner's phase convention.
+C
+C    The algorithm is suited to applications in which large quantum
+C  numbers arise, such as in molecular dynamics.
+C
+C***REFERENCES  1. Messiah, Albert., Quantum Mechanics, Volume II,
+C                  North-Holland Publishing Company, 1963.
+C               2. Schulten, Klaus and Gordon, Roy G., Exact recursive
+C                  evaluation of 3j and 6j coefficients for quantum-
+C                  mechanical coupling of angular momenta, J Math
+C                  Phys, v 16, no. 10, October 1975, pp. 1961-1970.
+C               3. Schulten, Klaus and Gordon, Roy G., Semiclassical
+C                  approximations to 3j and 6j coefficients for
+C                  quantum-mechanical coupling of angular momenta,
+C                  J Math Phys, v 16, no. 10, October 1975,
+C                  pp. 1971-1988.
+C               4. Schulten, Klaus and Gordon, Roy G., Recursive
+C                  evaluation of 3j and 6j coefficients, Computer
+C                  Phys Comm, v 11, 1976, pp. 269-278.
+C***ROUTINES CALLED  XERMSG
+C***REVISION HISTORY  (YYMMDD)
+C   750101  DATE WRITTEN
+C   880515  SLATEC prologue added by G. C. Nielson, NBS; parameters
+C           HUGE and TINY revised to depend on D1MACH.
+C   891229  Prologue description rewritten; other prologue sections
+C           revised; LMATCH (location of match point for recurrences)
+C           removed from argument list; argument IER changed to serve
+C           only as an error flag (previously, in cases without error,
+C           it returned the number of scalings); number of error codes
+C           increased to provide more precise error information;
+C           program comments revised; SLATEC error handler calls
+C           introduced to enable printing of error messages to meet
+C           SLATEC standards. These changes were done by D. W. Lozier,
+C           M. A. McClain and J. M. Smith of the National Institute
+C           of Standards and Technology, formerly NBS.
+C   910415  Mixed type expressions eliminated; variable C1 initialized;
+C           description of SIXCOF expanded. These changes were done by
+C           D. W. Lozier.
+C***END PROLOGUE  DRC6J
+C
+      INTEGER NDIM, IER
+      DOUBLE PRECISION L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF(NDIM)
+C
+      INTEGER I, INDEX, LSTEP, N, NFIN, NFINP1, NFINP2, NFINP3, NLIM,
+     +        NSTEP2
+      DOUBLE PRECISION A1, A1S, A2, A2S, C1, C1OLD, C2, CNORM, HVAL, D1MACH,
+     +                 DENOM, DV, EPS, HUGE, L1, NEWFAC, OLDFAC, ONE,
+     +                 RATIO, SIGN1, SIGN2, SRHUGE, SRTINY, SUM1, SUM2,
+     +                 SUMBAC, SUMFOR, SUMUNI, THREE, THRESH, TINY, TWO,
+     +                 X, X1, X2, X3, Y, Y1, Y2, Y3, ZERO
+C
+      DATA  ZERO,EPS,ONE,TWO,THREE /0.0D0,0.01D0,1.0D0,2.0D0,3.0D0/
+C
+C***FIRST EXECUTABLE STATEMENT  DRC6J
+      IER=0
+C  HUGE is the square root of one twentieth of the largest floating
+C  point number, approximately.
+      HVAL = SQRT(HUGE(1.D0)/20.0D0)
+      SRHUGE = SQRT(HVAL)
+      TINY = 1.0D0/HVAL
+      SRTINY = 1.0D0/SRHUGE
+C
+C     LMATCH = ZERO
+C
+C  Check error conditions 1, 2, and 3.
+C      IF((MOD(L2+L3+L5+L6+EPS,ONE).GE.EPS+EPS).OR.
+C     +   (MOD(L4+L2+L6+EPS,ONE).GE.EPS+EPS))THEN
+C         IER=1
+C         CALL XERMSG('SLATEC','DRC6J','L2+L3+L5+L6 or L4+L2+L6 not '//
+C     +      'integer.',IER,1)
+C         RETURN
+C      ELSEIF((L4+L2-L6.LT.ZERO).OR.(L4-L2+L6.LT.ZERO).OR.
+C     +   (-L4+L2+L6.LT.ZERO))THEN
+C         IER=2
+C         CALL XERMSG('SLATEC','DRC6J','L4, L2, L6 triangular '//
+C     +      'condition not satisfied.',IER,1)
+C         RETURN
+C      ELSEIF((L4-L5+L3.LT.ZERO).OR.(L4+L5-L3.LT.ZERO).OR.
+C     +   (-L4+L5+L3.LT.ZERO))THEN
+C         IER=3
+C         CALL XERMSG('SLATEC','DRC6J','L4, L5, L3 triangular '//
+C     +      'condition not satisfied.',IER,1)
+C         RETURN
+C      ENDIF
+C
+C  Limits for L1
+C
+      L1MIN = MAX(ABS(L2-L3),ABS(L5-L6))
+      L1MAX = MIN(L2+L3,L5+L6)
+C
+C  Check error condition 4.
+C      IF(MOD(L1MAX-L1MIN+EPS,ONE).GE.EPS+EPS)THEN
+C         IER=4
+C         CALL XERMSG('SLATEC','DRC6J','L1MAX-L1MIN not integer.',IER,1)
+C         RETURN
+C      ENDIF
+      IF(L1MIN.LT.L1MAX-EPS)   GO TO 20
+      IF(L1MIN.LT.L1MAX+EPS)   GO TO 10
+C
+C  Check error condition 5.
+C      IER=5
+C      CALL XERMSG('SLATEC','DRC6J','L1MIN greater than L1MAX.',IER,1)
+C      RETURN
+C
+C
+C  This is reached in case that L1 can take only one value
+C
+   10 CONTINUE
+C     LSCALE = 0
+      SIXCOF(1) = (-ONE) ** INT(L2+L3+L5+L6+EPS) /
+     1            SQRT((L1MIN+L1MIN+ONE)*(L4+L4+ONE))
+      RETURN
+C
+C
+C  This is reached in case that L1 can take more than one value.
+C
+   20 CONTINUE
+C     LSCALE = 0
+      NFIN = INT(L1MAX-L1MIN+ONE+EPS)
+C      IF(NDIM-NFIN)   21, 23, 23
+C
+C  Check error condition 6.
+C   21 IER = 6
+C      CALL XERMSG('SLATEC','DRC6J','Dimension of result array for 6j '//
+C     +            'coefficients too small.',IER,1)
+C      RETURN
+C
+C
+C  Start of forward recursion
+C
+   23 L1 = L1MIN
+      NEWFAC = 0.0D0
+      C1 = 0.0D0
+      SIXCOF(1) = SRTINY
+      SUM1 = (L1+L1+ONE) * TINY
+C
+      LSTEP = 1
+   30 LSTEP = LSTEP + 1
+      L1 = L1 + ONE
+C
+      OLDFAC = NEWFAC
+      A1 = (L1+L2+L3+ONE) * (L1-L2+L3) * (L1+L2-L3) * (-L1+L2+L3+ONE)
+      A2 = (L1+L5+L6+ONE) * (L1-L5+L6) * (L1+L5-L6) * (-L1+L5+L6+ONE)
+      NEWFAC = SQRT(A1*A2)
+C
+      IF(L1.LT.ONE+EPS)   GO TO 40
+C
+      DV = TWO * ( L2*(L2+ONE)*L5*(L5+ONE) + L3*(L3+ONE)*L6*(L6+ONE)
+     1           - L1*(L1-ONE)*L4*(L4+ONE) )
+     2                   - (L2*(L2+ONE) + L3*(L3+ONE) - L1*(L1-ONE))
+     3                   * (L5*(L5+ONE) + L6*(L6+ONE) - L1*(L1-ONE))
+C
+      DENOM = (L1-ONE) * NEWFAC
+C
+      IF(LSTEP-2)  32, 32, 31
+C
+   31 C1OLD = ABS(C1)
+   32 C1 = - (L1+L1-ONE) * DV / DENOM
+      GO TO 50
+C
+C  If L1 = 1, (L1 - 1) has to be factored out of DV, hence
+C
+   40 C1 = - TWO * ( L2*(L2+ONE) + L5*(L5+ONE) - L4*(L4+ONE) )
+     1 / NEWFAC
+C
+   50 IF(LSTEP.GT.2)   GO TO 60
+C
+C  If L1 = L1MIN + 1, the third term in recursion equation vanishes
+C
+      X = SRTINY * C1
+      SIXCOF(2) = X
+      SUM1 = SUM1 + TINY * (L1+L1+ONE) * C1 * C1
+C
+      IF(LSTEP.EQ.NFIN)   GO TO 220
+      GO TO 30
+C
+C
+   60 C2 = - L1 * OLDFAC / DENOM
+C
+C  Recursion to the next 6j coefficient X
+C
+      X = C1 * SIXCOF(LSTEP-1) + C2 * SIXCOF(LSTEP-2)
+      SIXCOF(LSTEP) = X
+C
+      SUMFOR = SUM1
+      SUM1 = SUM1 + (L1+L1+ONE) * X * X
+      IF(LSTEP.EQ.NFIN)   GO TO 100
+C
+C  See if last unnormalized 6j coefficient exceeds SRHUGE
+C
+      IF(ABS(X).LT.SRHUGE)   GO TO 80
+C
+C  This is reached if last 6j coefficient larger than SRHUGE,
+C  so that the recursion series SIXCOF(1), ... ,SIXCOF(LSTEP)
+C  has to be rescaled to prevent overflow
+C
+C     LSCALE = LSCALE + 1
+      DO 70 I=1,LSTEP
+      IF(ABS(SIXCOF(I)).LT.SRTINY)   SIXCOF(I) = ZERO
+   70 SIXCOF(I) = SIXCOF(I) / SRHUGE
+      SUM1 = SUM1 / HVAL
+      SUMFOR = SUMFOR / HVAL
+      X = X / SRHUGE
+C
+C
+C  As long as the coefficient ABS(C1) is decreasing, the recursion
+C  proceeds towards increasing 6j values and, hence, is numerically
+C  stable.  Once an increase of ABS(C1) is detected, the recursion
+C  direction is reversed.
+C
+   80 IF(C1OLD-ABS(C1))   100, 100, 30
+C
+C
+C  Keep three 6j coefficients around LMATCH for comparison later
+C  with backward recursion.
+C
+  100 CONTINUE
+C     LMATCH = L1 - 1
+      X1 = X
+      X2 = SIXCOF(LSTEP-1)
+      X3 = SIXCOF(LSTEP-2)
+C
+C
+C
+C  Starting backward recursion from L1MAX taking NSTEP2 steps, so
+C  that forward and backward recursion overlap at the three points
+C  L1 = LMATCH+1, LMATCH, LMATCH-1.
+C
+      NFINP1 = NFIN + 1
+      NFINP2 = NFIN + 2
+      NFINP3 = NFIN + 3
+      NSTEP2 = NFIN - LSTEP + 3
+      L1 = L1MAX
+C
+      SIXCOF(NFIN) = SRTINY
+      SUM2 = (L1+L1+ONE) * TINY
+C
+C
+      L1 = L1 + TWO
+      LSTEP = 1
+  110 LSTEP = LSTEP + 1
+      L1 = L1 - ONE
+C
+      OLDFAC = NEWFAC
+      A1S = (L1+L2+L3)*(L1-L2+L3-ONE)*(L1+L2-L3-ONE)*(-L1+L2+L3+TWO)
+      A2S = (L1+L5+L6)*(L1-L5+L6-ONE)*(L1+L5-L6-ONE)*(-L1+L5+L6+TWO)
+      NEWFAC = SQRT(A1S*A2S)
+C
+      DV = TWO * ( L2*(L2+ONE)*L5*(L5+ONE) + L3*(L3+ONE)*L6*(L6+ONE)
+     1           - L1*(L1-ONE)*L4*(L4+ONE) )
+     2                   - (L2*(L2+ONE) + L3*(L3+ONE) - L1*(L1-ONE))
+     3                   * (L5*(L5+ONE) + L6*(L6+ONE) - L1*(L1-ONE))
+C
+      DENOM = L1 * NEWFAC
+      C1 = - (L1+L1-ONE) * DV / DENOM
+      IF(LSTEP.GT.2)   GO TO 120
+C
+C  If L1 = L1MAX + 1 the third term in the recursion equation vanishes
+C
+      Y = SRTINY * C1
+      SIXCOF(NFIN-1) = Y
+      IF(LSTEP.EQ.NSTEP2)   GO TO 200
+      SUMBAC = SUM2
+      SUM2 = SUM2 + (L1+L1-THREE) * C1 * C1 * TINY
+      GO TO 110
+C
+C
+  120 C2 = - (L1-ONE) * OLDFAC / DENOM
+C
+C  Recursion to the next 6j coefficient Y
+C
+      Y = C1 * SIXCOF(NFINP2-LSTEP) + C2 * SIXCOF(NFINP3-LSTEP)
+      IF(LSTEP.EQ.NSTEP2)   GO TO 200
+      SIXCOF(NFINP1-LSTEP) = Y
+      SUMBAC = SUM2
+      SUM2 = SUM2 + (L1+L1-THREE) * Y * Y
+C
+C  See if last unnormalized 6j coefficient exceeds SRHUGE
+C
+      IF(ABS(Y).LT.SRHUGE)   GO TO 110
+C
+C  This is reached if last 6j coefficient larger than SRHUGE,
+C  so that the recursion series SIXCOF(NFIN), ... ,SIXCOF(NFIN-LSTEP+1)
+C  has to be rescaled to prevent overflow
+C
+C     LSCALE = LSCALE + 1
+      DO 130 I=1,LSTEP
+      INDEX = NFIN-I+1
+      IF(ABS(SIXCOF(INDEX)).LT.SRTINY)   SIXCOF(INDEX) = ZERO
+  130 SIXCOF(INDEX) = SIXCOF(INDEX) / SRHUGE
+      SUMBAC = SUMBAC / HVAL
+      SUM2 = SUM2 / HVAL
+C
+      GO TO 110
+C
+C
+C  The forward recursion 6j coefficients X1, X2, X3 are to be matched
+C  with the corresponding backward recursion values Y1, Y2, Y3.
+C
+  200 Y3 = Y
+      Y2 = SIXCOF(NFINP2-LSTEP)
+      Y1 = SIXCOF(NFINP3-LSTEP)
+C
+C
+C  Determine now RATIO such that YI = RATIO * XI  (I=1,2,3) holds
+C  with minimal error.
+C
+      RATIO = ( X1*Y1 + X2*Y2 + X3*Y3 ) / ( X1*X1 + X2*X2 + X3*X3 )
+      NLIM = NFIN - NSTEP2 + 1
+C
+      IF(ABS(RATIO).LT.ONE)   GO TO 211
+C
+      DO 210 N=1,NLIM
+  210 SIXCOF(N) = RATIO * SIXCOF(N)
+      SUMUNI = RATIO * RATIO * SUMFOR + SUMBAC
+      GO TO 230
+C
+  211 NLIM = NLIM + 1
+      RATIO = ONE / RATIO
+      DO 212 N=NLIM,NFIN
+  212 SIXCOF(N) = RATIO * SIXCOF(N)
+      SUMUNI = SUMFOR + RATIO*RATIO*SUMBAC
+      GO TO 230
+C
+  220 SUMUNI = SUM1
+C
+C
+C  Normalize 6j coefficients
+C
+  230 CNORM = ONE / SQRT((L4+L4+ONE)*SUMUNI)
+C
+C  Sign convention for last 6j coefficient determines overall phase
+C
+      SIGN1 = SIGN(ONE,SIXCOF(NFIN))
+      SIGN2 = (-ONE) ** INT(L2+L3+L5+L6+EPS)
+      IF(SIGN1*SIGN2) 235,235,236
+  235 CNORM = - CNORM
+C
+  236 IF(ABS(CNORM).LT.ONE)   GO TO 250
+C
+      DO 240 N=1,NFIN
+  240 SIXCOF(N) = CNORM * SIXCOF(N)
+      RETURN
+C
+  250 THRESH = TINY / ABS(CNORM)
+      DO 251 N=1,NFIN
+      IF(ABS(SIXCOF(N)).LT.THRESH)   SIXCOF(N) = ZERO
+  251 SIXCOF(N) = CNORM * SIXCOF(N)
+C
+      RETURN
+      END
diff --git a/wignerSymbols.h b/wignerSymbols.h
new file mode 100644
index 0000000..bacabb5
--- /dev/null
+++ b/wignerSymbols.h
@@ -0,0 +1,166 @@
+#ifndef WIGNER_SYMBOLS_H
+#define WIGNER_SYMBOLS_H
+
+/** \file wignerSymbols.h  
+ *
+ * 	\author Joey Dumont <joey.dumont@gmail.com>
+ *
+ * 	\since 2013-08-16
+ *
+ * 	\brief Defines utility functions for the evaluation of Wigner-3j and -6j symbols. 
+ *
+ * We use modified SLATEC (http://netlib.org/slatec) Fortran subroutines to compute 
+ * Wigner-3j and -6j symbols. We modified the subroutines so that they do not depend 
+ * d1mach, r1mach or i1mach as these are obsolete routines. They have been replaced 
+ * by intrinsics such as huge(), tiny(), epsilon() and spacing(), which are guaranteed
+ * to work. The files wignerSymbols.f90 and utils.f contain the subroutines and their
+ * dependencies. 
+ *
+ * We rely on the ISO C Binding to bind the Fortran subroutines to C++. This method
+ * of working insures that proper type casts are performed, among other things. We
+ * then use the same method as we did before (extern "C").
+ * 
+ */
+
+#include <cmath>
+#include <limits>
+#include <algorithm>
+#include <armadillo>
+
+namespace varPhase {
+/** @name Fortran Linkage
+ * We link the Fortran subroutines to our C++ code. 
+ * Because of the ISO C Binding (see file wigner_wrap.f90), 
+ * some arguments are passed by value rather than by reference. 
+ */
+///@{
+extern "C"
+{
+	/*! Wigner-3j symbols. */
+	extern void drc3jj_wrap(double,double,double,double,double*,double*,double*,int,int*);
+
+	/*! Wigner-6j symbols. */
+	extern void drc6j_wrap(double,double,double,double,double,double*,double*,double*,int,int*);
+}
+///@}
+
+/*! @name Evaluation of Wigner-3j and -6j symbols. 
+ * We implement the Fortran subroutines in C++. 
+ */
+///@{
+arma::colvec c_wigner3j(double l2, double l3,
+						double m1, double m2, double m3);
+
+double c_wigner3j(double l1, double l2, double l3,
+					double m1, double m2, double m3);
+
+double c_wigner3j_auxA(double l1, double l2, double l3,
+						double m1, double m2, double m3);
+double c_wigner3j_auxB(double l1, double l2, double l3,
+						double m1, double m2, double m3);
+
+arma::colvec c_wigner6j(double l2, double l3,
+						double l4, double l5, double l6);
+
+double c_wigner6j(double l1, double l2, double l3,
+					double l4, double l5, double l6);
+
+double c_wigner6j_auxA(double l1, double l2, double l3,
+						double l4, double l5, double l6);
+
+double c_wigner6j_auxB(double l1, double l2, double l3,
+						double l4, double l5, double l6);
+
+/*! Computes the Wigner-3j symbol for given l1,l2,l3,m1,m2,m3. We
+ * explicitly enforce the selection rules. */
+inline double wigner3j(double l1, double l2, double l3, 
+					 	double m1, double m2, double m3)
+{
+	// We enforce the selection rules.
+	bool select(true);
+	select = ( 
+		   m1+m2+m3==0.0 
+		&& std::floor(l1+l2+l3)==(l1+l2+l3)
+		&& l3 >= std::fabs(l1-l2) 
+		&& l3 <= l1+l2
+		&& std::fabs(m1) <= l1
+		&& std::fabs(m2) <= l2
+		&& std::fabs(m3) <= l3
+		);
+
+	if (!select) return 0.0;
+
+	// We compute the size of the resulting array.
+	int size = (int)std::ceil(l2+l3-std::max(std::fabs(l2-l3),std::fabs(m1)))+1;
+
+	// We prepare the output values.
+	double l1min, l1max;
+	double thrcof [size];
+	int ierr;
+
+	// External function call. 
+	drc3jj_wrap(l2,l3,m2,m3,&l1min,&l1max,thrcof,size,&ierr);
+
+	// We fetch and return the value with the proper l1 value.
+	int index = (int)l1-l1min; 
+	return thrcof[index];
+}
+
+/*! Computes the Clebsch-Gordan coefficient by relating it to the
+ * Wigner 3j symbol. It sometimes eases the notation to use the 
+ * Clebsch-Gordan coefficients directly. */
+inline double clebschGordan(double l1, double l2, double l3,
+							double m1, double m2, double m3)
+{
+	// We simply compute it via the 3j symbol. 
+	return (pow(-1.0,l1-l2+m3)*sqrt(2.0*l3+1.0)*wigner3j(l1,l2,l3,m1,m2,-m3));
+}
+
+/*! Computes the Wigner-6j symbol for given, l1, l2, l3, l4, l5, l6.
+ * We explicitly enforce the selection rules. */
+inline double wigner6j(double l1, double l2, double l3, 
+						double l4, double l5, double l6)
+{
+	// We enforce the selection rules. 
+	bool select(true);
+
+	// Triangle relations for the four tryads
+	select = ( 
+		   std::fabs(l1-l2) <= l3 && l3 <= l1+l2
+		&& std::fabs(l1-l5) <= l6 && l6 <= l1+l5
+		&& std::fabs(l4-l2) <= l6 && l6 <= l4+l2
+		&& std::fabs(l4-l5) <= l3 && l3 <= l4+l5
+		);
+
+	// Sum rule of the tryads
+	select = (
+		   std::floor(l1+l2+l3)==(l1+l2+l3)
+		&& std::floor(l1+l5+l6)==(l1+l5+l6)
+		&& std::floor(l4+l2+l6)==(l4+l2+l6)
+		&& std::floor(l4+l5+l3)==(l4+l5+l3)
+		);
+
+	if (!select) return 0.0;
+
+	// We compute the size of the resulting array. 
+	int size = (int)std::ceil(std::min(l2+l3,l5+l6)-std::max(std::fabs(l2-l3),std::fabs(l5-l6)))+1;
+
+	// We prepare the output values
+	double l1min, l1max;
+	double sixcof [size];
+	int ierr;
+
+	// External function call
+	drc6j_wrap(l2,l3,l4,l5,l6,&l1min,&l1max,sixcof,size,&ierr);
+
+	// We fetch and return the coefficient with the proper l1 value. 
+	int index = (int)l1-l1min;
+	return sixcof[index];
+}
+
+template <typename T> double sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
+}
+
+#endif  // WIGNER_SYMBOLS_H
\ No newline at end of file
diff --git a/wignerTest.f90 b/wignerTest.f90
new file mode 100644
index 0000000..81b161b
--- /dev/null
+++ b/wignerTest.f90
@@ -0,0 +1,11 @@
+	PROGRAM WIGNER
+! 	Declarations for the main program.
+	INTEGER IER, I1
+	DOUBLE PRECISION L1MIN, L1MAX, THRCOF(200), L2, L3, M1, M2, M3
+	READ *, L2, L3, M1, M2, M3
+	CALL DRC3JJ(L2,L3,M2,M3,L1MIN,L1MAX,THRCOF,200,IER)
+
+	DO 10 I = 1,15
+		PRINT *, THRCOF(I)
+10	CONTINUE
+	END
\ No newline at end of file
diff --git a/wigner_wrap.f90 b/wigner_wrap.f90
new file mode 100644
index 0000000..c6b32ed
--- /dev/null
+++ b/wigner_wrap.f90
@@ -0,0 +1,51 @@
+! ----------------------------------------------------------------
+! - Author: 		  Joey Dumont <joey.dumont@gmail.com>            -
+! - Date created:	2013-08-15							                       -
+! - Date modded:	2013-08-15                                     -
+! - Desc.:			  ISO C Binding wrapper for Wigner-3j and        -
+! - 				      -6j symbols. This will allow the use of        -
+! -					      Fortran subroutines drc3jj, drc3jm and         -
+! -					      drc6jj.                                        -
+! ----------------------------------------------------------------
+
+subroutine drc3jj_wrap(l2, l3, m2, m3, l1min, l1max, thrcof, ndim, ier) bind(C)
+
+   		use iso_c_binding
+    	implicit none
+
+    	real(c_double), value, intent(in)           :: l2, l3, m2, m3
+    	real(c_double), intent(out)                 :: l1min, l1max
+    	real(c_double), dimension(ndim), intent(out):: thrcof
+    	integer (c_int), value, intent(in)          :: ndim
+    	integer (c_int), intent(out)                :: ier
+
+    	interface
+          SUBROUTINE DRC3JJ (L2, L3, M2, M3, L1MIN, L1MAX, THRCOF, NDIM, IER)
+              INTEGER NDIM, IER
+              DOUBLE PRECISION L2, L3, M2, M3, L1MIN, L1MAX, THRCOF(NDIM)
+          end SUBROUTINE DRC3JJ
+    	end interface
+
+    	call DRC3JJ(l2, l3, m2, m3, l1min, l1max, thrcof, ndim, ier)
+
+end subroutine drc3jj_wrap
+
+subroutine drc6j_wrap(l2, l3, l4, l5, l6, l1min, l1max, sixcof, ndim, ier) bind(C)
+      use iso_c_binding
+      implicit none
+
+      real(c_double), value, intent(in)           :: l2, l3, l4, l5, l6
+      real(c_double), intent(out)                 :: l1min, l1max
+      real(c_double), dimension(ndim), intent(out):: sixcof
+      integer(c_int), value, intent(in)           :: ndim
+      integer(c_int), intent(out)                 :: ier
+
+      interface
+          SUBROUTINE DRC6J(L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF, NDIM, IER)
+              INTEGER NDIM, IER
+              DOUBLE PRECISION L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF(NDIM)
+          END SUBROUTINE DRC6J
+      end interface
+
+      call DRC6J(l2, l3, l4, l5, l6, l1min, l1max, sixcof, ndim, ier)
+end subroutine drc6j_wrap
\ No newline at end of file