Merge branch 'master' of ssh://git-weiss@git.bitwiz.org.uk/srv/git-private/crystfel

4 files changed, 104 insertions, 55 deletions

diff --git a/src/detector.c b/src/detector.c
index fbe2c5f2..23378e90 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -146,7 +146,7 @@ void record_image(struct image *image)
 	int x, y;
 	double fluence;
 	double ph_per_e;
-	double sa_per_pixel;
+	double pix_area, Lsq;
 	double area;
 	const int do_bloom = 1;
 
@@ -157,10 +157,9 @@ void record_image(struct image *image)
 	printf("Fluence = %5e photons / m^2, %5e photons total\n",
 	       fluence, fluence*area);
 
-	/* Solid angle subtended by a single pixel at the centre of the CCD */
-	sa_per_pixel = pow(1.0/(image->resolution * image->camera_len), 2.0);
-	printf("Solid angle of one pixel (at centre of CCD) = %e sr\n",
-	       sa_per_pixel);
+	/* Area of one pixel */
+	pix_area = pow(1.0/image->resolution, 2.0);
+	Lsq = pow(image->camera_len, 2.0);
 
 	image->hdr = malloc(image->width * image->height * sizeof(double));
 
@@ -169,6 +168,7 @@ void record_image(struct image *image)
 
 		double counts, intensity, sa, water;
 		double complex val;
+		double dsq, proj_area;
 
 		val = image->sfacs[x + image->width*y];
 		intensity = pow(cabs(val), 2.0);
@@ -181,13 +181,20 @@ void record_image(struct image *image)
 		//printf("%e, %e, ", intensity, water);
 		intensity += water;
 
-		/* What solid angle is subtended by this pixel? */
-		sa = sa_per_pixel * cos(image->twotheta[x + image->width*y]);
+		/* Area of pixel as seen from crystal (approximate) */
+		proj_area = pix_area * cos(image->twotheta[x + image->width*y]);
+
+		/* Calculate distance from crystal to pixel */
+		dsq = pow(((double)x - image->x_centre)/image->resolution, 2.0);
+		dsq += pow(((double)y - image->y_centre)/image->resolution, 2.0);
+
+		/* Projected area of pixel divided by distance squared */
+		sa = proj_area / (dsq + Lsq);
 
 		counts = intensity * ph_per_e * sa;
 		//printf("%e counts\n", counts);
 
-		counts += poisson_noise(counts);
+		counts = poisson_noise(counts);
 
 		image->hdr[x + image->width*y] = counts;
 
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 84956bc3..e9c56fb3 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -27,10 +27,6 @@
 #include "detector.h"
 
 
-/* Crystal size in metres */
-#define CRYSTAL_SIZE (500.0e-9)
-
-
 static void show_help(const char *s)
 {
 	printf("Syntax: %s\n\n", s);
diff --git a/src/sfac.c b/src/sfac.c
index 40ae8db4..cd76407f 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -119,6 +119,21 @@ static double complex get_f1f2(const char *n, double en)
 }
 
 
+struct waas_kirf_factors {
+	char *n;
+	float a1;
+	float a2;
+	float a3;
+	float a4;
+	float a5;
+	float b1;
+	float b2;
+	float b3;
+	float b4;
+	float b5;
+	float c;
+};
+
 /* s = sin(theta)/lambda in metres^-1*/
 static double get_waas_kirf(const char *n, double s)
 {
@@ -127,66 +142,97 @@ 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];
+	static struct waas_kirf_factors waaskirfcache[N_MEMO];
+	int found;
+	static int n_waaskirfcache = 0;
+
+	found = 0;
+	for ( i=0; i<n_waaskirfcache; i++ ) {
+		if ( strcmp(waaskirfcache[i].n, n) == 0 ) {
+
+			found = 1;
+
+			a1 = waaskirfcache[n_waaskirfcache].a1;
+			a2 = waaskirfcache[n_waaskirfcache].a2;
+			a3 = waaskirfcache[n_waaskirfcache].a3;
+			a4 = waaskirfcache[n_waaskirfcache].a4;
+			a5 = waaskirfcache[n_waaskirfcache].a5;
+			b1 = waaskirfcache[n_waaskirfcache].b1;
+			b2 = waaskirfcache[n_waaskirfcache].b2;
+			b3 = waaskirfcache[n_waaskirfcache].b3;
+			b4 = waaskirfcache[n_waaskirfcache].b4;
+			b5 = waaskirfcache[n_waaskirfcache].b5;
+			c = waaskirfcache[n_waaskirfcache].c;
+
+			break;
 		}
 	}
 
