More new geometry fixes

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];