Add memoization to scattering factor calculation

1 files changed, 40 insertions, 1 deletions

diff --git a/src/sfac.c b/src/sfac.c
index 342c9615..45514992 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -20,6 +20,9 @@
 #include "sfac.h"
 
 
+#define N_MEMO 1024
+
+
 /* Look up f1 and f2 for this atom at this energy (in J/photon) */
 static double complex get_f1f2(const char *n, double en)
 {
@@ -28,6 +31,17 @@ static double complex get_f1f2(const char *n, double en)
 	char line[1024];
 	char *rval;
 	double last_E, last_f1, last_f2;
+	static char *memo_n[N_MEMO];
+	static double memo_en[N_MEMO];
+	static double complex memo_res[N_MEMO];
+	static int n_memo = 0;
+	int i;
+	
+	for ( i=0; i<n_memo; i++ ) {
+		if ( (memo_en[i] == en) && (strcmp(memo_n[i], n) == 0) ) {
+			return memo_res[i];
+		}
+	}
 
 	snprintf(filename, 63, "scattering-factors/%s.nff", n);
 	fh = fopen(filename, "r");
@@ -71,6 +85,7 @@ static double complex get_f1f2(const char *n, double en)
 			/* Perform (linear) interpolation */
 			double f;
 			double actual_f1, actual_f2;
+			double complex res;
 
 			f = (en - last_E) / (E - last_E);
 
@@ -78,7 +93,15 @@ static double complex get_f1f2(const char *n, double en)
 			actual_f2 = last_f2 + f * (f2 - last_f2);
 
 			fclose(fh);
-			return actual_f1 + I*actual_f2;
+			
+			res = actual_f1 + I*actual_f2;
+			
+			memo_n[n_memo] = strdup(n);
+			memo_en[n_memo] = en;
+			memo_res[n_memo++] = res;
+			n_memo = n_memo % N_MEMO;
+			
+			return res;
 
 		}
 
@@ -101,6 +124,17 @@ static double get_waas_kirf(const char *n, double s)
 	double f;
 	float a1, a2, a3, a4, a5, c, b1, b2, b3, b4, b5;
 	double s2;
+	static char *memo_n[N_MEMO];
+	static double memo_s[N_MEMO];
+	static double memo_res[N_MEMO];
+	static int n_memo = 0;
+	int i;
+	
+	for ( i=0; i<n_memo; i++ ) {
+		if ( (memo_s[i] == s) && (strcmp(memo_n[i], n) == 0) ) {
+			return memo_res[i];
+		}
+	}
 	
 	fh = fopen("scattering-factors/f0_WaasKirf.dat", "r");
 	if ( fh == NULL ) {
@@ -145,6 +179,11 @@ static double get_waas_kirf(const char *n, double s)
 	f = c + a1*exp(-b1*s2) + a2*exp(-b2*s2) + a3*exp(-b3*s2)
 	      + a4*exp(-b4*s2) + a5*exp(-b5*s2);
 
+	memo_n[n_memo] = strdup(n);
+	memo_s[n_memo] = s;
+	memo_res[n_memo++] = f;
+	n_memo = n_memo % N_MEMO;
+
 	return f;
 }