Split f_lattice out

1 files changed, 52 insertions, 33 deletions

diff --git a/src/diffraction.c b/src/diffraction.c
index b2da4622..902ac008 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -20,16 +20,58 @@
 #include "ewald.h"
 
 
-void get_diffraction(struct image *image, UnitCell *cell)
+static double lattice_factor(struct threevec q, double ax, double ay, double az,
+                                                double bx, double by, double bz,
+                                                double cx, double cy, double cz)
 {
-	int x, y;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;	
+	struct threevec Udotq;
+	double f1, f2, f3;
 	int na = 64;
 	int nb = 64;
 	int nc = 64;
 
+	Udotq.u = (ax*q.u + ay*q.v + az*q.w)/2.0;
+	Udotq.v = (bx*q.u + by*q.v + bz*q.w)/2.0;
+	Udotq.w = (cx*q.u + cy*q.v + cz*q.w)/2.0;
+	
+	if ( na > 1 ) {
+		f1 = sin(2.0*M_PI*(double)na*Udotq.u)
+		       / sin(2.0*M_PI*Udotq.u);
+	} else {
+		f1 = 1.0;
+	}
+	
+	if ( nb > 1 ) {
+		f2 = sin(2.0*M_PI*(double)nb*Udotq.v)
+		       / sin(2.0*M_PI*Udotq.v);
+	} else {
+		f2 = 1.0;
+	}
+	
+	if ( nc > 1 ) {
+		f3 = sin(2.0*M_PI*(double)nc*Udotq.w)
+		       / sin(2.0*M_PI*Udotq.w);
+	} else {
+		f3 = 1.0;
+	}
+	
+	return f1 * f2 * f3;
+}
+
+
+static double molecule_factor(struct threevec q)
+{
+	return 1e5;
+}
+
+
+void get_diffraction(struct image *image, UnitCell *cell)
+{
+	int x, y;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+	
 	/* Generate the array of reciprocal space vectors in image->qvecs */
 	get_ewald(image);
 	
@@ -42,38 +84,15 @@ void get_diffraction(struct image *image, UnitCell *cell)
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
 	
+		double f_lattice, f_molecule;
 		struct threevec q;
-		struct threevec Udotq;
-		double f1, f2, f3;
-		                
-		q = image->qvecs[x + image->width*y];
 		
-		Udotq.u = (ax*q.u + ay*q.v + az*q.w)/2.0;
-		Udotq.v = (bx*q.u + by*q.v + bz*q.w)/2.0;
-		Udotq.w = (cx*q.u + cy*q.v + cz*q.w)/2.0;
-		
-		if ( na > 1 ) {
-			f1 = sin(2.0*M_PI*(double)na*Udotq.u)
-			       / sin(2.0*M_PI*Udotq.u);
-		} else {
-			f1 = 1.0;
-		}
-		
-		if ( nb > 1 ) {
-			f2 = sin(2.0*M_PI*(double)nb*Udotq.v)
-			       / sin(2.0*M_PI*Udotq.v);
-		} else {
-			f2 = 1.0;
-		}
+		q = image->qvecs[x + image->width*y];
 		
-		if ( nc > 1 ) {
-			f3 = sin(2.0*M_PI*(double)nc*Udotq.w)
-			       / sin(2.0*M_PI*Udotq.w);
-		} else {
-			f3 = 1.0;
-		}
+		f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz);
+		f_molecule = molecule_factor(q);
 		
-		image->sfacs[x + image->width*y] = f1 * f2 * f3;
+		image->sfacs[x + image->width*y] = f_lattice * f_molecule;
 
 	}
 	}