From b6af09d78d49d552a05b81bd20a7854a180373f5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 23 Aug 2019 17:05:06 +0200 Subject: Formatting/fussiness --- libcrystfel/src/geometry.c | 137 ++++++++++++++++++++++++++++++--------------- 1 file changed, 92 insertions(+), 45 deletions(-) (limited to 'libcrystfel/src/geometry.c') diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 366464fa..4d25b89f 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -238,14 +238,16 @@ static double random_partiality(signed int h, signed int k, signed int l, return p; } + static inline double safe_khalf(const double xl, const double yl, const double zl) { - if (zl>0) return 0.0/0.0; + if ( zl > 0.0 ) return 0.0/0.0; return -(xl*xl+yl*yl+zl*zl) / (2.0*zl); } + static Reflection *check_reflection(struct image *image, Crystal *cryst, signed int h, signed int k, signed int l, double xl, double yl, double zl, @@ -253,73 +255,118 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, { Reflection *refl; double R; - double knom=0, khalf=0, kpred=0; + int i, n; + double knom, khalf; double dcs, exerr; - double partiality = 0, mean_kpred=0, M2_kpred=0; - double sumw_k = 0, mean_k = 0, M2_k=0; - struct gaussian g; - /* this arbitrary value is there to mimic previous behavoir */ + double partiality, mean_kpred, M2_kpred; + double sumw_k, mean_k, M2_k; + + /* This arbitrary value is there to mimic previous behaviour */ const double min_partiality = exp(-0.5*1.7*1.7); /* Don't predict 000 */ if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL; R = fabs(crystal_get_profile_radius(cryst)); - assert(spectrum_get_num_gaussians(image->spectrum)); - for ( int i = 0 ; i!=spectrum_get_num_gaussians(image->spectrum) ; ++i ) { + n = spectrum_get_num_gaussians(image->spectrum); + assert(n > 0); + + partiality = 0.0; + mean_kpred = 0.0; + M2_kpred = 0.0; + sumw_k = 0.0; + mean_k = 0.0; + M2_k = 0.0; + for ( i=0; ispectrum, i); - double exerr2 = 0; - /* project lattice point onto ewald sphere */ - double x = xl, y = yl, z = zl + g.kcen; - double norm = 1.0/sqrt(x*x+y*y+z*z); - x*=norm; y*=norm; z*=norm; - z-=1; - /* width of ewald sphere in the direction of the projection */ - const double sigma_proj = sqrt(x*x+y*y+z*z)*g.sigma; - mean_variance(g.kcen,g.area,&sumw_k,&mean_k,&M2_k); - M2_k+=g.area*g.sigma*g.sigma; - const double w0 = 1.0/(R*R); - const double w1 = 1.0/(sigma_proj*sigma_proj); - x*=g.kcen;y*=g.kcen;z*=g.kcen; - /* three because the general case fails in extreme cases */ - if ( w0 / w1 <= DBL_MIN ) { /* "Laue" corner case */ + /* Project lattice point onto Ewald sphere */ + x = xl; + y = yl; + z = zl + g.kcen; + norm = 1.0/sqrt(x*x+y*y+z*z); + x *= norm; + y *= norm; + z *= norm; + z -= 1.0; + + /* Width of Ewald sphere in the direction of the projection */ + sigma_proj = sqrt(x*x+y*y+z*z)*g.sigma; + + mean_variance(g.kcen, g.area, &sumw_k, &mean_k, &M2_k); + M2_k += g.area * g.sigma * g.sigma; + w0 = 1.0/(R*R); + w1 = 1.0/(sigma_proj*sigma_proj); + + x *= g.kcen; + y *= g.kcen; + z *= g.kcen; + + /* Three because the general case fails in extreme cases */ + if ( w0 / w1 <= DBL_MIN ) { + + /* 'Laue' corner case */ kpred = g.kcen; - exerr2 = g.kcen - safe_khalf(xl,yl,zl); - exerr2*= exerr2; - } else if ( w1 / w0 <= DBL_MIN ) { /* monochromatic corner case */ + exerr2 = g.kcen - safe_khalf(xl, yl, zl); + exerr2 *= exerr2; + + } else if ( w1 / w0 <= DBL_MIN ) { + + /* 'Monochromatic' corner case */ kpred = safe_khalf(xl,yl,zl); exerr2 = g.kcen - kpred; exerr2*= exerr2; - } else { /* general case */ - /* closest point on ewald sphere */ - const double zlp0 = zl>0?zl:0; /* project zl to 0, bit of a hack... */ + + } else { + + /* General case */ + + /* Closest point on Ewald sphere. + * Project zl to 0, bit of a hack... */ + const double zlp0 = zl>0?zl:0; exerr2 = (x-xl)*(x-xl) + (y-yl)*(y-yl) + (z-zl)*(z-zl); - /* weighted average between projected lattice point and ewald sphere */ + + /* Weighted average between projected lattice point + * and Ewald sphere */ x = ( xl *w0 + x*w1 ) / ( w0 + w1 ); y = ( yl *w0 + y*w1 ) / ( w0 + w1 ); z = ( zlp0*w0 + z*w1 ) / ( w0 + w1 ); kpred = safe_khalf(x,y,z); + + } + sigma2 = R*R + sigma_proj*sigma_proj; + exponent = - 0.5 * exerr2 / sigma2; + if ( exponent > -700.0 ) { + overlap_integral = exp(exponent) * sqrt(2*M_PI*R*R) / sqrt(2*M_PI*sigma2); + } else { + exponent = 0.0; } - const double sigma2 = R*R + sigma_proj*sigma_proj; - const double exponent = - 0.5 * exerr2 / sigma2; - const double overlap_integral = exponent < -700.0 ? 0.0 : - exp( exponent ) - * sqrt( 2 * M_PI * R * R ) - / sqrt( 2 * M_PI * sigma2 ); - mean_variance(kpred,g.area*overlap_integral,&partiality,&mean_kpred,&M2_kpred); + + mean_variance(kpred, g.area*overlap_integral, + &partiality, &mean_kpred, &M2_kpred); + } - partiality *= sqrt( ( R*R + M2_k/sumw_k) / ( R*R ) ); /* reverting the lorentz factor */ + + /* Revert the 'Lorentz' factor */ + partiality *= sqrt( ( R*R + M2_k/sumw_k) / ( R*R ) ); + if ( (updateme == NULL) && ( partiality < min_partiality ) ) return NULL; - kpred = mean_kpred; /* Calculate excitation error */ knom = 1.0/image->lambda; dcs = distance3d(0.0, 0.0, -knom, xl, yl, zl); exerr = 1.0/image->lambda - dcs; /* Loss of precision */ - /* to estimate excitation error from partiality but loosing sign */ - /* that seems to be a problem */ - /* exerr = R * sqrt( - 2 * log( partiality ) ) ; */ + + /* Could also estimate excitation error like this, but it loses the sign + * (which is needed elsewhere): + * exerr = R * sqrt(-2*log(partiality)); */ if ( updateme == NULL ) { refl = reflection_new(h, k, l); @@ -333,7 +380,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, if ( (image->det != NULL) && (updateme != NULL) ) { double fs, ss; - locate_peak_on_panel(xl, yl, zl, kpred, + locate_peak_on_panel(xl, yl, zl, mean_kpred, get_panel(updateme), &fs, &ss); set_detector_pos(refl, fs, ss); @@ -345,7 +392,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, double fs, ss; /* Position on detector */ signed int p; /* Panel number */ - p = locate_peak(xl, yl, zl, kpred, + p = locate_peak(xl, yl, zl, mean_kpred, image->det, &fs, &ss); if ( p == -1 ) { reflection_free(refl); @@ -357,7 +404,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, } khalf = (- xl*xl - yl*yl - zl*zl) / (2.0*zl); - set_kpred(refl, kpred); + set_kpred(refl, mean_kpred); set_khalf(refl, khalf); set_exerr(refl, exerr); set_lorentz(refl, 1.0); -- cgit v1.2.3