diff --git a/qpms/gaunt.c b/qpms/gaunt.c
new file mode 100644
index 0000000..449c690
--- /dev/null
+++ b/qpms/gaunt.c
@@ -0,0 +1,62 @@
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+
+
+// logarithm of factorial (from basicsubs.f90)
+double lnf (double z) {
+	// expansion parameters
+	static const double v_c0[] = {
+		0.16427423239836267e5, -0.48589401600331902e5, 0.55557391003815523e5, -0.30964901015912058e5,
+	       	0.87287202992571788e4, -0.11714474574532352e4, 0.63103078123601037e2, -0.93060589791758878e0, 
+	       	0.13919002438227877e-2,-0.45006835613027859e-8, 0.13069587914063262e-9
+	};
+	static const double cp =  2.5066282746310005;
+	double a = 1.; 
+	double b = z + 10.5;
+	b = (z + 0.5) * log(b) - b;
+	for (int i = 0; i < (sizeof(v_c0) / sizeof(double)); i++) {
+		z += 1.;
+		a += v_c0[i] / z;
+	}
+
+	return b+log(cp*a);
+}
+
+// logarithm of Pochhammer function (from basicsubs.f90)
+double lpoch(double x, double n) {
+	if(fabs(n) < 1e-5) // ???
+		return 1.;
+	double sum = x+n;
+	return lnf(sum-1.) - lnf(x-1.);
+}
+
+double f_a0 (int m, int n, int mu, int nu) {
+	double logw = lnf(n+nu-m-mu) - lnf(n-m) - lnf(nu-mu);
+	double logp = lpoch(n+1, n) + lpoch(nu+1, nu) - lpoch(n+nu+1, n+nu);
+	return exp(logw+logp);
+}
+
+
+/*
+double gaunt_gevero_direct_convert(int m, int n, int mu, int nu, int qmax, double *v_aq, int *error) {
+	int v_zero[qmax] = {0};
+	*error = 0;
+
+	if(abs(m)>n || abs(mu)=nu) {
+		*error = 1;
+		fprintf(stderr, "invalid values for m, n, mu or nu\n")
+		return NAN;
+	}
+
+	switch(qmax) {
+		case 0:
+			v_aq[0]  = f_a0(m,n,mu,nu);
+			break;
+		case 1:
+			v_aq[0] = f_a0(m,n,mu,nu);
+			!!!!!!!!!TODO CONTINUE HERE
+
+
+
+*/