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 | |
| 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
| -rw-r--r-- | src/reproject.c | 11 | ||||
| -rw-r--r-- | src/utils.c | 5 | ||||
| -rw-r--r-- | src/utils.h | 1 |
3 files changed, 13 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; diff --git a/src/utils.c b/src/utils.c index 3fe2e23..30a6e34 100644 --- a/src/utils.c +++ b/src/utils.c @@ -112,4 +112,9 @@ void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v, const char *pr } +int sign(double a) { + if ( a < 0 ) return -1; + if ( a > 0 ) return +1; + return 0; +} diff --git a/src/utils.h b/src/utils.h index 5736739..5c2e586 100644 --- a/src/utils.h +++ b/src/utils.h @@ -30,6 +30,7 @@ extern double distance3d(double x1, double y1, double z1, double x2, double y2, extern size_t skipspace(const char *s); extern void chomp(char *s); extern void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v, const char *prefix); +extern int sign(double a); #define rad2deg(a) ((a)*180/M_PI) #define deg2rad(a) ((a)*M_PI/180) |