Add Henke table lookup

3 files changed, 76 insertions, 14 deletions

diff --git a/src/detector.c b/src/detector.c
index 9afe2cf3..bfc78b6e 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -44,15 +44,16 @@ void record_image(struct image *image)
 	for ( y=0; y<image->height; y++ ) {
 	
 		double counts;
-		double val, intensity;
+		double intensity;
 		double sa;
-		
+		double complex val;
+
 		val = image->sfacs[x + image->width*y];
-		
+		intensity = val * conj(val);
+
 		/* What solid angle is subtended by this pixel? */
 		sa = sa_per_pixel * cos(image->twotheta[x + image->width*y]);
-		
-		intensity = pow(val, 2.0);
+
 		counts = intensity * ph_per_e * sa;
 		
 		image->data[x + image->width*y] = counts;
diff --git a/src/diffraction.c b/src/diffraction.c
index 91154b05..4f5a1d1d 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -61,9 +61,70 @@ static double lattice_factor(struct threevec q, double ax, double ay, double az,
 }
 
 
-static complex get_sfac(const char *n, double s, double en)
+/* Look up f1 and f2 for this atom at this energy (in J/photon) */
+static complex get_sfac(const char *n, double en)
 {
-	return 1.0;
+	FILE *fh;
+	char filename[64];
+	char line[1024];
+	char *rval;
+	float last_E, last_f1, last_f2;
+
+	snprintf(filename, 63, "scattering-factors/%s.nff", n);
+	fh = fopen(filename, "r");
+	if ( fh == NULL ) {
+		fprintf(stderr, "Couldn't open file '%s'\n", filename);
+		return 0.0;
+	}
+
+	en = J_to_eV(en);
+
+	/* Discard first line */
+	fgets(line, 1023, fh);
+
+	last_E = 0.0;
+	last_f1 = 0.0;
+	last_f2 = 0.0;
+	do {
+
+		int r;
+		float E, f1, f2;
+
+		rval = fgets(line, 1023, fh);
+
+		r = sscanf(line, "%f %f %f", &E, &f1, &f2);
+		if ( r != 3 ) {
+			fprintf(stderr, "WTF?\n");
+			abort();
+		}
+
+		if ( E < en ) {
+			last_E = E;
+			last_f1 = f1;
+			last_f2 = f2;
+		} else {
+
+			double f;
+			float actual_f1, actual_f2;
+
+			f = (en - last_E) / (E - last_E);
+
+			actual_f1 = last_f1 + f * (f1 - last_f1);
+			actual_f2 = last_f2 + f * (f2 - last_f2);
+
+			fclose(fh);
+			return actual_f1 + I*actual_f2;
+
+		}
+
+	} while ( rval != NULL );
+
+	fclose(fh);
+
+	fprintf(stderr, "Couldn't find scattering factors for '%s' at %f eV!\n",
+	        n, en);
+
+	return 0.0;
 }
 
 
@@ -77,9 +138,9 @@ static complex molecule_factor(struct molecule *mol, struct threevec q,
 	s = modulus(q.u, q.v, q.w);
 	
 	for ( i=0; i<mol->n_species; i++ ) {
-		
-		complex sfac;
-		complex contrib = 0.0;
+
+		double complex sfac;
+		double complex contrib = 0.0;
 		struct mol_species *spec;
 		int j;
 		
@@ -90,10 +151,10 @@ static complex molecule_factor(struct molecule *mol, struct threevec q,
 			double ph;
 			
 			ph= q.u*spec->x[j] + q.v*spec->y[j] + q.w*spec->z[j];
-			contrib += cos(ph) + I*sin(ph);
+			contrib += cos(2.0*M_PI*ph) + I*sin(2.0*M_PI*ph);
 		}
-		
-		sfac = get_sfac(spec->species, s, en);
+
+		sfac = get_sfac(spec->species, en);
 		F += sfac * contrib * exp(-2.0 * spec->B[j] * s);
 	}
 	
diff --git a/src/image.h b/src/image.h
index 3f67e973..7c64c6d6 100644
--- a/src/image.h
+++ b/src/image.h
@@ -81,7 +81,7 @@ struct molecule
 struct image {
 
 	uint16_t		*data;
-	complex			*sfacs;
+	double complex		*sfacs;
 	struct threevec		*qvecs;
 	double			*twotheta;
 	struct molecule		*molecule;