diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2008-01-16 17:50:50 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2008-01-16 17:50:50 +0000 |
commit | 762a08d3ddba51df373f443e9499c5b850d9f50c (patch) | |
tree | 8a858369f7ebba18cad1f3edce49c03f4d334804 /src | |
parent | 1d42a788bcf53de056cd050559ec338ef5840525 (diff) |
More refinement work...
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@254 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src')
-rw-r--r-- | src/iprtest.c | 4 | ||||
-rw-r--r-- | src/refine.c | 35 | ||||
-rw-r--r-- | src/utils.c | 6 | ||||
-rw-r--r-- | src/utils.h | 2 |
4 files changed, 23 insertions, 24 deletions
diff --git a/src/iprtest.c b/src/iprtest.c index cb2fb10..7f9541d 100644 --- a/src/iprtest.c +++ b/src/iprtest.c @@ -117,8 +117,8 @@ int main(int argc, char *argv[]) { ctx->image = malloc(sizeof(ImageRecord)); ctx->image->image = NULL; - ctx->image->tilt = 20.0; - ctx->image->omega = 0.0; + ctx->image->tilt = 10.0; + ctx->image->omega = 90.0; ctx->image->slop = 0.0; ctx->image->fmode = FORMULATION_PIXELSIZE; ctx->image->pixel_size = 0.6e8; diff --git a/src/refine.c b/src/refine.c index 14681c9..90a02b1 100644 --- a/src/refine.c +++ b/src/refine.c @@ -66,7 +66,7 @@ static void refine_fix_unconstrained(gsl_matrix *M) { int zero2 = 1; - if ( row_altered ) break; + if ( row_altered ) break; /* This row is now fixed, so move on to the next */ for ( i=0; i<M->size1; i++ ) { if ( gsl_matrix_get(M, i, col) != 0.0 ) { @@ -79,7 +79,6 @@ static void refine_fix_unconstrained(gsl_matrix *M) { gsl_matrix_set(M, row, col, 1.0); modified = 1; row_altered = 1; - continue; } } @@ -116,6 +115,9 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { int s; double det; + if ( axis == AXIS_X ) printf("RF: ------------------------------------------------------------------ Refining x coordinates -----\n"); + if ( axis == AXIS_Y ) printf("RF: ------------------------------------------------------------------ Refining y coordinates -----\n"); + gsl_matrix_set_zero(M); gsl_vector_set_zero(p); @@ -197,11 +199,11 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { /* Do the fitting */ refine_fix_unconstrained(M); - matrix_vector_show(M, p); + matrix_vector_show(M, p, "RF: "); perm = gsl_permutation_alloc(M->size1); gsl_linalg_LU_decomp(M, perm, &s); det = gsl_linalg_LU_det(M, s); - printf("Determinant of M = %f\n", det); + printf("RF: Determinant of M = %f\n", det); q = gsl_vector_alloc(3); /* This is the "answer" */ gsl_vector_set_zero(q); @@ -230,31 +232,26 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { gsl_vector_free(p); gsl_vector_free(q); - printf("a should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dax/1e9, day/1e9); - printf("b should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dbx/1e9, dby/1e9); - printf("c should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dcx/1e9, dcy/1e9); + printf("RF: ------------------------------------------------------------------ Refinement results ---------\n"); + printf("RF: a should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dax/1e9, day/1e9); + printf("RF: b should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dbx/1e9, dby/1e9); + printf("RF: c should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dcx/1e9, dcy/1e9); /* Update the cell */ mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("a changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, + printf("RF: a changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, 100*dlx/cell->a.x, 100*dly/cell->a.y, 100*dlz/cell->a.z); - if ( dlx/cell->a.x < 0.1 ) cell->a.x += dlx; - if ( dly/cell->a.y < 0.1 ) cell->a.y += dly; - if ( dlz/cell->a.z < 0.1 ) cell->a.z += dlz; + cell->a.x += dlx; cell->a.y += dly; cell->a.z += dlz; mapping_rotate(dbx, dby, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("b changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, + printf("RF: b changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, 100*dlx/cell->b.x, 100*dly/cell->b.y, 100*dlz/cell->b.z); - if ( dlx/cell->b.x < 0.1 ) cell->b.x += dlx; - if ( dly/cell->b.y < 0.1 ) cell->b.y += dly; - if ( dlz/cell->b.z < 0.1 ) cell->b.z += dlz; + cell->b.x += dlx; cell->b.y += dly; cell->b.z += dlz; mapping_rotate(dcx, dcy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("c changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, + printf("RF: c changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, 100*dlx/cell->c.x, 100*dly/cell->c.y, 100*dlz/cell->c.z); - if ( dlx/cell->c.x < 0.1 ) cell->c.x += dlx; - if ( dly/cell->c.y < 0.1 ) cell->c.y += dly; - if ( dlz/cell->c.z < 0.1 ) cell->c.z += dlz; + cell->c.x += dlx; cell->c.y += dly; cell->c.z += dlz; return NULL; diff --git a/src/utils.c b/src/utils.c index 8fc3cf6..1e82abb 100644 --- a/src/utils.c +++ b/src/utils.c @@ -96,12 +96,14 @@ void chomp(char *s) { } -void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v) { +void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v, const char *prefix) { int i, j; + if ( prefix == NULL ) prefix = ""; + for ( i=0; i<m->size1; i++ ) { - printf("[ "); + printf("%s[ ", prefix); for ( j=0; j<m->size2; j++ ) { printf("%12.8f ", gsl_matrix_get(m, i, j)); } diff --git a/src/utils.h b/src/utils.h index 47e283a..804f1af 100644 --- a/src/utils.h +++ b/src/utils.h @@ -30,7 +30,7 @@ extern double lambda(double voltage); extern double distance3d(double x1, double y1, double z1, double x2, double y2, double z2); 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); +extern void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v, const char *prefix); #define rad2deg(a) ((a)*180/M_PI) #define deg2rad(a) ((a)*M_PI/180) |