aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2020-06-03 17:16:59 +0200
committerThomas White <taw@physics.org>2020-07-29 18:53:33 +0200
commit6e5d500a2786131c7304b6c001142e524791bc12 (patch)
tree0225c8f3dea1545332ef9294693ffc6b3263097b /src
parent2037afa4693d4032606d98eb67522fabb58a42ea (diff)
Remove "struct detector" completely, part I
record_image has been moved to pattern_sim.c
Diffstat (limited to 'src')
-rw-r--r--src/pattern_sim.c126
1 files changed, 126 insertions, 0 deletions
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 1fce0f60..9a144a36 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -116,6 +116,132 @@ static void show_help(const char *s)
}
+static void record_panel(struct panel *p, float *dp, int do_poisson,
+ gsl_rng *rng, double ph_per_e, double background,
+ double lambda,
+ int *n_neg1, int *n_inf1, int *n_nan1,
+ int *n_neg2, int *n_inf2, int *n_nan2,
+ double *max_tt)
+{
+ int fs, ss;
+
+ for ( ss=0; ss<p->h; ss++ ) {
+ for ( fs=0; fs<p->w; fs++ ) {
+
+ double counts;
+ double cf;
+ double intensity, sa;
+ double pix_area, Lsq;
+ double xs, ys, rx, ry;
+ double dsq, proj_area;
+ float dval;
+ double twotheta;
+
+ intensity = (double)dp[fs + p->w*ss];
+ if ( isinf(intensity) ) (*n_inf1)++;
+ if ( intensity < 0.0 ) (*n_neg1)++;
+ if ( isnan(intensity) ) (*n_nan1)++;
+
+ /* Area of one pixel */
+ pix_area = pow(1.0/p->res, 2.0);
+ Lsq = pow(p->clen, 2.0);
+
+ /* Calculate distance from crystal to pixel */
+ xs = fs*p->fsx + ss*p->ssx;
+ ys = ss*p->fsy + ss*p->ssy;
+ rx = (xs + p->cnx) / p->res;
+ ry = (ys + p->cny) / p->res;
+ dsq = pow(rx, 2.0) + pow(ry, 2.0);
+ twotheta = atan2(sqrt(dsq), p->clen);
+
+ /* Area of pixel as seen from crystal (approximate) */
+ proj_area = pix_area * cos(twotheta);
+
+ /* Projected area of pixel divided by distance squared */
+ sa = proj_area / (dsq + Lsq);
+
+ if ( do_poisson ) {
+ counts = poisson_noise(rng, intensity * ph_per_e * sa);
+ } else {
+ cf = intensity * ph_per_e * sa;
+ counts = cf;
+ }
+
+ /* Number of photons in pixel */
+ dval = counts + poisson_noise(rng, background);
+
+ /* Convert to ADU */
+ dval *= p->adu_per_photon;
+
+ /* Saturation */
+ if ( dval > p->max_adu ) dval = p->max_adu;
+
+ dp[fs + p->w*ss] = dval;
+
+ /* Sanity checks */
+ if ( isinf(dp[fs + p->w*ss]) ) n_inf2++;
+ if ( isnan(dp[fs + p->w*ss]) ) n_nan2++;
+ if ( dp[fs + p->w*ss] < 0.0 ) n_neg2++;
+
+ if ( twotheta > *max_tt ) *max_tt = twotheta;
+
+ }
+ }
+}
+
+
+void record_image(struct image *image, int do_poisson, double background,
+ gsl_rng *rng, double beam_radius, double nphotons)
+{
+ double total_energy, energy_density;
+ double ph_per_e;
+ double area;
+ double max_tt = 0.0;
+ int n_inf1 = 0;
+ int n_neg1 = 0;
+ int n_nan1 = 0;
+ int n_inf2 = 0;
+ int n_neg2 = 0;
+ int n_nan2 = 0;
+ int pn;
+
+ /* How many photons are scattered per electron? */
+ area = M_PI*pow(beam_radius, 2.0);
+ total_energy = nphotons * ph_lambda_to_en(image->lambda);
+ energy_density = total_energy / area;
+ ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0);
+ STATUS("Fluence = %8.2e photons, "
+ "Energy density = %5.3f kJ/cm^2, "
+ "Total energy = %5.3f microJ\n",
+ nphotons, energy_density/1e7, total_energy*1e6);
+
+ fill_in_adu(image);
+
+ for ( pn=0; pn<image->det->n_panels; pn++ ) {
+
+ record_panel(&image->det->panels[pn], image->dp[pn],
+ do_poisson, rng, ph_per_e, background,
+ image->lambda,
+ &n_neg1, &n_inf1, &n_nan1,
+ &n_neg2, &n_inf2, &n_nan2, &max_tt);
+ }
+
+ STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n",
+ rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9);
+
+ STATUS("Halve the d values to get the voxel size for a synthesis.\n");
+
+ if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) {
+ ERROR("WARNING: The raw calculation produced %i negative"
+ " values, %i infinities and %i NaNs.\n",
+ n_neg1, n_inf1, n_nan1);
+ ERROR("WARNING: After processing, there were %i negative"
+ " values, %i infinities and %i NaNs.\n",
+ n_neg2, n_inf2, n_nan2);
+ }
+}
+
+
static double *intensities_from_list(RefList *list, SymOpList *sym)
{
Reflection *refl;