aboutsummaryrefslogtreecommitdiff
path: root/src/peaks.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-06-07 17:48:47 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:27 +0100
commitb7050c2de9a6e990fac39554d65018b641397dee (patch)
treed94970d78b065a5ccc24afcdfff971ec1f3b89cc /src/peaks.c
parent6c06f74f4ab6d8089f55355abb4b137267f8ac6f (diff)
Get rid of integrate_pixels()
Diffstat (limited to 'src/peaks.c')
-rw-r--r--src/peaks.c165
1 files changed, 0 insertions, 165 deletions
diff --git a/src/peaks.c b/src/peaks.c
index e4a7cfb1..c542e44c 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -572,168 +572,3 @@ void integrate_reflections(struct image *image, int polar, int use_closer,
}
}
-
-
-RefList *integrate_pixels(struct image *image, int circular_domain,
- double domain_r, int do_polar)
-{
- int i;
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- int fs, ss;
- double aslen, bslen, cslen;
- double *intensities;
- double *xmom;
- double *ymom;
- ReflItemList *obs;
- RefList *reflections;
-
- obs = new_items();
- intensities = new_list_intensity();
- xmom = new_list_intensity();
- ymom = new_list_intensity();
- reflections = reflist_new();
-
- /* "Borrow" direction values to get reciprocal lengths */
- cell_get_reciprocal(image->indexed_cell, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz);
- aslen = modulus(ax, ay, az);
- bslen = modulus(bx, by, bz);
- cslen = modulus(cx, cy, cz);
-
- cell_get_cartesian(image->indexed_cell, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz);
- /* For each pixel */
- fesetround(1); /* Round towards nearest */
- 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 */
- signed int h, k, l;
- struct rvec q;
- double dist;
- struct panel *p;
- double twotheta;
-
- p = find_panel(image->det, fs, ss);
- if ( p == NULL ) continue;
- if ( p->no_index ) continue;
-
- 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;
- ld = q.u * cx + q.v * cy + q.w * cz;
-
- h = lrint(hd);
- k = lrint(kd);
- l = lrint(ld);
-
- dh = hd - h; dk = kd - k; dl = ld - l;
-
- if ( circular_domain ) {
-
- /* Circular integration domain */
- dist = sqrt(pow(dh*aslen, 2.0) + pow(dk*bslen, 2.0)
- + pow(dl*cslen, 2.0));
-
- } else {
-
- /* "Crystallographic" integration domain */
- dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0));
- }
-
- if ( dist < domain_r ) {
-
- double val;
- double pix_area, Lsq, proj_area, dsq, sa;
- double phi, pa, pb, pol;
- double xs, ys, rx, ry;
-
- /* Veto if we want to integrate a bad region */
- if ( image->flags != NULL ) {
- int flags;
- flags = image->flags[fs+image->width*ss];
- if ( !(flags & 0x01) ) 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 */
- proj_area = pix_area * cos(twotheta);
-
- /* Calculate distance from crystal to pixel */
- 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));
-
- /* Projected area of pixel / distance squared */
- sa = 1.0e7 * proj_area / (dsq + Lsq);
-
- /* Solid angle correction is needed in this case */
- val /= sa;
-
- if ( do_polar ) {
-
- 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;
-
- }
-
- /* Add value to sum */
- integrate_intensity(intensities, h, k, l, val);
-
- 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);
- }
-
- }
-
- }
- }
-
- for ( i=0; i<num_items(obs); i++ ) {
-
- struct refl_item *it;
- double intensity, xmomv, ymomv;
- double xp, yp;
- Reflection *refl;
-
- it = get_item(obs, i);
- intensity = lookup_intensity(intensities, it->h, it->k, it->l);
- xmomv = lookup_intensity(xmom, it->h, it->k, it->l);
- ymomv = lookup_intensity(ymom, it->h, it->k, it->l);
-
- xp = xmomv / (double)intensity;
- yp = ymomv / (double)intensity;
-
- refl = add_refl(reflections, it->h, it->k, it->l);
- set_int(refl, intensity);
- set_detector_pos(refl, 0.0, xp, yp);
-
- }
-
- free(xmom);
- free(ymom);
- free(intensities);
- delete_items(obs);
-
- return reflections;
-}