-	fh = fopen("scattering-factors/f0_WaasKirf.dat", "r");
-	if ( fh == NULL ) {
-		fprintf(stderr, "Couldn't open f0_WaasKirf.dat\n");
-		return 0.0;
-	}
+	if ( !found ) {
 
-	do {
+		fh = fopen("scattering-factors/f0_WaasKirf.dat", "r");
+		if ( fh == NULL ) {
+			fprintf(stderr, "Couldn't open f0_WaasKirf.dat\n");
+			return 0.0;
+		}
 
-		int r;
-		char line[1024];
-		char sp[1024];
-		int Z;
+		do {
 
-		rval = fgets(line, 1023, fh);
+			int r;
+			char line[1024];
+			char sp[1024];
+			int Z;
 
-		if ( (line[0] != '#') || (line[1] != 'S') ) continue;
+			rval = fgets(line, 1023, fh);
 
-		r = sscanf(line, "#S  %i  %s", &Z, sp);
-		if ( (r != 2) || (strcmp(sp, n) != 0) ) continue;
+			if ( (line[0] != '#') || (line[1] != 'S') ) continue;
 
-		/* Skip two lines */
-		fgets(line, 1023, fh);
-		fgets(line, 1023, fh);
+			r = sscanf(line, "#S  %i  %s", &Z, sp);
+			if ( (r != 2) || (strcmp(sp, n) != 0) ) continue;
 
-		/* Read scattering coefficients */
-		rval = fgets(line, 1023, fh);
-		r = sscanf(line, "  %f  %f  %f  %f  %f  %f  %f %f %f %f %f",
-		          &a1, &a2, &a3, &a4, &a5, &c, &b1, &b2, &b3, &b4, &b5);
-		if ( r != 11 ) {
-			fprintf(stderr, "Couldn't read scattering factors\n");
-			return 0.0;
-		}
+			/* Skip two lines */
+			fgets(line, 1023, fh);
+			fgets(line, 1023, fh);
 
-		break;
+			/* Read scattering coefficients */
+			rval = fgets(line, 1023, fh);
+			r = sscanf(line,
+			             "  %f  %f  %f  %f  %f  %f  %f %f %f %f %f",
+			            &a1, &a2, &a3, &a4, &a5, &c,
+			            &b1, &b2, &b3, &b4, &b5);
+			if ( r != 11 ) {
+				fprintf(stderr, "Couldn't read scattering "
+				                "factors (WaasKirf)\n");
+				return 0.0;
+			}
 
-	} while ( rval != NULL );
+			break;
 
-	fclose(fh);
+		} while ( rval != NULL );
+
+		fclose(fh);
+
+		waaskirfcache[n_waaskirfcache].a1 = a1;
+		waaskirfcache[n_waaskirfcache].a2 = a2;
+		waaskirfcache[n_waaskirfcache].a3 = a3;
+		waaskirfcache[n_waaskirfcache].a4 = a4;
+		waaskirfcache[n_waaskirfcache].a5 = a5;
+		waaskirfcache[n_waaskirfcache].c = c;
+		waaskirfcache[n_waaskirfcache].b1 = b1;
+		waaskirfcache[n_waaskirfcache].b2 = b2;
+		waaskirfcache[n_waaskirfcache].b3 = b3;
+		waaskirfcache[n_waaskirfcache].b4 = b4;
+		waaskirfcache[n_waaskirfcache].b5 = b5;
+		waaskirfcache[n_waaskirfcache++].n = strdup(n);
+		n_waaskirfcache = n_waaskirfcache % N_MEMO;  /* Unlikely */
+
+	}
 
 	s2 = pow(s/1e10, 2.0);  /* s2 is s squared in Angstroms^-2 */
 	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;
 }
 
@@ -480,7 +526,7 @@ void get_reflections_cached(struct molecule *mol, double en)
 
 	mol->reflections = new_list_sfac();
 	r = fread(mol->reflections, sizeof(double complex), IDIM*IDIM*IDIM, fh);
-	if ( r <  IDIM*IDIM*IDIM ) {
+	if ( r < IDIM*IDIM*IDIM ) {
 		printf("Found cache file (%s), but failed to read.\n", s);
 		goto calc;
 	}
diff --git a/src/utils.h b/src/utils.h
index 39988ae1..bb0dd77e 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -126,7 +126,7 @@ static inline double distance3d(double x1, double y1, double z1,
 /* -------------------- Indexed lists for specified types ------------------- */
 
 /* Maxmimum index to hold values up to (can be increased if necessary) */
-#define INDMAX 40
+#define INDMAX 70
 
 /* Array size */
 #define IDIM (INDMAX*2 +1)