aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/calibrate_detector.c2
-rw-r--r--src/cubeit.c2
-rw-r--r--src/detector.c2
-rw-r--r--src/detector.h2
-rw-r--r--src/diffraction.c29
-rw-r--r--src/index.c2
-rw-r--r--src/peaks.c57
-rw-r--r--src/powder_plot.c2
8 files changed, 47 insertions, 51 deletions
diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c
index 20503388..ab8b3ecb 100644
--- a/src/calibrate_detector.c
+++ b/src/calibrate_detector.c
@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
int intensity;
struct rvec r;
- r = get_q(&image, x, y, 1, NULL, 1.0/image.lambda);
+ r = get_q(&image, x, y, NULL, 1.0/image.lambda);
q = modulus(r.u, r.v, r.w);
intensity = image.data[x + image.width*y];
diff --git a/src/cubeit.c b/src/cubeit.c
index d5a0a539..e412c2c1 100644
--- a/src/cubeit.c
+++ b/src/cubeit.c
@@ -247,7 +247,7 @@ static void sum_image(void *pg, int cookie)
struct rvec q;
signed int ha, ka, la;
- q = get_q(&image, x, y, 1, NULL, 1.0/image.lambda);
+ q = get_q(&image, x, y, NULL, 1.0/image.lambda);
hd = q.u * ax + q.v * ay + q.w * az;
kd = q.u * bx + q.v * by + q.w * bz;
diff --git a/src/detector.c b/src/detector.c
index 47a3d51a..eb5d1083 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -63,7 +63,7 @@ static int dir_conv(const char *a, signed int *sx, signed int *sy)
struct rvec get_q(struct image *image, double fs, double ss,
- unsigned int sampling, float *ttp, float k)
+ double *ttp, double k)
{
struct rvec q;
double twotheta, r, az;
diff --git a/src/detector.h b/src/detector.h
index 6821ad27..028bd589 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -58,7 +58,7 @@ struct detector
};
extern struct rvec get_q(struct image *image, double xs, double ys,
- unsigned int sampling, float *ttp, float k);
+ double *ttp, double k);
extern double get_tt(struct image *image, double xs, double ys);
diff --git a/src/diffraction.c b/src/diffraction.c
index 5b949594..88be4ffe 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -16,6 +16,7 @@
#include <string.h>
#include <complex.h>
#include <assert.h>
+#include <fenv.h>
#include "image.h"
#include "utils.h"
@@ -319,13 +320,12 @@ static double molecule_factor(const double *intensities, const double *phases,
}
-
void get_diffraction(struct image *image, int na, int nb, int nc,
const double *intensities, const double *phases,
const unsigned char *flags, UnitCell *cell,
GradientMethod m, const char *sym)
{
- unsigned int xs, ys;
+ unsigned int fs, ss;
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
@@ -343,31 +343,36 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0));
bwstep = (khigh-klow) / BWSAMPLING;
- for ( xs=0; xs<image->width*SAMPLING; xs++ ) {
- for ( ys=0; ys<image->height*SAMPLING; ys++ ) {
+ for ( fs=0; fs<image->width*SAMPLING; fs++ ) {
+ for ( ss=0; ss<image->height*SAMPLING; ss++ ) {
struct rvec q;
- float twotheta;
+ double twotheta;
double sw = 1.0/(SAMPLING*SAMPLING); /* Sample weight */
- const unsigned int x = xs / SAMPLING;
- const unsigned int y = ys / SAMPLING; /* Integer part only */
+ const double dfs = fs / SAMPLING;
+ const double dss = ss / SAMPLING;
int kstep;
for ( kstep=0; kstep<BWSAMPLING; kstep++ ) {
- float k;
+ double k;
double kw = 1.0/BWSAMPLING;
double intensity;
double f_lattice, I_lattice;
double I_molecule;
+ int sfs, sss;
/* Calculate k this time round */
k = klow + kstep * bwstep;
- q = get_q(image, xs, ys, SAMPLING, &twotheta, k);
- image->twotheta[x + image->width*y] = twotheta;
+ q = get_q(image, dfs, dss, &twotheta, k);
+
+ /* Determine which pixel this subsample belongs to */
+ fesetround(1); /* Round to nearest */
+ sfs = round(dfs); sss = round(dss);
+ image->twotheta[sfs + image->width*sss] = twotheta;
f_lattice = lattice_factor(q, ax, ay, az,
bx, by, bz,
@@ -387,12 +392,12 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
I_lattice = pow(f_lattice, 2.0);
intensity = sw * kw * I_lattice * I_molecule;
- image->data[x + image->width*y] += intensity;
+ image->data[sfs + image->width*sss] += intensity;
}
}
- progress_bar(xs, SAMPLING*image->width-1, "Calculating diffraction");
+ progress_bar(fs, SAMPLING*image->width-1, "Calculating diffraction");
}
}
diff --git a/src/index.c b/src/index.c
index 5675eb0f..15875295 100644
--- a/src/index.c
+++ b/src/index.c
@@ -121,7 +121,7 @@ void map_all_peaks(struct image *image)
f = image_get_feature(image->features, i);
if ( f == NULL ) continue;
- r = get_q(image, f->x, f->y, 1, NULL, 1.0/image->lambda);
+ r = get_q(image, f->x, f->y, NULL, 1.0/image->lambda);
f->rx = r.u; f->ry = r.v; f->rz = r.w;
}
diff --git a/src/peaks.c b/src/peaks.c
index 640188e7..24bdca5e 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -463,7 +463,7 @@ void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex)
f = image_get_feature(image->features, i);
if ( f == NULL ) continue;
- r = get_q(image, f->x, f->y, 1, NULL, 1.0/image->lambda);
+ r = get_q(image, f->x, f->y, NULL, 1.0/image->lambda);
q = modulus(r.u, r.v, r.w);
fprintf(ofh, "%8.3f %8.3f %8.3f %12.3f\n",
@@ -480,7 +480,7 @@ void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex)
RefList *find_projected_peaks(struct image *image, UnitCell *cell,
int circular_domain, double domain_r)
{
- int x, y;
+ int fs, ss;
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
@@ -499,8 +499,8 @@ RefList *find_projected_peaks(struct image *image, UnitCell *cell,
cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
fesetround(1); /* Round towards nearest */
- for ( x=0; x<image->width; x++ ) {
- for ( y=0; y<image->height; y++ ) {
+ for ( fs=0; fs<image->width; fs++ ) {
+ for ( ss=0; ss<image->height; ss++ ) {
double hd, kd, ld; /* Indices with decimal places */
double dh, dk, dl; /* Distances in h,k,l directions */
@@ -510,7 +510,7 @@ RefList *find_projected_peaks(struct image *image, UnitCell *cell,
Reflection *refl;
double cur_dist;
- q = get_q(image, x, y, 1, NULL, 1.0/image->lambda);
+ q = get_q(image, fs, ss, NULL, 1.0/image->lambda);
hd = q.u * ax + q.v * ay + q.w * az;
kd = q.u * bx + q.v * by + q.w * bz;
@@ -539,12 +539,12 @@ RefList *find_projected_peaks(struct image *image, UnitCell *cell,
if ( refl != NULL ) {
cur_dist = get_excitation_error(refl);
if ( dist < cur_dist ) {
- set_detector_pos(refl, dist, x, y);
+ set_detector_pos(refl, dist, fs, ss);
}
} else {
Reflection *new;
new = add_refl(reflections, h, k, l);
- set_detector_pos(new, dist, x, y);
+ set_detector_pos(new, dist, fs, ss);
n_reflections++;
}
@@ -592,7 +592,7 @@ int peak_sanity_check(struct image *image, UnitCell *cell,
n_feat++;
/* Get closest hkl */
- q = get_q(image, f->x, f->y, 1, NULL, 1.0/image->lambda);
+ q = get_q(image, f->x, f->y, NULL, 1.0/image->lambda);
hd = q.u * ax + q.v * ay + q.w * az;
kd = q.u * bx + q.v * by + q.w * bz;
@@ -784,7 +784,7 @@ void output_pixels(struct image *image, UnitCell *cell,
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
- int x, y;
+ int fs, ss;
double aslen, bslen, cslen;
double *intensities;
double *xmom;
@@ -812,8 +812,8 @@ void output_pixels(struct image *image, UnitCell *cell,
&cx, &cy, &cz);
/* For each pixel */
fesetround(1); /* Round towards nearest */
- for ( x=0; x<image->width; x++ ) {
- for ( y=0; y<image->height; y++ ) {
+ for ( fs=0; fs<image->width; fs++ ) {
+ for ( ss=0; ss<image->height; ss++ ) {
double hd, kd, ld; /* Indices with decimal places */
double dh, dk, dl; /* Distances in h,k,l directions */
@@ -821,12 +821,13 @@ void output_pixels(struct image *image, UnitCell *cell,
struct rvec q;
double dist;
struct panel *p;
+ double twotheta;
- p = find_panel(image->det, x, y);
+ p = find_panel(image->det, fs, ss);
if ( p == NULL ) continue;
if ( p->no_index ) continue;
- q = get_q(image, x, y, 1, NULL, 1.0/image->lambda);
+ q = get_q(image, fs, ss, &twotheta, 1.0/image->lambda);
hd = q.u * ax + q.v * ay + q.w * az;
kd = q.u * bx + q.v * by + q.w * bz;
@@ -853,37 +854,29 @@ void output_pixels(struct image *image, UnitCell *cell,
if ( dist < domain_r ) {
double val;
- struct panel *p;
double pix_area, Lsq, proj_area, dsq, sa;
double phi, pa, pb, pol;
- float tt = 0.0;
double xs, ys, rx, ry;
/* Veto if we want to integrate a bad region */
if ( image->flags != NULL ) {
int flags;
- flags = image->flags[x+image->width*y];
+ flags = image->flags[fs+image->width*ss];
if ( !(flags & 0x01) ) continue;
}
- val = image->data[x+image->width*y];
-
- p = find_panel(image->det, x, y);
- if ( p == NULL ) continue;
-
- if ( p->no_index ) continue;
+ val = image->data[fs+image->width*ss];
/* Area of one pixel */
pix_area = pow(1.0/p->res, 2.0);
Lsq = pow(p->clen, 2.0);
/* Area of pixel as seen from crystal */
- tt = get_tt(image, x, y);
- proj_area = pix_area * cos(tt);
+ proj_area = pix_area * cos(twotheta);
/* Calculate distance from crystal to pixel */
- xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx;
- ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy;
+ xs = (fs-p->min_fs)*p->fsx + (ss-p->min_ss)*p->ssx;
+ ys = (fs-p->min_fs)*p->fsy + (ss-p->min_ss)*p->ssy;
rx = (xs + p->cnx) / p->res;
ry = (ys + p->cny) / p->res;
dsq = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
@@ -896,11 +889,9 @@ void output_pixels(struct image *image, UnitCell *cell,
if ( do_polar ) {
- tt = get_tt(image, x, y);
-
- phi = atan2(y, x);
- pa = pow(sin(phi)*sin(tt), 2.0);
- pb = pow(cos(tt), 2.0);
+ phi = atan2(ry, rx);
+ pa = pow(sin(phi)*sin(twotheta), 2.0);
+ pb = pow(cos(twotheta), 2.0);
pol = 1.0 - 2.0*POL*(1-pa) + POL*(1.0+pb);
val /= pol;
@@ -910,8 +901,8 @@ void output_pixels(struct image *image, UnitCell *cell,
/* Add value to sum */
integrate_intensity(intensities, h, k, l, val);
- integrate_intensity(xmom, h, k, l, val*x);
- integrate_intensity(ymom, h, k, l, val*y);
+ integrate_intensity(xmom, h, k, l, val*fs);
+ integrate_intensity(ymom, h, k, l, val*ss);
if ( !find_item(obs, h, k, l) ) {
add_item(obs, h, k, l);
diff --git a/src/powder_plot.c b/src/powder_plot.c
index b079a194..b1beb414 100644
--- a/src/powder_plot.c
+++ b/src/powder_plot.c
@@ -111,7 +111,7 @@ int main(int argc, char *argv[])
int intensity;
struct rvec r;
- r = get_q(&image, x, y, 1, NULL, 1.0/image.lambda);
+ r = get_q(&image, x, y, NULL, 1.0/image.lambda);
q = modulus(r.u, r.v, r.w);
intensity = image.data[x + image.width*y];
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 03:16:16 +0000