Work in progress on photon correctness

5 files changed, 56 insertions, 3 deletions

diff --git a/src/detector.c b/src/detector.c
index ebcbb814..9e4e5649 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -29,11 +29,13 @@ void record_image(struct image *image)
 		uint16_t counts;
 		double val;
 		double intensity;
+		double phac;
 		
 		val = image->sfacs[x + image->width*y];
+		phac = image->phactors[x + image->width*y];
 		
 		intensity = pow(val, 2.0);
-		counts = intensity*16;
+		counts = intensity * phac;
 		
 		image->data[x + image->width*y] = counts;
 	
diff --git a/src/ewald.c b/src/ewald.c
index 9d8a0450..e9f36596 100644
--- a/src/ewald.c
+++ b/src/ewald.c
@@ -19,16 +19,30 @@
 #include "cell.h"
 
 
+/* Pulse energy density in J/m^2 */
+#define PULSE_ENERGY_DENSITY (30.0e7)
+
+
 void get_ewald(struct image *image)
 {
 	int x, y;
 	double k;  /* Wavenumber */
+	double i0fac;
 	
 	k = 1/image->lambda;
 	
+	/* How many photons are scattered per electron? */
+	i0fac = PULSE_ENERGY_DENSITY * pow(THOMSON_LENGTH, 2.0)
+	          / image->xray_energy;
+	
+	printf("%e photons are scattered per electron\n", i0fac);
+	
 	image->qvecs = malloc(image->width * image->height
 	                       * sizeof(struct threevec));
 	
+	image->phactors = malloc(image->width * image->height
+	                          * sizeof(double));
+	
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
 	
@@ -36,6 +50,7 @@ void get_ewald(struct image *image)
 		double twothetax, twothetay, twotheta;
 		double qx, qy, qz;
 		
+		/* Calculate q vectors for Ewald sphere */
 		rx = ((double)x - image->x_centre) / image->resolution;
 		ry = ((double)y - image->y_centre) / image->resolution;
 		r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
@@ -51,6 +66,9 @@ void get_ewald(struct image *image)
 		image->qvecs[x + image->width*y].u = qx;
 		image->qvecs[x + image->width*y].v = qy;
 		image->qvecs[x + image->width*y].w = qz;
+		
+		/* Calculate photon factor ("phactor") */
+		
 	
 	}
 	}
diff --git a/src/image.h b/src/image.h
index 35bd5da9..d5d5bb9c 100644
--- a/src/image.h
+++ b/src/image.h
@@ -62,6 +62,7 @@ struct image {
 	uint16_t		*data;
 	double			*sfacs;
 	struct threevec		*qvecs;
+	double			*phactors;
 
 	/* Radians.  Defines where the pattern lies in reciprocal space */
 	double			tilt;
@@ -78,7 +79,9 @@ struct image {
 	double			resolution;	/* pixels per metre */
 
 	/* Wavelength must always be given */
-	double			lambda;
+	double			lambda;		/* Wavelength in m */
+	double			xray_energy;	/* X-ray energy
+	                                         * in J (per photon) */
 
 	int			width;
 	int			height;
diff --git a/src/main.c b/src/main.c
index c101fc87..27ac1f9c 100644
--- a/src/main.c
+++ b/src/main.c
@@ -75,11 +75,15 @@ int main(int argc, char *argv[])
 	image.y_centre = 255.5;
 	image.camera_len = 0.2;  /* 20 cm */
 	image.resolution = 5120; /* 512 pixels in 10 cm */
-	image.lambda = 0.2e-9;   /* LCLS wavelength */
+	image.xray_energy = eV_to_J(2.0e3); /* 2 keV energy */
+	image.lambda = ph_en_to_lambda(image.xray_energy);  /* Wavelength */
 	image.qvecs = NULL;
 	image.sfacs = NULL;
 	image.data = NULL;
 	
+	/* Splurge a few useful numbers */
+	printf("Wavelength is %f nm\n", image.lambda/1.0e-9);
+	
 	get_diffraction(&image, cell);
 	record_image(&image);
 
diff --git a/src/utils.h b/src/utils.h
index f1b72af8..1de0ea09 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -18,6 +18,20 @@
 
 #include <math.h>
 
+
+/* Electron charge in C */
+#define ELECTRON_CHARGE (1.6021773e-19)
+
+/* Planck's constant (Js) */
+#define PLANCK (6.62606896e-34)
+
+/* Speed of light in vacuo (m/s) */
+#define C_VACUO (299792458)
+
+/* Thomson scattering length (m) */
+#define THOMSON_LENGTH (2.81794e-15)
+
+
 extern unsigned int biggest(signed int a, signed int b);
 extern unsigned int smallest(signed int a, signed int b);
 extern double distance(double x1, double y1, double x2, double y2);
@@ -42,4 +56,16 @@ extern void mapping_rotate(double x, double y, double z,
 
 #define is_odd(a) ((a)%2==1)
 
+/* Photon energy (J) to wavelength (m) */
+#define ph_en_to_lambda(a) ((PLANCK*C_VACUO)/(a))
+
+/* Photon wavelength (m) to energy (J) */
+#define ph_lambda_to_en(a) ((PLANCK*C_VACUO)/(a))
+
+/* eV to Joules */
+#define eV_to_J(a) ((a)*ELECTRON_CHARGE)
+
+/* Joules to eV */
+#define J_to_eV(a) ((a)/ELECTRON_CHARGE)
+
 #endif	/* UTILS_H */