From 26bf7df5817d1df238698dc8f63e4543e746c7c1 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Fri, 13 Nov 2009 17:13:11 +0100
Subject: Solid angle stuff

---
 src/detector.c | 12 +++++++++---
 src/ewald.c    |  6 ++++++
 src/image.h    |  2 ++
 src/main.c     |  1 +
 4 files changed, 18 insertions(+), 3 deletions(-)

(limited to 'src')

diff --git a/src/detector.c b/src/detector.c
index 2b0c12d4..339ff626 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -26,26 +26,32 @@ void record_image(struct image *image)
 {
 	int x, y;
 	double ph_per_e;
+	double twotheta_max, np, sa_per_pixel;
 	
 	/* How many photons are scattered per electron? */
 	ph_per_e = PULSE_ENERGY_DENSITY * pow(THOMSON_LENGTH, 2.0)
 	          / image->xray_energy;
-	
 	printf("%e photons are scattered per electron\n", ph_per_e);
 	
+	twotheta_max = image->twotheta[0];
+	
+	np = sqrt(pow(image->x_centre, 2.0) + pow(image->y_centre, 2.0));
+	sa_per_pixel = pow(2.0 * twotheta_max / np, 2.0);
+	printf("sa per pixel=%e\n", sa_per_pixel);
+	
 	image->data = malloc(image->width * image->height * sizeof(uint16_t));
 	
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
 	
-		uint16_t counts;
+		double counts;
 		double val, intensity;
 		double sa;
 		
 		val = image->sfacs[x + image->width*y];
 		
 		/* What solid angle is subtended by this pixel? */
-		sa = 1.0;
+		sa = sa_per_pixel * cos(image->twotheta[x + image->width*y]);
 		
 		intensity = pow(val, 2.0);
 		counts = intensity * ph_per_e * sa;
diff --git a/src/ewald.c b/src/ewald.c
index d543d98a..5be65644 100644
--- a/src/ewald.c
+++ b/src/ewald.c
@@ -30,6 +30,9 @@ void get_ewald(struct image *image)
 	image->qvecs = malloc(image->width * image->height
 	                       * sizeof(struct threevec));
 	
+	image->twotheta = malloc(image->width * image->height
+	                       * sizeof(double));
+	
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
 	
@@ -50,9 +53,12 @@ void get_ewald(struct image *image)
 		qy = k * sin(twothetay);
 		qz = k - k * cos(twotheta);
 		
+		/* FIXME: Rotate vector here */
+		
 		image->qvecs[x + image->width*y].u = qx;
 		image->qvecs[x + image->width*y].v = qy;
 		image->qvecs[x + image->width*y].w = qz;
+		image->twotheta[x + image->width*y] = twotheta;
 	
 	}
 	}
diff --git a/src/image.h b/src/image.h
index b675ab41..dd6f7c6b 100644
--- a/src/image.h
+++ b/src/image.h
@@ -62,6 +62,8 @@ struct image {
 	uint16_t		*data;
 	double			*sfacs;
 	struct threevec		*qvecs;
+	double			*twotheta;
+	
 
 	/* Radians.  Defines where the pattern lies in reciprocal space */
 	double			tilt;
diff --git a/src/main.c b/src/main.c
index 27ac1f9c..86cf93d6 100644
--- a/src/main.c
+++ b/src/main.c
@@ -80,6 +80,7 @@ int main(int argc, char *argv[])
 	image.qvecs = NULL;
 	image.sfacs = NULL;
 	image.data = NULL;
+	image.twotheta = NULL;
 	
 	/* Splurge a few useful numbers */
 	printf("Wavelength is %f nm\n", image.lambda/1.0e-9);
-- 
cgit v1.2.3