Add (slow) bandwidth and multisampling

6 files changed, 125 insertions, 29 deletions

diff --git a/src/detector.c b/src/detector.c
index 2c313cdc..ccf37024 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -182,7 +182,7 @@ void record_image(struct image *image, int do_water, int do_poisson,
 		if ( do_water ) {
 
 			/* Add intensity contribution from water */
-			water = water_intensity(image->qvecs[x + image->width*y],
+			water = water_intensity(image->qvecs[0][x + image->width*y],
 			                        ph_lambda_to_en(image->lambda),
 			                        BEAM_RADIUS, WATER_RADIUS);
 			intensity += water;
diff --git a/src/diffraction.c b/src/diffraction.c
index 27037839..b9430ba7 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -174,20 +174,25 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
 		double complex f_molecule;
 		struct rvec q;
 		double complex val;
+		int i;
 
-		q = image->qvecs[x + image->width*y];
+		for ( i=0; i<image->nspheres; i++ ) {
 
-		f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz,
-		                           na, nb, nc);
-		if ( no_sfac ) {
-			f_molecule = 1.0;
-		} else {
-			f_molecule = molecule_factor(image->molecule, q,
-		                             ax,ay,az,bx,by,bz,cx,cy,cz);
-		}
+			q = image->qvecs[i][x + image->width*y];
+
+			f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz,
+			                           na, nb, nc);
+			if ( no_sfac ) {
+				f_molecule = 10000.0;
+			} else {
+				f_molecule = molecule_factor(image->molecule, q,
+			                            ax,ay,az,bx,by,bz,cx,cy,cz);
+			}
 
-		val = f_molecule * f_lattice;
-		image->sfacs[x + image->width*y] = val;
+			val = f_molecule * f_lattice;
+			image->sfacs[x + image->width*y] = val;
+
+		}
 
 	}
 	progress_bar(x, image->width-1, "Calculating lattice factors");
diff --git a/src/ewald.c b/src/ewald.c
index 94fb4ed6..f24e9273 100644
--- a/src/ewald.c
+++ b/src/ewald.c
@@ -21,6 +21,11 @@
 #include "detector.h"
 
 
+#define SAMPLING (4)
+#define BWSAMPLING (10)
+#define BANDWIDTH (0.05)
+
+
 static struct rvec quat_rot(struct rvec q, struct quaternion z)
 {
 	struct rvec res;
@@ -52,18 +57,9 @@ static struct rvec quat_rot(struct rvec q, struct quaternion z)
 }
 
 
-void get_ewald(struct image *image)
+static void add_sphere(struct image *image, double k)
 {
 	int x, y;
-	double k;  /* Wavenumber */
-
-	k = 1/image->lambda;
-
-	image->qvecs = malloc(image->width * image->height
-	                       * sizeof(struct rvec));
-
-	image->twotheta = malloc(image->width * image->height
-	                       * sizeof(double));
 
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
@@ -73,8 +69,8 @@ void get_ewald(struct image *image)
 		double r;
 		double twothetax, twothetay, twotheta;
 		double qx, qy, qz;
-		struct rvec q;
-		int p;
+		struct rvec q1, q2, q3, q4;
+		int p, sx, sy, i;
 
 		/* Calculate q vectors for Ewald sphere */
 		for ( p=0; p<image->det.n_panels; p++ ) {
@@ -86,23 +82,80 @@ void get_ewald(struct image *image)
 				                            / image->resolution;
 				ry = ((double)y - image->det.panels[p].cy)
 				                            / image->resolution;
+				break;
 			}
 		}
+
+		/* Bottom left corner */
 		r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+		twothetax = atan2(rx, image->camera_len);
+		twothetay = atan2(ry, image->camera_len);
+		twotheta = atan2(r, image->camera_len);
+		qx = k * sin(twothetax);
+		qy = k * sin(twothetay);
+		qz = k - k * cos(twotheta);
+		q1.u = qx;  q1.v = qy;  q1.w = qz;
+		/* 2theta value is calculated at the bottom left only */
+		image->twotheta[x + image->width*y] = twotheta;
 
+		/* Bottom right corner (using the same panel configuration!) */
+		rx = ((double)(x+1) - image->det.panels[p].cx)
+		                            / image->resolution;
+		ry = ((double)y - image->det.panels[p].cy)
+		                            / image->resolution;
 		twothetax = atan2(rx, image->camera_len);
 		twothetay = atan2(ry, image->camera_len);
 		twotheta = atan2(r, image->camera_len);
+		qx = k * sin(twothetax);
+		qy = k * sin(twothetay);
+		qz = k - k * cos(twotheta);
+		q2.u = qx;  q2.v = qy;  q2.w = qz;
 
