diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/refine.c | 59 | ||||
| -rw-r--r-- | src/reflections.c | 2 |
2 files changed, 31 insertions, 30 deletions
diff --git a/src/refine.c b/src/refine.c index 90a02b1..4e5a3ab 100644 --- a/src/refine.c +++ b/src/refine.c @@ -95,9 +95,6 @@ static void refine_fix_unconstrained(gsl_matrix *M) { ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { ImageFeatureList *flist; - gsl_matrix *M; - gsl_vector *q; - gsl_vector *p; int j, ns; Axis axis; double dax = 0.0, dbx = 0.0, dcx = 0.0; @@ -105,19 +102,24 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { double dlx = 0.0, dly = 0.0, dlz = 0.0; flist = image->features; - - M = gsl_matrix_alloc(3, 3); - p = gsl_vector_alloc(3); - + for ( axis = AXIS_X; axis <= AXIS_Y; axis++ ) { gsl_permutation *perm; int s; double det; + gsl_matrix *M; + gsl_vector *q; + gsl_vector *p; + gsl_vector *scratch; + gsl_matrix *V; + gsl_vector *S; if ( axis == AXIS_X ) printf("RF: ------------------------------------------------------------------ Refining x coordinates -----\n"); if ( axis == AXIS_Y ) printf("RF: ------------------------------------------------------------------ Refining y coordinates -----\n"); + M = gsl_matrix_alloc(3, 3); + p = gsl_vector_alloc(3); gsl_matrix_set_zero(M); gsl_vector_set_zero(p); @@ -144,15 +146,6 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { diy = rf->partner->y - rf->y; printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, h, k, l, dix, diy); - //if ( (h<0) || (k<0) || (l<0) ) { - // h = -h; - // k = -k; - // l = -l; - // dix = -dix; - // diy = -diy; - // printf("\n Fudge: %3i %3i %3i dev = %+9.5f %+9.5f px ", h, k, l, dix, diy); - //} - old_x = rf->partner->x; old_y = rf->partner->y; rf->partner->x = dix + rf->partner->parent->x_centre; @@ -197,17 +190,30 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { return NULL; } - /* Do the fitting */ + /* Prepare for solving */ refine_fix_unconstrained(M); matrix_vector_show(M, p, "RF: "); + + /* Calculate and display the determinant */ perm = gsl_permutation_alloc(M->size1); gsl_linalg_LU_decomp(M, perm, &s); det = gsl_linalg_LU_det(M, s); printf("RF: Determinant of M = %f\n", det); + gsl_permutation_free(perm); + /* Solve the equations */ + V = gsl_matrix_alloc(M->size2, M->size2); + S = gsl_vector_alloc(M->size2); + scratch = gsl_vector_alloc(M->size2); + gsl_linalg_SV_decomp(M, V, S, scratch); q = gsl_vector_alloc(3); /* This is the "answer" */ gsl_vector_set_zero(q); - gsl_linalg_LU_solve(M, perm, p, q); + gsl_linalg_SV_solve(M, V, S, p, q); + gsl_vector_free(scratch); + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_matrix_free(M); + gsl_vector_free(p); switch ( axis ) { case AXIS_X : { @@ -224,14 +230,10 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { } } - gsl_permutation_free(perm); - + gsl_vector_free(q); + } - gsl_matrix_free(M); - gsl_vector_free(p); - gsl_vector_free(q); - 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); @@ -239,18 +241,15 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { /* Update the cell */ mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - 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); + printf("RF: a changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9); 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("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); + printf("RF: b changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9); 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("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); + printf("RF: c changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9); cell->c.x += dlx; cell->c.y += dly; cell->c.z += dlz; return NULL; diff --git a/src/reflections.c b/src/reflections.c index 8bfa7ff..58a7992 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -235,6 +235,8 @@ ReflectionList *reflection_list_from_cell(Basis *basis) { double x, y, z; + if ( h != 0 ) continue; + x = h*basis->a.x + k*basis->b.x + l*basis->c.x; y = h*basis->a.y + k*basis->b.y + l*basis->c.y; z = h*basis->a.z + k*basis->b.z + l*basis->c.z; |