/* * detector.c * * Detector properties * * (c) 2007-2009 Thomas White * * pattern_sim - Simulate diffraction patterns from small crystals * */ #include #include #include #include #include "image.h" #include "utils.h" /* Pulse energy density in J/m^2 */ #define PULSE_ENERGY_DENSITY (30.0e7) /* Detector's quantum efficiency */ #define DQE (0.9) /* Detector's saturation value */ #define SATURATION (60000) /* Bleed excess intensity into neighbouring pixels */ static void bloom_values(double *tmp, int x, int y, int width, int height, double val) { double overflow; overflow = val - SATURATION; /* Intensity which bleeds off the edge of the detector is lost */ if ( x > 0 ) { tmp[x-1 + width*y] += overflow / 6.0; if ( y > 0 ) { tmp[x-1 + width*(y-1)] += overflow / 12.0; } if ( y < height-1 ) { tmp[x-1 + width*(y+1)] += overflow / 12.0; } } if ( x < width-1 ) { tmp[x+1 + width*y] += overflow / 6.0; if ( y > 0 ) { tmp[x+1 + width*(y-1)] += overflow / 12.0; } if ( y < height-1 ) { tmp[x+1 + width*(y+1)] += overflow / 12.0; } } if ( y > 0 ) { tmp[x + width*(y-1)] += overflow / 6.0; } if ( y < height-1 ) { tmp[x + width*(y+1)] += overflow / 6.0; } } static uint16_t *bloom(double *hdr_in, int width, int height) { int x, y; uint16_t *data; double *tmp; double *hdr; int did_something; data = malloc(width * height * sizeof(uint16_t)); tmp = malloc(width * height * sizeof(double)); hdr = malloc(width * height * sizeof(double)); memcpy(hdr, hdr_in, width*height*sizeof(double)); /* Apply DQE (once only) */ for ( x=0; xxray_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->hdr = malloc(image->width * image->height * sizeof(double)); for ( x=0; xwidth; x++ ) { for ( y=0; yheight; y++ ) { double counts; double intensity; double sa; double complex val; val = image->sfacs[x + image->width*y]; intensity = (double)(val * conj(val)); /* What solid angle is subtended by this pixel? */ sa = sa_per_pixel * cos(image->twotheta[x + image->width*y]); counts = intensity * ph_per_e * sa; image->hdr[x + image->width*y] = counts; } } image->data = bloom(image->hdr, image->width, image->height); }