diff options
| author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2008-08-12 16:30:06 +0000 |
|---|---|---|
| committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2008-08-12 16:30:06 +0000 |
| commit | 4ea397e11edfb5ef711e52a0d4088943df052d7b (patch) | |
| tree | 3d0e44b6c798a156dd85728d9c583d30ed191688 /src/reproject.c | |
| parent | d963362e6cb208ca730cf45b2b72ab3e86489907 (diff) | |
Improve numerical stability of reprojection
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@285 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src/reproject.c')
| -rw-r--r-- | src/reproject.c | 11 |
1 files changed, 7 insertions, 4 deletions
diff --git a/src/reproject.c b/src/reproject.c index e1ce13e..d071f9c 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -85,7 +85,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * double xl, yl, zl; double a, b, c; - double A1, A2, s1, s2, s; + double A1, A2, s1, s2, s, temp; /* Get the coordinates of the reciprocal lattice point */ xl = reflection->x; @@ -95,9 +95,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * /* Next, solve the relrod equation to calculate the excitation error */ a = 1.0; b = -2.0*(xl*nx + yl*ny + zl*nz); - c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda); - A1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a); - A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a); + c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda); /* FIXME: Don't think this is stable */ + printf("%e %e\n", xl*xl + yl*yl + zl*zl, 1.0/(image->lambda*image->lambda)); + /* Solve the quadratic equation */ + temp = -0.5 * (b + sign(b)*sqrt(b*b - 4*a*c)); + A1 = temp / a; + A2 = c / temp; s1 = 1.0/image->lambda - A1; s2 = 1.0/image->lambda - A2; if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2; |