diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/refine.c | 8 | ||||
| -rw-r--r-- | src/refinetest2d.c | 28 | ||||
| -rw-r--r-- | src/reproject.c | 10 |
3 files changed, 22 insertions, 24 deletions
diff --git a/src/refine.c b/src/refine.c index 987fba1..5de8df3 100644 --- a/src/refine.c +++ b/src/refine.c @@ -33,10 +33,10 @@ #define DISPFACTOR 1.0e9 /* Number of parameters */ -#define NUM_PARAMS 4 +#define NUM_PARAMS 9 /* Refine debug */ -#define REFINE_DEBUG 1 +#define REFINE_DEBUG 0 /* A simplex is an array of ten of these */ typedef struct { @@ -408,8 +408,6 @@ double refine_do_cell(ControlContext *ctx) { d[f].dy = dly; d[f].dz = dlz; - printf("Feature %3i: %8f %8f %8f\n", f, dlx/DISPFACTOR, dly/DISPFACTOR, dlz/DISPFACTOR); - f++; } @@ -434,7 +432,7 @@ double refine_do_cell(ControlContext *ctx) { memcpy(&s[9], &s[0], sizeof(SimplexVertex)); s[9].dcz = delta; /* Iterate */ - maxiter = 5000; + maxiter = 500; for ( it=0; it<maxiter; it++ ) { double conv; diff --git a/src/refinetest2d.c b/src/refinetest2d.c index 03b65d6..cdc0c8d 100644 --- a/src/refinetest2d.c +++ b/src/refinetest2d.c @@ -30,17 +30,17 @@ static int check_cell(Basis *cell, Basis *cell_real) { int fail; -// printf(" Calculated Actual\n"); -// printf("---------------------------\n"); -// printf("ax %+8f %+8f\n", cell->a.x/1e9, cell_real->a.x/1e9); -// printf("ay %+8f %+8f\n", cell->a.y/1e9, cell_real->a.y/1e9); -// printf("az %+8f %+8f\n", cell->a.z/1e9, cell_real->a.z/1e9); -// printf("bx %+8f %+8f\n", cell->b.x/1e9, cell_real->b.x/1e9); -// printf("by %+8f %+8f\n", cell->b.y/1e9, cell_real->b.y/1e9); -// printf("bz %+8f %+8f\n", cell->b.z/1e9, cell_real->b.z/1e9); -// printf("cx %+8f %+8f\n", cell->c.x/1e9, cell_real->c.x/1e9); -// printf("cy %+8f %+8f\n", cell->c.y/1e9, cell_real->c.y/1e9); -// printf("cz %+8f %+8f\n", cell->c.z/1e9, cell_real->c.z/1e9); + printf(" Refinement Actual\n"); + printf("---------------------------\n"); + printf("ax %+8f %+8f\n", cell->a.x/1e9, cell_real->a.x/1e9); + printf("ay %+8f %+8f\n", cell->a.y/1e9, cell_real->a.y/1e9); + printf("az %+8f %+8f\n", cell->a.z/1e9, cell_real->a.z/1e9); + printf("bx %+8f %+8f\n", cell->b.x/1e9, cell_real->b.x/1e9); + printf("by %+8f %+8f\n", cell->b.y/1e9, cell_real->b.y/1e9); + printf("bz %+8f %+8f\n", cell->b.z/1e9, cell_real->b.z/1e9); + printf("cx %+8f %+8f\n", cell->c.x/1e9, cell_real->c.x/1e9); + printf("cy %+8f %+8f\n", cell->c.y/1e9, cell_real->c.y/1e9); + printf("cz %+8f %+8f\n", cell->c.z/1e9, cell_real->c.z/1e9); fail = 0; if ( fabs(cell->a.x - cell_real->a.x) > 0.01e9 ) { @@ -107,7 +107,7 @@ int main(int argc, char *argv[]) { ctx->x_centre = 256; ctx->y_centre = 256; ctx->pixel_size = 5e7; - image_add(ctx->images, NULL, 512, 512, deg2rad(90.0), ctx); + image_add(ctx->images, NULL, 512, 512, deg2rad(0.0), ctx); ctx->images->images[0].features = image_feature_list_new(); /* Fudge to avoid horrifying pointer-related death */ @@ -122,7 +122,7 @@ int main(int argc, char *argv[]) { /* The "real" reflections */ reflections_real = reflection_list_from_cell(cell_real); ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real); - printf("RT: %i test features generated.\n", ctx->images->images[0].features->n_features); +// printf("RT: %i test features generated.\n", ctx->images->images[0].features->n_features); /* The "model" cell to be refined */ ctx->cell = malloc(sizeof(Basis)); @@ -145,6 +145,8 @@ int main(int argc, char *argv[]) { if ( fail ) return 1; + printf("\n2D refinement test OK.\n"); + return 0; } diff --git a/src/reproject.c b/src/reproject.c index 439c2f7..e1ce13e 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, temp; + double A1, A2, s1, s2, s; /* Get the coordinates of the reciprocal lattice point */ xl = reflection->x; @@ -95,11 +95,9 @@ 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); /* FIXME: Don't think this is stable */ - /* Solve the quadratic equation */ - temp = -0.5 * (b + sign(b)*sqrt(b*b - 4*a*c)); - A1 = temp / a; - A2 = c / temp; + 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); s1 = 1.0/image->lambda - A1; s2 = 1.0/image->lambda - A2; if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2; |