+		/* Top left corner (using the same panel configuration!) */
+		rx = ((double)x - image->det.panels[p].cx)
+		                            / image->resolution;
+		ry = ((double)(y+1) - image->det.panels[p].cy)
+		                            / image->resolution;
+		twothetax = atan2(rx, image->camera_len);
+		twothetay = atan2(ry, image->camera_len);
+		twotheta = atan2(r, image->camera_len);
 		qx = k * sin(twothetax);
 		qy = k * sin(twothetay);
 		qz = k - k * cos(twotheta);
+		q3.u = qx;  q3.v = qy;  q3.w = qz;
 
-		q.u = qx;  q.v = qy;  q.w = qz;
-		image->qvecs[x + image->width*y] = quat_rot(q,
+		/* Top right corner (using the same panel configuration!) */
+		rx = ((double)(x+1) - image->det.panels[p].cx)
+		                            / image->resolution;
+		ry = ((double)(y+1) - image->det.panels[p].cy)
+		                            / image->resolution;
+		twothetax = atan2(rx, image->camera_len);
+		twothetay = atan2(ry, image->camera_len);
+		twotheta = atan2(r, image->camera_len);
+		qx = k * sin(twothetax);
+		qy = k * sin(twothetay);
+		qz = k - k * cos(twotheta);
+		q4.u = qx;  q4.v = qy;  q4.w = qz;
+
+		/* Now interpolate between the values to get
+		 * the sampling points */
+		i = 0;
+		for ( sx=0; sx<SAMPLING; sx++ ) {
+		for ( sy=0; sy<SAMPLING; sy++ ) {
+
+			struct rvec q;
+
+			q.u = q1.u + ((q2.u - q1.u)/SAMPLING)*sx
+			           + ((q3.u - q1.u)/SAMPLING)*sy;
+			q.v = q1.v + ((q2.v - q1.v)/SAMPLING)*sx
+			           + ((q3.v - q1.v)/SAMPLING)*sy;
+			q.w = q1.w + ((q2.w - q1.w)/SAMPLING)*sx
+			           + ((q3.w - q1.w)/SAMPLING)*sy;
+			image->qvecs[i++][x + image->width*y] = quat_rot(q,
 		                                            image->orientation);
 
-		image->twotheta[x + image->width*y] = twotheta;
+		}
+		}
 
 		if ( (x==0) && (y==(int)image->y_centre) ) {
 			double s;
@@ -123,4 +176,37 @@ void get_ewald(struct image *image)
 
 	}
 	}
+
+}
+
+
+void get_ewald(struct image *image)
+{
+	double kc;  /* Wavenumber */
+	int i, kstep;
+
+	kc = 1/image->lambda;  /* Centre */
+
+	image->qvecs = malloc(image->width * image->height
+	                       * sizeof(struct rvec *));
+
+	image->twotheta = malloc(image->width * image->height
+	                       * sizeof(double));
+
+	/* Create the spheres */
+	image->nspheres = SAMPLING*SAMPLING*BWSAMPLING;
+	for ( i=0; i<image->nspheres; i++ ) {
+		image->qvecs[i] = malloc(image->width * image->height
+                                                 * sizeof(struct rvec));
+	}
+
+	for ( kstep=0; kstep<BWSAMPLING; kstep++ ) {
+
+		double k;
+
+		k = kc + (kstep-(BWSAMPLING/2))*kc*(BANDWIDTH/BWSAMPLING);
+
+		add_sphere(image, k);
+
+	}
 }
diff --git a/src/image.h b/src/image.h
index 166b20df..ff87b1bc 100644
--- a/src/image.h
+++ b/src/image.h
@@ -76,7 +76,8 @@ struct image {
 	int			*hdr;      /* Actual counts */
 	int16_t			*data;     /* Integer counts after bloom */
 	double complex		*sfacs;
-	struct rvec		*qvecs;
+	struct rvec		**qvecs;
+	int			nspheres;
 	double			*twotheta;
 	struct molecule		*molecule;
 	UnitCell		*indexed_cell;
diff --git a/src/intensities.c b/src/intensities.c
index d5b3604c..edc15910 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -72,7 +72,7 @@ void output_intensities(struct image *image, UnitCell *cell)
 		int found = 0;
 		int j;
 
-		q = image->qvecs[x + image->width*y];
+		q = image->qvecs[0][x + image->width*y];
 
 		hd = q.u * ax + q.v * ay + q.w * az;
 		kd = q.u * bx + q.v * by + q.w * bz;
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index bfc6be37..70384250 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -247,6 +247,7 @@ int main(int argc, char *argv[])
 	do {
 
 		int na, nb, nc;
+		int i;
 
 		na = 8*random()/RAND_MAX + 4;
 		nb = 8*random()/RAND_MAX + 4;
@@ -302,6 +303,9 @@ int main(int argc, char *argv[])
 
 		/* Clean up */
 		free(image.data);
+		for ( i=0; i<image.nspheres; i++ ) {
+			free(image.qvecs[i]);
+		}
 		free(image.qvecs);
 		free(image.hdr);
 		free(image.sfacs);