aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/geometry.c27
-rw-r--r--src/partialator.c29
-rw-r--r--src/post-refinement.c22
3 files changed, 47 insertions, 31 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index ffb17133..4f9515c5 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -522,7 +522,7 @@ static void set_random_partialities(Crystal *cryst)
static double do_integral(double q2, double zl, double R,
- double lambda, double sig, char *verbose)
+ Spectrum *spectrum, char *verbose)
{
int i;
double kmin, kmax, kstart, kfinis;
@@ -530,7 +530,6 @@ static double do_integral(double q2, double zl, double R,
double total = 0.0;
double k0, k1;
const int SAMPLES = 50; /* Number of samples for integration */
- const double N = 1.5; /* Pointiness of spectrum */
FILE *fh = NULL;
assert(R*R < q2);
@@ -545,8 +544,7 @@ static double do_integral(double q2, double zl, double R,
if ( k1 < 0.0 ) k1 = +INFINITY;
/* Range over which E is significantly different from zero */
- kmin = 1.0 / (lambda + 3.0*sig);
- kmax = 1.0 / (lambda - 3.0*sig);
+ spectrum_get_range(spectrum, &kmin, &kmax);
/* Calculate range over which E*P is different from zero */
if ( k0 < k1 ) {
@@ -592,9 +590,6 @@ static double do_integral(double q2, double zl, double R,
snprintf(fn, 63, "partial%s.graph", verbose);
fh = fopen(fn, "w");
fprintf(fh, " n p wavelength E P\n");
- STATUS("Nominal k = %e m^-1\n", 1.0/lambda);
- STATUS(" (wavelength %e m)\n", lambda);
- STATUS("Bandwidth %e m\n", sig);
STATUS("k1/2 = %e m^-1\n", -q2/(2.0*zl));
STATUS(" (wavelength %e m)\n", 1.0/(-q2/(2.0*zl)));
STATUS("Reflection k goes from %e to %e m^-1\n", k1, k0);
@@ -607,7 +602,7 @@ static double do_integral(double q2, double zl, double R,
for ( i=0; i<SAMPLES; i++ ) {
- double p, kp, lrel;
+ double p, kp;
double E, P;
kp = kstart + i*inc;
@@ -615,9 +610,7 @@ static double do_integral(double q2, double zl, double R,
p = pref + 0.5 - kp/(2.0*R);
/* Spectral energy term */
- lrel = fabs(1.0/kp - lambda);
- E = exp(-0.5 * pow(lrel / sig, N));
- E /= sqrt(2.0 * M_PI * sig);
+ E = spectrum_get_density_at_k(spectrum, kp);
/* RLP profile term */
P = 4.0*p * (1.0 - p);
@@ -633,9 +626,6 @@ static double do_integral(double q2, double zl, double R,
if ( isnan(total) ) {
ERROR("NaN total!\n");
- STATUS("Nominal k = %e m^-1\n", 1.0/lambda);
- STATUS(" (wavelength %e m)\n", lambda);
- STATUS("Bandwidth %e m\n", sig);
STATUS("k1/2 = %e m^-1\n", -q2/(2.0*zl));
STATUS(" (wavelength %e m)\n", 1.0/(-q2/(2.0*zl)));
STATUS("Reflection k goes from %e to %e m^-1\n", k1, k0);
@@ -655,7 +645,7 @@ static void ginn_spectrum_partialities(Crystal *cryst)
RefList *list;
Reflection *refl;
RefListIterator *iter;
- double r0, m, lambda, sig;
+ double r0, m;
struct image *image;
UnitCell *cell;
double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
@@ -683,8 +673,6 @@ static void ginn_spectrum_partialities(Crystal *cryst)
r0 = fabs(crystal_get_profile_radius(cryst));
m = crystal_get_mosaicity(cryst);
- lambda = image->lambda;
- sig = image->bw * lambda / 2.0;
for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
refl != NULL;
@@ -710,8 +698,9 @@ static void ginn_spectrum_partialities(Crystal *cryst)
//snprintf(tmp, 255, "-%i,%i,%i", h, k, l);
char *tmp = NULL;
- total = do_integral(q2, zl, R, lambda, sig, tmp);
- norm = do_integral(q2, -0.5*q2*lambda, R, lambda, sig, NULL);
+ total = do_integral(q2, zl, R, image->spectrum, tmp);
+ norm = do_integral(q2, -0.5*q2*image->lambda, R,
+ image->spectrum, NULL);
set_partiality(refl, total/norm);
set_lorentz(refl, 1.0);
diff --git a/src/partialator.c b/src/partialator.c
index 92fcfbe4..de7eeb4b 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -570,8 +570,11 @@ static void skip_to_end(FILE *fh)
}
-static int set_initial_params(Crystal *cr, FILE *fh)
+static int set_initial_params(Crystal *cr, FILE *fh, double force_bandwidth,
+ double force_radius, double force_lambda)
{
+ struct image *image = crystal_get_image(cr);
+
if ( fh != NULL ) {
int err;
@@ -596,6 +599,17 @@ static int set_initial_params(Crystal *cr, FILE *fh)
}
+ if ( force_bandwidth > 0.0 ) {
+ image->bw = force_bandwidth;
+ }
+ if ( force_radius > 0.0 ) {
+ crystal_set_profile_radius(cr, force_radius);
+ }
+ if ( force_lambda > 0.0 ) {
+ image->lambda = force_lambda;
+ }
+ image->spectrum = spectrum_generate_gaussian(image->lambda, image->bw);
+
return 0;
}
@@ -1420,7 +1434,9 @@ int main(int argc, char *argv[])
crystal_set_user_flag(cr, PRFLAG_OK);
reflist_free(cr_refl);
- if ( set_initial_params(cr, sparams_fh) ) {
+ if ( set_initial_params(cr, sparams_fh, force_bandwidth,
+ force_radius, force_lambda) )
+ {
ERROR("Failed to set initial parameters\n");
return 1;
}
@@ -1451,15 +1467,6 @@ int main(int argc, char *argv[])
STATUS("Initial partiality calculation...\n");
for ( i=0; i<n_crystals; i++ ) {
Crystal *cr = crystals[i];
- if ( force_bandwidth > 0.0 ) {
- crystal_get_image(cr)->bw = force_bandwidth;
- }
- if ( force_radius > 0.0 ) {
- crystal_set_profile_radius(cr, force_radius);
- }
- if ( force_lambda > 0.0 ) {
- crystal_get_image(cr)->lambda = force_lambda;
- }
update_predictions(cr);
calculate_partialities(cr, pmodel);
}
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 59d44425..6db1a48a 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -153,14 +153,22 @@ static void apply_parameters(const Crystal *cr_orig, Crystal *cr_tgt,
struct rf_alteration alter)
{
double R, lambda;
+ struct gaussian g;
+ struct image *image;
+ image = crystal_get_image(cr_tgt);
R = crystal_get_profile_radius(cr_orig);
lambda = crystal_get_image_const(cr_orig)->lambda;
rotate_cell_xy(crystal_get_cell_const(cr_orig), crystal_get_cell(cr_tgt),
alter.rot_x, alter.rot_y);
crystal_set_profile_radius(cr_tgt, R+alter.delta_R);
- crystal_get_image(cr_tgt)->lambda = lambda+alter.delta_wave;
+ image->lambda = lambda+alter.delta_wave;
+
+ g.kcen = 1.0/image->lambda;
+ g.sigma = image->bw/image->lambda;
+ g.height = 1;
+ spectrum_set_gaussians(image->spectrum, &g, 1);
}
@@ -442,6 +450,7 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
struct rf_priv priv;
RefList *list;
UnitCell *cell;
+ Spectrum *spectrum;
priv.cr = cr;
priv.full = full;
@@ -449,6 +458,8 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
priv.scaleflags = scaleflags;
priv.cr_tgt = crystal_copy(cr);
priv.image_tgt = *crystal_get_image(cr);
+ spectrum = spectrum_new();
+ priv.image_tgt.spectrum = spectrum;
crystal_set_image(priv.cr_tgt, &priv.image_tgt);
list = copy_reflist(crystal_get_reflections(cr));
crystal_set_reflections(priv.cr_tgt, list);
@@ -490,6 +501,7 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
reflist_free(crystal_get_reflections(priv.cr_tgt));
crystal_free(priv.cr_tgt);
+ spectrum_free(spectrum);
}
@@ -502,6 +514,7 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
struct rf_priv priv;
RefList *list;
UnitCell *cell;
+ Spectrum *spectrum;
priv.cr = cr;
priv.full = full;
@@ -509,6 +522,8 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
priv.scaleflags = scaleflags;
priv.cr_tgt = crystal_copy(cr);
priv.image_tgt = *crystal_get_image(cr);
+ spectrum = spectrum_new();
+ priv.image_tgt.spectrum = spectrum;
crystal_set_image(priv.cr_tgt, &priv.image_tgt);
list = copy_reflist(crystal_get_reflections(cr));
crystal_set_reflections(priv.cr_tgt, list);
@@ -550,6 +565,7 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
reflist_free(crystal_get_reflections(priv.cr_tgt));
crystal_free(priv.cr_tgt);
+ spectrum_free(spectrum);
}
@@ -635,6 +651,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
RefList *list;
FILE *fh = NULL;
UnitCell *cell;
+ Spectrum *spectrum;
try_reindex(cr, full, sym, amb, scaleflags);
@@ -651,6 +668,8 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
priv.scaleflags = scaleflags;
priv.cr_tgt = crystal_copy(cr);
priv.image_tgt = *crystal_get_image(cr);
+ spectrum = spectrum_new();
+ priv.image_tgt.spectrum = spectrum;
crystal_set_image(priv.cr_tgt, &priv.image_tgt);
list = copy_reflist(crystal_get_reflections(cr));
crystal_set_reflections(priv.cr_tgt, list);
@@ -718,6 +737,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
reflist_free(crystal_get_reflections(priv.cr_tgt));
crystal_free(priv.cr_tgt);
+ spectrum_free(spectrum);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-24 12:31:00 +0000