From e44890d0a13815f7fd61d777637c5bd09cb09dbc Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Tue, 21 Sep 2010 15:39:53 +0200
Subject: Add the possibility to use a circular (instead of 'crystallographic')
 integration domain

---
 src/peaks.c | 31 ++++++++++++++++++++++++-------
 1 file changed, 24 insertions(+), 7 deletions(-)

(limited to 'src/peaks.c')

diff --git a/src/peaks.c b/src/peaks.c
index fdf5b961..5cfdbdde 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -451,7 +451,8 @@ void dump_peaks(struct image *image, pthread_mutex_t *mutex)
 }
 
 
-int find_projected_peaks(struct image *image, UnitCell *cell)
+int find_projected_peaks(struct image *image, UnitCell *cell,
+                         int circular_domain, double domain_r)
 {
 	int x, y;
 	double ax, ay, az;
@@ -459,11 +460,15 @@ int find_projected_peaks(struct image *image, UnitCell *cell)
 	double cx, cy, cz;
 	struct reflhit *hits;
 	int n_hits = 0;
+	double alen, blen, clen;
 
 	hits = malloc(sizeof(struct reflhit)*MAX_HITS);
 	if ( hits == NULL ) return 0;
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	alen = modulus(ax, ay, az);
+	blen = modulus(bx, by, bz);
+	clen = modulus(cx, cy, cz);
 
 	fesetround(1);  /* Round towards nearest */
 	for ( x=0; x<image->width; x++ ) {
@@ -490,8 +495,17 @@ int find_projected_peaks(struct image *image, UnitCell *cell)
 		dh = hd - h;
 		dk = kd - k;
 		dl = ld - l;
-		dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0));
-		if ( dist > 0.1 ) continue;
+
+		if ( circular_domain ) {
+			/* Circular integration domain */
+			dist = sqrt(pow(dh*alen, 2.0) + pow(dk*blen, 2.0)
+			                              + pow(dl*clen, 2.0));
+			if ( dist > domain_r ) continue;
+		} else {
+			/* "Crystallographic" integration domain */
+			dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0));
+			if ( dist > domain_r ) continue;
+		}
 
 		for ( j=0; j<n_hits; j++ ) {
 			if ( (hits[j].h == h) && (hits[j].k == k)
@@ -530,12 +544,13 @@ int find_projected_peaks(struct image *image, UnitCell *cell)
 }
 
 
-int peak_sanity_check(struct image *image, UnitCell *cell)
+int peak_sanity_check(struct image *image, UnitCell *cell,
+                      int circular_domain, double domain_r)
 {
 	int i;
 	int n_sane = 0;
 
-	find_projected_peaks(image, cell);
+	find_projected_peaks(image, cell, circular_domain, domain_r);
 	if ( image->n_hits == 0 ) return 0;  /* Failed sanity check: no peaks */
 
 	for ( i=0; i<image->n_hits; i++ ) {
@@ -564,7 +579,7 @@ int peak_sanity_check(struct image *image, UnitCell *cell)
 
 void output_intensities(struct image *image, UnitCell *cell,
                         pthread_mutex_t *mutex, int polar, int sa,
-                        int use_closer)
+                        int use_closer, int circular_domain, double domain_r)
 {
 	int i;
 	int n_found;
@@ -577,7 +592,9 @@ void output_intensities(struct image *image, UnitCell *cell,
 	double csx, csy, csz;
 	double a, b, c, al, be, ga;
 
-	if ( image->n_hits == 0 ) find_projected_peaks(image, cell);
+	if ( image->n_hits == 0 ) {
+		find_projected_peaks(image, cell, circular_domain, domain_r);
+	}
 	if ( image->n_hits == 0 ) return;
 
 	/* Get exclusive access to the output stream if necessary */
-- 
cgit v1.2.3