From b235ba3004e7581fe543d9ce89d643a58bc0ae2f Mon Sep 17 00:00:00 2001
From: Joey Dumont <joey.dumont@gmail.com>
Date: Sun, 10 Nov 2013 14:07:48 -0500
Subject: [PATCH] Remove test files.

---
 src/fortranLinkage.h  |  10 -
 src/testWigner.cpp    | 221 -----------
 src/wignerSymbols.cpp | 463 ----------------------
 src/wignerSymbols.f   | 867 ------------------------------------------
 src/wignerSymbols.h   | 166 --------
 src/wignerTest.f90    |  11 -
 src/wigner_wrap.f90   |  51 ---
 7 files changed, 1789 deletions(-)
 delete mode 100644 src/fortranLinkage.h
 delete mode 100644 src/testWigner.cpp
 delete mode 100644 src/wignerSymbols.cpp
 delete mode 100644 src/wignerSymbols.f
 delete mode 100644 src/wignerSymbols.h
 delete mode 100644 src/wignerTest.f90
 delete mode 100644 src/wigner_wrap.f90

diff --git a/src/fortranLinkage.h b/src/fortranLinkage.h
deleted file mode 100644
index be4973e..0000000
--- a/src/fortranLinkage.h
+++ /dev/null
@@ -1,10 +0,0 @@
-#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/src/testWigner.cpp b/src/testWigner.cpp
deleted file mode 100644
index 5cb7d94..0000000
--- a/src/testWigner.cpp
+++ /dev/null
@@ -1,221 +0,0 @@
-#include "wignerSymbols.h"
-
-#include <iostream>
-#include <iomanip>
-#include <armadillo>
-#include <time.h>
-
-#define EPS 1.0e-10
-
-using namespace std;
-using namespace arma;
-
-// This is a test.
-// 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/src/wignerSymbols.cpp b/src/wignerSymbols.cpp
deleted file mode 100644
index da23b45..0000000
--- a/src/wignerSymbols.cpp
+++ /dev/null
@@ -1,463 +0,0 @@
-#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/src/wignerSymbols.f b/src/wignerSymbols.f
deleted file mode 100644
index 80add9b..0000000
--- a/src/wignerSymbols.f
+++ /dev/null
@@ -1,867 +0,0 @@
-*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/src/wignerSymbols.h b/src/wignerSymbols.h
deleted file mode 100644
index bacabb5..0000000
--- a/src/wignerSymbols.h
+++ /dev/null
@@ -1,166 +0,0 @@
-#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/src/wignerTest.f90 b/src/wignerTest.f90
deleted file mode 100644
index 81b161b..0000000
--- a/src/wignerTest.f90
+++ /dev/null
@@ -1,11 +0,0 @@
-	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/src/wigner_wrap.f90 b/src/wigner_wrap.f90
deleted file mode 100644
index c6b32ed..0000000
--- a/src/wigner_wrap.f90
+++ /dev/null
@@ -1,51 +0,0 @@
-! ----------------------------------------------------------------
-! - 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