aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/detector.c23
-rw-r--r--src/pattern_sim.c4
-rw-r--r--src/sfac.c130
-rw-r--r--src/utils.h2
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)
generated by cgit v1.2.3 (git 2.25.1) at 2024-08-18 18:29:05 +0000