From 7601a2956156569461498a76105267129e848fb0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 27 Jun 2019 16:34:45 +0200 Subject: partialator: Use Spectrum API --- libcrystfel/src/geometry.c | 27 ++++++++------------------- src/partialator.c | 29 ++++++++++++++++++----------- src/post-refinement.c | 22 +++++++++++++++++++++- 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; ilambda; - 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 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); } -- cgit v1.2.3