diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/refine.c | 87 |
1 files changed, 57 insertions, 30 deletions
diff --git a/src/refine.c b/src/refine.c index 45aec8e..7871dc1 100644 --- a/src/refine.c +++ b/src/refine.c @@ -43,11 +43,38 @@ void refine_do_sequence(ControlContext *ctx) { } +static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) { + + double fom = 0.0; + int f; + + for ( f=0; f<nf; f++ ) { + + double xdf, ydf, zdf; + + xdf = d[f].h*s[i].dax + d[f].k*s[i].dbx + d[f].l*s[i].dcx; + ydf = d[f].h*s[i].day + d[f].k*s[i].dby + d[f].l*s[i].dcy; + zdf = d[f].h*s[i].daz + d[f].k*s[i].dbz + d[f].l*s[i].dcz; + xdf -= d[f].dx; + ydf -= d[f].dy; + zdf -= d[f].dz; + + fom += sqrt(xdf*xdf + ydf*ydf + zdf*zdf); + + } + + return fom; + +} + +static void refine_simplex_reflect(SimplexVertex *s, int v_worst) { +} + void refine_do_cell(ControlContext *ctx) { SimplexVertex s[10]; Deviation *d; - const double delta = 0.1e9; + double delta; int i, nf, f, v_worst; double fom_worst; @@ -61,20 +88,6 @@ void refine_do_cell(ControlContext *ctx) { return; } - /* Initialise the simplex */ - s[0].dax = 0.0; s[0].dbx = 0.0; s[0].dcx = 0.0; - s[0].day = 0.0; s[0].dby = 0.0; s[0].dcy = 0.0; - s[0].daz = 0.0; s[0].dbz = 0.0; s[0].dcz = 0.0; - memcpy(&s[1], &s[0], sizeof(SimplexVertex)); s[1].dax = delta; - memcpy(&s[2], &s[0], sizeof(SimplexVertex)); s[2].dbx = delta; - memcpy(&s[3], &s[0], sizeof(SimplexVertex)); s[3].dcx = delta; - memcpy(&s[4], &s[0], sizeof(SimplexVertex)); s[4].day = delta; - memcpy(&s[5], &s[0], sizeof(SimplexVertex)); s[5].dby = delta; - memcpy(&s[6], &s[0], sizeof(SimplexVertex)); s[6].dcy = delta; - memcpy(&s[7], &s[0], sizeof(SimplexVertex)); s[7].daz = delta; - memcpy(&s[8], &s[0], sizeof(SimplexVertex)); s[8].dbz = delta; - memcpy(&s[9], &s[0], sizeof(SimplexVertex)); s[9].dcz = delta; - /* Create the table of indicies and deviations */ nf = 0; for ( i=0; i<ctx->images->n_images; i++ ) { @@ -100,7 +113,7 @@ void refine_do_cell(ControlContext *ctx) { image = &ctx->images->images[i]; - for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) { + for ( j=0; j<ctx->images->images[i].rflist->n_features; j++ ) { ImageFeature *rf; double dix, diy, dx, dy; @@ -117,7 +130,7 @@ void refine_do_cell(ControlContext *ctx) { /* Determine the difference vector */ dix = rf->partner->x - rf->x; diy = rf->partner->y - rf->y; - printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, d[f].h, d[f].k, d[f].l, dix, diy); + //printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, d[f].h, d[f].k, d[f].l, dix, diy); old_x = rf->partner->x; old_y = rf->partner->y; @@ -127,7 +140,7 @@ void refine_do_cell(ControlContext *ctx) { mapping_rotate(dx, dy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); rf->partner->x = old_x; rf->partner->y = old_y; - printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9); + //printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9); d[f].dx = dlx; d[f].dy = dly; @@ -140,27 +153,41 @@ void refine_do_cell(ControlContext *ctx) { } assert( f == nf ); + /* Initialise the simplex */ + delta = 0.01e9; + s[0].dax = 0.0; s[0].dbx = 0.0; s[0].dcx = 0.0; + s[0].day = 0.0; s[0].dby = 0.0; s[0].dcy = 0.0; + s[0].daz = 0.0; s[0].dbz = 0.0; s[0].dcz = 0.0; + memcpy(&s[1], &s[0], sizeof(SimplexVertex)); s[1].dax = delta; + memcpy(&s[2], &s[0], sizeof(SimplexVertex)); s[2].dbx = delta; + memcpy(&s[3], &s[0], sizeof(SimplexVertex)); s[3].dcx = delta; + memcpy(&s[4], &s[0], sizeof(SimplexVertex)); s[4].day = delta; + memcpy(&s[5], &s[0], sizeof(SimplexVertex)); s[5].dby = delta; + memcpy(&s[6], &s[0], sizeof(SimplexVertex)); s[6].dcy = delta; + memcpy(&s[7], &s[0], sizeof(SimplexVertex)); s[7].daz = delta; + memcpy(&s[8], &s[0], sizeof(SimplexVertex)); s[8].dbz = delta; + memcpy(&s[9], &s[0], sizeof(SimplexVertex)); s[9].dcz = delta; + /* Find the least favourable vertex of the simplex */ v_worst = 0; - fom_worst = 0; + fom_worst = 0.0; for ( i=0; i<10; i++ ) { - double fom = 0; - int j; + double fom; - for ( j=0; j<nf; j++ ) { + fom = refine_mean_dev(d, nf, s, i); - double xdev, ydev, zdev; - - xdev = 0; - ydev = 0; - zdev = 0; - - fom += sqrt(xdev*xdev + ydev*ydev + zdev*zdev); - + printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, (fom/nf)/1e9); + if ( fom > fom_worst ) { + fom_worst = fom; + v_worst = i; } } + printf("The worst vertex is number %i\n", v_worst); + + /* Reflect this vertex across the opposite face */ + refine_simplex_reflect(s, v_worst); ctx->images->images[ctx->dw->cur_image].rflist = NULL; reproject_lattice_changed(ctx); |