aboutsummaryrefslogtreecommitdiff
path: root/src/hrs-scaling.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-01-24 18:05:37 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:11 +0100
commite007d9e7cd0f943ffe17c2b069a9e0efceed0d2c (patch)
tree5d2f82ebaf6ba8e045bf33217b93f7892cbfb4df /src/hrs-scaling.c
parent384bf68d6946a62ebc33f632b924869fcb4e2e43 (diff)
Revert "Remove row-crossing stuff since it should not be needed"
This reverts commit a043b3cb7cd57da3f5f7feb9d470d1e771f5b664. Conflicts: src/hrs-scaling.c
Diffstat (limited to 'src/hrs-scaling.c')
-rw-r--r--src/hrs-scaling.c32
1 files changed, 24 insertions, 8 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 86a5ffab..eee459a2 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -137,19 +137,23 @@ static double iterate_scale(struct image *images, int n,
gsl_vector *shifts;
int l;
double max_shift;
+ int crossed = 0;
int n_ref;
- M = gsl_matrix_calloc(n, n);
- v = gsl_vector_calloc(n);
+ M = gsl_matrix_calloc(n-1, n-1);
+ v = gsl_vector_calloc(n-1);
n_ref = num_items(obs);
for ( l=0; l<n; l++ ) { /* "Equation number": one equation per frame */
int m; /* Frame (scale factor) number */
int h;
+ int lcomp;
double vc_tot = 0.0;
struct image *imagel = &images[l];
+ if ( l == crossed ) continue;
+
/* Determine the "solution" vector component */
for ( h=0; h<n_ref; h++ ) {
@@ -172,15 +176,19 @@ static double iterate_scale(struct image *images, int n,
}
- gsl_vector_set(v, l, vc_tot);
+ lcomp = l;
+ if ( l > crossed ) lcomp--;
+ gsl_vector_set(v, lcomp, vc_tot);
/* Now fill in the matrix components */
for ( m=0; m<n; m++ ) {
double mc_tot = 0.0;
+ int mcomp;
struct image *imagem = &images[m];
if ( m > l ) continue; /* Matrix is symmetric */
+ if ( m == crossed ) continue;
for ( h=0; h<n_ref; h++ ) {
@@ -211,14 +219,17 @@ static double iterate_scale(struct image *images, int n,
}
- gsl_matrix_set(M, l, m, mc_tot);
- gsl_matrix_set(M, m, l, mc_tot);
+ mcomp = m;
+ if ( m > crossed ) mcomp--;
+ gsl_matrix_set(M, lcomp, mcomp, mc_tot);
+ gsl_matrix_set(M, mcomp, lcomp, mc_tot);
}
}
- show_matrix_eqn(M, v, n);
+ show_matrix_eqn(M, v, n-1);
+
gsl_eigen_symmv_workspace *work;
gsl_vector *e_val;
@@ -235,14 +246,19 @@ static double iterate_scale(struct image *images, int n,
show_eigen(e_vec, e_val, n);
#if 0 /* HRS method */
- shifts = gsl_vector_alloc(n);
+ shifts = gsl_vector_alloc(n-1);
+
gsl_linalg_HH_solve(M, v, shifts);
max_shift = 0.0;
for ( l=0; l<n-1; l++ ) {
double shift = gsl_vector_get(shifts, l);
+ int limg;
+
+ limg = l;
+ if ( limg >= crossed ) limg++;
- images[l].osf += shift;
+ images[limg].osf += shift;
if ( fabs(shift) > fabs(max_shift) ) {
max_shift = fabs(shift);