diff options
author | Thomas White <taw@bitwiz.org.uk> | 2009-11-26 23:55:47 +0100 |
---|---|---|
committer | Thomas White <taw@bitwiz.org.uk> | 2009-11-26 23:55:47 +0100 |
commit | 9e3ebf0fd97a3a2d47093a8da9fa2d55f57394a9 (patch) | |
tree | aa860c288f8579aaae12f49b40e6e1e9d492e240 | |
parent | 6570641ec8250465e9a91bf4ec363a61b5f06ba3 (diff) | |
parent | f5b2f73b8d4775f1853c2a5002974a0aa1605b5d (diff) |
Merge branch 'master' of ssh://git-weiss@git.bitwiz.org.uk/srv/git-private/crystfel
-rw-r--r-- | src/detector.c | 23 | ||||
-rw-r--r-- | src/pattern_sim.c | 4 | ||||
-rw-r--r-- | src/sfac.c | 130 | ||||
-rw-r--r-- | src/utils.h | 2 |
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); @@ -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) |