aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2009-11-13 12:10:12 +0100
committerThomas White <taw@physics.org>2009-11-13 12:10:12 +0100
commitef6a971cf432321ceb057b8c355c8e6814d5aff6 (patch)
tree1d8d7373a1b7c1591cf28fcf61230122e19c1ca4
parent857cc8c991a5ee4e717d01822da271ae34f5adaa (diff)
Work in progress on photon correctness
-rw-r--r--src/detector.c4
-rw-r--r--src/ewald.c18
-rw-r--r--src/image.h5
-rw-r--r--src/main.c6
-rw-r--r--src/utils.h26
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 */