asdf.c: codestyle fixed

1 files changed, 469 insertions, 441 deletions

diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
index 0735a14f..6ad29645 100644
--- a/libcrystfel/src/asdf.c
+++ b/libcrystfel/src/asdf.c
@@ -26,24 +26,24 @@ struct asdf_private {
 
 /* Possible direct vector */
 struct tvector {
-    gsl_vector *t;
-    int n; // number of fitting reflections
-    int *fits; 
+	gsl_vector *t;
+	int n; // number of fitting reflections
+	int *fits; 
 };
 
 struct asdf_cell {
-    gsl_vector *axes[3];
-    gsl_vector *reciprocal[3];
-    
-    int n; // number of fitting reflections
-    double volume; 
-    
-    int N_refls; // total number of reflections
-    int *reflections; // reflections[i] = 1 if reflections fits
-    double **indices; // indices[i] = [h, k, l] for all reflections (not rounded)
-    
-    int acl; // minimum number of reflections fitting to one of the axes[]
-    int n_max; // maximum number of reflections fitting to some t-vector	
+	gsl_vector *axes[3];
+	gsl_vector *reciprocal[3];
+	
+	int n; // number of fitting reflections
+	double volume; 
+	
+	int N_refls; // total number of reflections
+	int *reflections; // reflections[i] = 1 if reflections fits
+	double **indices; // indices[i] = [h, k, l] for all reflections (not rounded)
+	
+	int acl; // minimum number of reflections fitting to one of the axes[]
+	int n_max; // maximum number of reflections fitting to some t-vector	
 };
 
 struct tvector tvector_new(int n) {
@@ -73,13 +73,13 @@ static int tvector_memcpy(struct tvector *dest, struct tvector *src, int n) {
 
 static int asdf_cell_free(struct asdf_cell c) {
 	int i;
-	for (i = 0; i < 3; i ++) {
+	for ( i = 0; i < 3; i++ ) {
 		gsl_vector_free(c.axes[i]);
 		gsl_vector_free(c.reciprocal[i]);
 	}
 	
 	free(c.reflections);
-	for (i = 0; i < c.N_refls; i ++) {
+	for ( i = 0; i < c.N_refls; i++ ) {
 		free(c.indices[i]);
 	}
 	free(c.indices);
@@ -92,7 +92,7 @@ static struct asdf_cell asdf_cell_new(int n) {
 	struct asdf_cell c;
 
 	int i;
-	for (i = 0; i < 3; i ++) {
+	for ( i = 0; i < 3; i++ ) {
 		c.axes[i] = gsl_vector_alloc(3);
 		c.reciprocal[i] = gsl_vector_alloc(3);
 	}
@@ -101,7 +101,7 @@ static struct asdf_cell asdf_cell_new(int n) {
 	c.reflections = malloc(sizeof(int) * n);
 	
 	c.indices = malloc(sizeof(int *) * n);
-	for (i = 0; i < n; i ++) {
+	for ( i = 0; i < n; i++ ) {
 		c.indices[i] = malloc(sizeof(int) * 3);
 	}
 	
@@ -115,7 +115,7 @@ static struct asdf_cell asdf_cell_new(int n) {
 
 static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) {
 	int i;
-	for (i = 0; i < 3; i ++) {
+	for ( i = 0; i < 3; i++ ) {
 		gsl_vector_memcpy(dest->axes[i], src->axes[i]);
 		gsl_vector_memcpy(dest->reciprocal[i], src->reciprocal[i]);
 	}
@@ -130,7 +130,7 @@ static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) {
 	
 	memcpy(dest->indices, src->indices, sizeof(int *) * n);
 	
-	for (i = 0; i < n; i ++) {
+	for (i  = 0; i < n; i++ ) {
 		memcpy(dest->indices[i], src->indices[i], sizeof(int) * 3);
 	}
 	
@@ -140,51 +140,52 @@ static int asdf_cell_memcpy(struct asdf_cell *dest, struct asdf_cell *src) {
 }
 
 /* result = vec1 cross vec2 */
-static int cross_product(gsl_vector *vec1, gsl_vector *vec2, gsl_vector **result) {
-    double c1[3], c2[3], p[3];
-    int i;
-    for (i = 0; i < 3; i += 1) {
-        c1[i] = gsl_vector_get(vec1, i);
-        c2[i] = gsl_vector_get(vec2, i);
-    }
-    
-    p[0] = c1[1] * c2[2] - c1[2] * c2[1];
-    p[1] = - c1[0] * c2[2] + c1[2] * c2[0];
-    p[2] = c1[0] * c2[1] - c1[1] * c2[0];
-    
-    for (i = 0; i < 3; i += 1) {
-        gsl_vector_set(*result, i, p[i]);
-    }
-    
-    return 1;
+static int cross_product(gsl_vector *vec1, gsl_vector *vec2, 
+                         gsl_vector **result) {
+	double c1[3], c2[3], p[3];
+	int i;
+	for ( i = 0; i < 3; i++ ) {
+		c1[i] = gsl_vector_get(vec1, i);
+		c2[i] = gsl_vector_get(vec2, i);
+	}
+
+	p[0] = c1[1] * c2[2] - c1[2] * c2[1];
+	p[1] = - c1[0] * c2[2] + c1[2] * c2[0];
+	p[2] = c1[0] * c2[1] - c1[1] * c2[0];
+
+	for ( i = 0; i < 3; i++ ) {
+		gsl_vector_set(*result, i, p[i]);
+	}
+
+	return 1;
 }
 
 /* Returns triple product of three gsl_vectors */ 
 static double calc_volume(gsl_vector *vec1, gsl_vector *vec2, gsl_vector *vec3) {
-    double volume;
-    gsl_vector *cross = gsl_vector_alloc(3);    
-    
-    cross_product(vec1, vec2, &cross);
-    gsl_blas_ddot(vec3, cross, &volume);
-    
-    gsl_vector_free(cross);
-    return volume;
+	double volume;
+	gsl_vector *cross = gsl_vector_alloc(3);    
+	
+	cross_product(vec1, vec2, &cross);
+	gsl_blas_ddot(vec3, cross, &volume);
+	
+	gsl_vector_free(cross);
+	return volume;
 }
 
 static int calc_reciprocal(gsl_vector **direct, gsl_vector **reciprocal) {
-    double volume;
-    
-    cross_product(direct[1], direct[2], &reciprocal[0]);
-    gsl_blas_ddot(direct[0], reciprocal[0], &volume);
-    gsl_vector_scale(reciprocal[0], 1/volume);
-    
-    cross_product(direct[2], direct[0], &reciprocal[1]);
-    gsl_vector_scale(reciprocal[1], 1/volume);
-    
-    cross_product(direct[0], direct[1], &reciprocal[2]);
-    gsl_vector_scale(reciprocal[2], 1/volume);
-    
-    return 1;
+	double volume;
+	
+	cross_product(direct[1], direct[2], &reciprocal[0]);
+	gsl_blas_ddot(direct[0], reciprocal[0], &volume);
+	gsl_vector_scale(reciprocal[0], 1/volume);
+	
+	cross_product(direct[2], direct[0], &reciprocal[1]);
+	gsl_vector_scale(reciprocal[1], 1/volume);
+	
+	cross_product(direct[0], direct[1], &reciprocal[2]);
+	gsl_vector_scale(reciprocal[2], 1/volume);
+	
+	return 1;
 }
 
 static int check_cell(struct asdf_private *dp, struct image *image,
@@ -230,21 +231,21 @@ static int check_cell(struct asdf_private *dp, struct image *image,
 
 static int compare_doubles (const void *a, const void *b)
 {
-  const double *da = (const double *) a;
-  const double *db = (const double *) b;
-
-  return (*da > *db) - (*da < *db);
+	const double *da = (const double *) a;
+	const double *db = (const double *) b;
+	
+	return (*da > *db) - (*da < *db);
 }
 
 /* Compares tvectors by length */
 static int compare_tvectors (const void *a, const void *b)
 {
-    struct tvector *ta = (struct tvector *) a;
-    struct tvector *tb = (struct tvector *) b;
+	struct tvector *ta = (struct tvector *) a;
+	struct tvector *tb = (struct tvector *) b;
 	
 	//~ if (ta->n == tb->n) {
-		return (gsl_blas_dnrm2(ta->t) > gsl_blas_dnrm2(tb->t)) - 
-			   (gsl_blas_dnrm2(ta->t) < gsl_blas_dnrm2(tb->t));
+	return (gsl_blas_dnrm2(ta->t) > gsl_blas_dnrm2(tb->t)) - 
+	       (gsl_blas_dnrm2(ta->t) < gsl_blas_dnrm2(tb->t));
 	//~ } 
 	//~ 
 	//~ return (ta->n > tb->n) - (ta->n < tb->n);
@@ -253,75 +254,75 @@ static int compare_tvectors (const void *a, const void *b)
 /* Calculates normal to a triplet c1, c2, c3. Returns 0 if reflections are on
  * the same line */
 static int calc_normal(gsl_vector *c1, gsl_vector *c2, gsl_vector *c3, 
-                gsl_vector *normal) {
+                       gsl_vector *normal) {
     
-    gsl_vector *c12 = gsl_vector_alloc(3);
-    gsl_vector *c23 = gsl_vector_alloc(3);
-    gsl_vector *c31 = gsl_vector_alloc(3);
-    gsl_vector *res = gsl_vector_alloc(3);
-    
-    cross_product(c1, c2, &c12);
-    cross_product(c2, c3, &c23);
-    cross_product(c3, c1, &c31);
-    
-    int i;
-    for (i = 0; i < 3; i += 1) {
-        gsl_vector_set(res, i, gsl_vector_get(c12, i) +
-                               gsl_vector_get(c23, i) +
-                               gsl_vector_get(c31, i));
-    }  
-    
-    gsl_vector_free(c12);
-    gsl_vector_free(c23);
-    gsl_vector_free(c31);
-    
-    double norm = gsl_blas_dnrm2(res);
-    if (norm < 0.0001) {
-        gsl_vector_free(res);
-        return 0;
-    } else {
-        gsl_vector_scale(res, 1/norm);
-        gsl_vector_memcpy(normal, res);
-        gsl_vector_free(res);
-    }
-    
-    return 1;
+	gsl_vector *c12 = gsl_vector_alloc(3);
+	gsl_vector *c23 = gsl_vector_alloc(3);
+	gsl_vector *c31 = gsl_vector_alloc(3);
+	gsl_vector *res = gsl_vector_alloc(3);
+	
+	cross_product(c1, c2, &c12);
+	cross_product(c2, c3, &c23);
+	cross_product(c3, c1, &c31);
+	
+	int i;
+	for ( i = 0; i < 3; i++ ) {
+		gsl_vector_set(res, i, gsl_vector_get(c12, i) +
+		                       gsl_vector_get(c23, i) +
+		                       gsl_vector_get(c31, i));
+	}  
+	
+	gsl_vector_free(c12);
+	gsl_vector_free(c23);
+	gsl_vector_free(c31);
+	
+	double norm = gsl_blas_dnrm2(res);
+	if ( norm < 0.0001 ) {
+		gsl_vector_free(res);
+		return 0;
+	} else {
+		gsl_vector_scale(res, 1/norm);
+		gsl_vector_memcpy(normal, res);
+		gsl_vector_free(res);
+	}
+	
+	return 1;
 }
 
 static float find_ds_fft(double *projections, int N_projections, double d_max) {    
  
-    int n = N_projections; 
-    double projections_sorted[n];
-    memcpy(projections_sorted, projections, sizeof(double) * n);    
-    qsort(projections_sorted, n, sizeof(double), compare_doubles);
-       
-    int i;
-    
-    int N = 1024; // number of points in fft calculation
-    
-    double *in = malloc(sizeof(double) * N);
-    fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N);
-    
-    for (i = 0; i < N; i++) {
-        in[i] = 0;
-    }
-    
-    for (i = 0; i < n; i++) {
-        in[(int)((projections_sorted[i] - projections_sorted[0]) / 
-                 (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++;
-    }
+	int n = N_projections; 
+	double projections_sorted[n];
+	memcpy(projections_sorted, projections, sizeof(double) * n);    
+	qsort(projections_sorted, n, sizeof(double), compare_doubles);
+	
+	int i;
+	
+	int N = 1024; // number of points in fft calculation
+	
+	double *in = malloc(sizeof(double) * N);
+	fftw_complex *out = fftw_malloc(sizeof (fftw_complex) * N);
+	
+	for ( i = 0; i < N; i++ ) {
+		in[i] = 0;
+	}
+	
+	for ( i = 0; i < n; i++ ) {
+		in[(int)((projections_sorted[i] - projections_sorted[0]) / 
+		   (projections_sorted[n - 1] - projections_sorted[0]) * N)] ++;
+	}
 		
-    fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);
-    fftw_execute(p);
-    fftw_destroy_plan(p);
-    
-    int i_max = (int)(d_max * (projections_sorted[n - 1] - 
-                               projections_sorted[0]));
-    
-    int d = 1;
-    double max = 0;
-    double a;
-    for (i = 1; i <= i_max; i ++) {
+	fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);
+	fftw_execute(p);
+	fftw_destroy_plan(p);
+	
+	int i_max = (int)(d_max * (projections_sorted[n - 1] - 
+	                           projections_sorted[0]));
+	
+	int d = 1;
+	double max = 0;
+	double a;
+	for ( i = 1; i <= i_max; i++ ) {
 		a = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]);
 		if (a > max) {
 			max = a;
@@ -343,17 +344,19 @@ static float find_ds_fft(double *projections, int N_projections, double d_max) {
 static int check_refl_fitting_ds(double *projections, int N_projections, 
                                  double ds, double LevelFit) {
 	
-	if (ds == 0) return 0;
+	if ( ds == 0 ) return 0;
 	
-    int i;
-    int n = 0;
-	for (i = 0; i < N_projections; i += 1) {
-		if (fabs(projections[i] - ds * round(projections[i]/ds)) < LevelFit) {
+	int i;
+	int n = 0;
+	for ( i = 0; i < N_projections; i++ ) {
+		if ( fabs(projections[i] - 
+		     ds * round(projections[i]/ds)) < LevelFit ) 
+		{
 			n += 1;
 		}
 	}
 	
-    return n; 
+	return n; 
 }
 
 /* Refines d*, writes 1 to fits[i] if the i-th projection fits d* */ 
@@ -370,11 +373,12 @@ static float refine_ds(double *projections, int N_projections, double ds,
 	
 	double c1, cov11, sumsq;
 	double ds_new = ds;
-	while (N_fits < N_new) {
+	while ( N_fits < N_new ) {
 		N_fits = 0;
-		for (i = 0; i < N_projections; i ++) {
-			if (fabs(projections[i] - ds_new * round(projections[i]/ds_new)) < 
-			                                                         LevelFit) {
+		for ( i = 0; i < N_projections; i++ ) {
+			if ( fabs(projections[i] - ds_new * 
+			     round(projections[i] / ds_new)) < LevelFit ) 
+			{
 				fit_refls[N_fits] = projections[i];
 				indices[N_fits] = round(projections[i]/ds_new);
 				N_fits ++;
@@ -385,9 +389,11 @@ static float refine_ds(double *projections, int N_projections, double ds,
 		}
 		
 					
-		gsl_fit_mul(indices, 1, fit_refls, 1, N_fits, &c1, &cov11, &sumsq);
-		N_new = check_refl_fitting_ds(projections, N_projections, c1, LevelFit);
-		if (N_new >= N_fits) {
+		gsl_fit_mul(indices, 1, fit_refls, 1, N_fits, &c1, &cov11, 
+		            &sumsq);
+		N_new = check_refl_fitting_ds(projections, N_projections, c1, 
+		                              LevelFit);
+		if ( N_new >= N_fits ) {
 			ds_new = c1;
 		}		
 	}
@@ -397,142 +403,146 @@ static float refine_ds(double *projections, int N_projections, double ds,
 
 static int check_refl_fitting_cell(struct asdf_cell *c, 
                                    gsl_vector **reflections, 
-                                   int N_reflections, double IndexFit) {    
-                                
-    double dist[3];
-    
-    calc_reciprocal(c->axes, c->reciprocal);
-    c->n = 0;
-    int i, j, k;
-    for(i = 0; i < N_reflections; i += 1) { 
-        
-        for (j = 0; j < 3; j += 1) dist[j] = 0;
-        
-        for (j = 0; j < 3; j += 1) {
-            gsl_blas_ddot(reflections[i], c->axes[j], &c->indices[i][j]);
-    
-            for (k = 0; k < 3; k += 1) {
-                dist[k] += gsl_vector_get(c->reciprocal[j], k) * 
-                           (c->indices[i][j] - round(c->indices[i][j]));
-            }
-        }
-        
-        /* A reflection fits if the distance (in reciprocal space) between the 
-         * observed and calculated reflection position is less than Indexfit */
-         
-        if (dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2] < 
-                                                         IndexFit * IndexFit ) {
-            c->reflections[i] = 1;
-            c->n += 1;
-        } else {
-            c->reflections[i] = 0;
-        }
-    }
-    
-    return 1;
+                                   int N_reflections, double IndexFit) {
+	      
+	double dist[3];
+	
+	calc_reciprocal(c->axes, c->reciprocal);
+	c->n = 0;
+	int i, j, k;
+	for( i = 0; i < N_reflections; i += 1 ) { 
+	
+		for ( j = 0; j < 3; j++ ) dist[j] = 0;
+	
+		for ( j = 0; j < 3; j++ ) {
+			gsl_blas_ddot(reflections[i], c->axes[j],
+			              &c->indices[i][j]);
+	
+			for ( k = 0; k < 3; k++ ) {
+				dist[k] += gsl_vector_get(c->reciprocal[j], k) * 
+				(c->indices[i][j] - round(c->indices[i][j]));
+			}
+		}
+		
+		/* A reflection fits if the distance (in reciprocal space) 
+		 * between the observed and calculated reflection position 
+		 * is less than Indexfit */
+		
+		if ( dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2] < 
+							 IndexFit * IndexFit ) {
+			c->reflections[i] = 1;
+			c->n++;
+		} else {
+			c->reflections[i] = 0;
+		}
+	}
+	
+	return 1;
 }
 
 static void print_asdf_cell(struct asdf_cell cc) {
-    double a, b, c, alpha, beta, gamma, ab, bc, ca;
-    
-    a = gsl_blas_dnrm2(cc.axes[0]);
-    b = gsl_blas_dnrm2(cc.axes[1]);
-    c = gsl_blas_dnrm2(cc.axes[2]);
-    
-    gsl_blas_ddot(cc.axes[0], cc.axes[1], &ab);
-    gsl_blas_ddot(cc.axes[1], cc.axes[2], &bc);
-    gsl_blas_ddot(cc.axes[0], cc.axes[2], &ca);
-    
-    alpha = acos(bc/b/c)/M_PI*180;
-    beta = acos(ca/a/c)/M_PI*180;
-    gamma = acos(ab/a/b)/M_PI*180;
-    
-    //~ int i, j;
-    //~ for (i = 0; i < 3; i ++) {
+	double a, b, c, alpha, beta, gamma, ab, bc, ca;
+	
+	a = gsl_blas_dnrm2(cc.axes[0]);
+	b = gsl_blas_dnrm2(cc.axes[1]);
+	c = gsl_blas_dnrm2(cc.axes[2]);
+	
+	gsl_blas_ddot(cc.axes[0], cc.axes[1], &ab);
+	gsl_blas_ddot(cc.axes[1], cc.axes[2], &bc);
+	gsl_blas_ddot(cc.axes[0], cc.axes[2], &ca);
+	
+	alpha = acos(bc/b/c)/M_PI*180;
+	beta = acos(ca/a/c)/M_PI*180;
+	gamma = acos(ab/a/b)/M_PI*180;
+	
+	//~ int i, j;
+	//~ for (i = 0; i < 3; i ++) {
 		//~ for (j = 0; j < 3; j ++) {
 			//~ printf("%f ", gsl_vector_get(cc.axes[i], j));
 		//~ }
 		//~ printf("\n");
 	//~ }
 	
-    printf("%.2f %.2f %.2f %.2f %.2f %.2f %.0f %d \n", a, b, c, 
-                                                       alpha, beta, gamma, 
-                                                       cc.volume, cc.n);
-    
+	printf("%.2f %.2f %.2f %.2f %.2f %.2f %.0f %d \n", a, b, c, 
+						           alpha, beta, gamma, 
+						           cc.volume, cc.n);
+
 }
 
+/* Returns 0 when refinement doesn't converge (i.e. all fitting reflections  
+ * lie in the same plane) */
 static int refine_asdf_cell(struct asdf_cell c, gsl_vector **reflections, 
                             int N_reflections, double IndexFit) {
     
-    gsl_matrix *X = gsl_matrix_alloc(c.n, 3);
-    
-    gsl_vector *r[] = {gsl_vector_alloc(c.n), 
-		               gsl_vector_alloc(c.n), 
-		               gsl_vector_alloc(c.n)};
-    
-    gsl_vector *res = gsl_vector_alloc(3);
-    gsl_matrix *cov = gsl_matrix_alloc (3, 3);
-    double chisq;
-    
-    int i, j;
-    int n = 0;
-    for (i = 0; i < N_reflections; i += 1) {
-        if (c.reflections[i] == 1) {
-            for (j = 0; j < 3; j += 1) {
-                gsl_matrix_set(X, n, j, round(c.indices[i][j]));
-                gsl_vector_set(r[j], n, gsl_vector_get(reflections[i], j));
-            }
-            n += 1;
-        }
-    }
-    
-    gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(c.n, 3);
-    
-    for (i = 0; i < 3; i += 1) {
-        gsl_multifit_linear (X, r[i], res, cov, &chisq, work);
-        
-        for (j = 0; j < 3; j += 1) {
-            gsl_vector_set(c.reciprocal[j], i, gsl_vector_get(res, j));
-        }
-    }
-
-    calc_reciprocal(c.reciprocal, c.axes);
+	gsl_matrix *X = gsl_matrix_alloc(c.n, 3);
+	
+	gsl_vector *r[] = {gsl_vector_alloc(c.n), 
+	gsl_vector_alloc(c.n), 
+	gsl_vector_alloc(c.n)};
+	
+	gsl_vector *res = gsl_vector_alloc(3);
+	gsl_matrix *cov = gsl_matrix_alloc (3, 3);
+	double chisq;
+	
+	int i, j;
+	int n = 0;
+	for ( i = 0; i < N_reflections; i++ ) if ( c.reflections[i] == 1 ) 
+	{
+		for ( j = 0; j < 3; j++ ) {
+			gsl_matrix_set(X, n, j, round(c.indices[i][j]));
+			gsl_vector_set(r[j], n, 
+			               gsl_vector_get(reflections[i], j));
+		}
+		n++;
+	}
+	
+	gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(c.n, 3);
+	
+	for ( i = 0; i < 3; i++ ) {
+		gsl_multifit_linear (X, r[i], res, cov, &chisq, work);
+		
+		for (j = 0; j < 3; j++ ) {
+			gsl_vector_set(c.reciprocal[j], i, 
+			               gsl_vector_get(res, j));
+		}
+	}
+	
+	calc_reciprocal(c.reciprocal, c.axes);
 	
 	double a[3];
-	for (i = 0; i < 3; i ++) {
+	for ( i = 0; i < 3; i++ ) {
 		a[i] = gsl_blas_dnrm2(c.axes[i]);
 	}
-	    
-    gsl_multifit_linear_free(work);
-    gsl_vector_free(res);
-    gsl_matrix_free(cov);
-    gsl_matrix_free(X);
-    for (i = 0; i < 3; i += 1) {
-        gsl_vector_free(r[i]);
-    }
-    
-    if (fabs(a[0]) > 10000 || fabs(a[1]) > 10000 || fabs(a[2]) > 10000 
-                                                 || isnan(a[0])) {
+	
+	gsl_multifit_linear_free(work);
+	gsl_vector_free(res);
+	gsl_matrix_free(cov);
+	gsl_matrix_free(X);
+	for ( i = 0; i < 3; i++ ) {
+		gsl_vector_free(r[i]);
+	}
+	
+	if ( fabs(a[0]) > 10000 || fabs(a[1]) > 10000 || 
+	     fabs(a[2]) > 10000 || isnan(a[0]) ) 
+	{
 		return 0;
 	}
 	
-    return 1;
+	return 1;
 }
 
 static int reduce_asdf_cell(struct asdf_cell *cl) {
 	double a, b, c, alpha, beta, gamma, ab, bc, ca, bb, cc;
    
-    gsl_vector *va = gsl_vector_alloc(3);
-    
-    gsl_vector *vb = gsl_vector_alloc(3);
-    
-    gsl_vector *vc = gsl_vector_alloc(3);
+	gsl_vector *va = gsl_vector_alloc(3);	
+	gsl_vector *vb = gsl_vector_alloc(3);	
+	gsl_vector *vc = gsl_vector_alloc(3);
 	
 	int changed = 1;
-	while (changed) {
+	while ( changed ) {
 		changed = 0;
     
-        gsl_vector_memcpy(va, cl->axes[0]);
+		gsl_vector_memcpy(va, cl->axes[0]);
 		gsl_vector_memcpy(vb, cl->axes[1]);
 		gsl_vector_memcpy(vc, cl->axes[2]);
 	
@@ -548,9 +558,9 @@ static int reduce_asdf_cell(struct asdf_cell *cl) {
 		beta = acos(ca/a/c)/M_PI*180;
 		gamma = acos(ab/a/b)/M_PI*180;
 	
-		if (changed == 0) {
+		if ( changed == 0 ) {
 
-			if (gamma < 90) {
+			if ( gamma < 90 ) {
 				gsl_vector_scale(vb, -1);
 				gamma = 180 - gamma;
 				alpha = 180 - alpha;
@@ -558,9 +568,9 @@ static int reduce_asdf_cell(struct asdf_cell *cl) {
 			
 			gsl_vector_add(vb, va);
 			bb = gsl_blas_dnrm2(vb);
-			if (bb < b) {
+			if ( bb < b ) {
 				b = bb;
-				if (a < b) {
+				if ( a < b ) {
 					gsl_vector_memcpy(cl->axes[1], vb);
 				} else {
 					gsl_vector_memcpy(cl->axes[1], va);
@@ -570,9 +580,9 @@ static int reduce_asdf_cell(struct asdf_cell *cl) {
 			} 
 		}
 		
-		if (changed == 0) {
+		if ( changed == 0 ) {
 			
-			if (beta < 90) {
+			if ( beta < 90 ) {
 				gsl_vector_scale(vc, -1);				
 				beta = 180 - beta;
 				alpha = 180 - alpha;
@@ -580,11 +590,11 @@ static int reduce_asdf_cell(struct asdf_cell *cl) {
 			
 			gsl_vector_add(vc, va);
 			cc = gsl_blas_dnrm2(vc);
-			if (cc < c) {
+			if ( cc < c ) {
 				c = cc;
-				if (b < c) {
+				if ( b < c ) {
 					gsl_vector_memcpy(cl->axes[2], vc);
-				} else if (a < c) {
+				} else if ( a < c ) {
 					gsl_vector_memcpy(cl->axes[1], vc);
 					gsl_vector_memcpy(cl->axes[2], vb);
 				} else {
@@ -596,8 +606,8 @@ static int reduce_asdf_cell(struct asdf_cell *cl) {
 			} 
 		}
 		
-		if (changed == 0) {
-			if (alpha < 90) {
+		if ( changed == 0 ) {
+			if ( alpha < 90 ) {
 				gsl_vector_scale(vc, -1);			
 				beta = 180 - beta;
 				alpha = 180 - alpha;
@@ -605,11 +615,11 @@ static int reduce_asdf_cell(struct asdf_cell *cl) {
 			
 			gsl_vector_add(vc, vb);
 			cc = gsl_blas_dnrm2(vc);
-			if (cc < c) {
+			if ( cc < c ) {
 				c = cc;
-				if (b < c) {
+				if ( b < c ) {
 					gsl_vector_memcpy(cl->axes[2], vc);
-				} else if (a < c) {
+				} else if ( a < c ) {
 					gsl_vector_memcpy(cl->axes[1], vc);
 					gsl_vector_memcpy(cl->axes[2], vb);
 				} else {
@@ -624,7 +634,7 @@ static int reduce_asdf_cell(struct asdf_cell *cl) {
 		
 	cross_product(cl->axes[0], cl->axes[1], &vc);
 	gsl_blas_ddot(vc, cl->axes[2], &cl->volume);
-	if (cl->volume < 0) {
+	if ( cl->volume < 0 ) {
 		gsl_vector_scale(cl->axes[2], -1);
 		cl->volume *= -1;
 	}
@@ -637,34 +647,36 @@ static int reduce_asdf_cell(struct asdf_cell *cl) {
 }
 
 static int check_cell_angles(gsl_vector *va, gsl_vector *vb, gsl_vector *vc, 
-                             double max_cos) {
+                             double max_cos) 
+{
 	double a, b, c, cosa, cosb, cosg, ab, bc, ca;
 	
-    a = gsl_blas_dnrm2(va);
-    b = gsl_blas_dnrm2(vb);
-    c = gsl_blas_dnrm2(vc);
-    
-    gsl_blas_ddot(va, vb, &ab);
-    gsl_blas_ddot(vb, vc, &bc);
-    gsl_blas_ddot(vc, va, &ca);
-    
-    cosa = bc/b/c;
-    cosb = ca/a/c;
-    cosg = ab/a/b;
-    
-    if (fabs(cosa) > max_cos || fabs(cosb) > max_cos || fabs(cosg) > max_cos) {
+	a = gsl_blas_dnrm2(va);
+	b = gsl_blas_dnrm2(vb);
+	c = gsl_blas_dnrm2(vc);
+	
+	gsl_blas_ddot(va, vb, &ab);
+	gsl_blas_ddot(vb, vc, &bc);
+	gsl_blas_ddot(vc, va, &ca);
+	
+	cosa = bc/b/c;
+	cosb = ca/a/c;
+	cosg = ab/a/b;
+	
+	if ( fabs(cosa) > max_cos || fabs(cosb) > max_cos || 
+	                            fabs(cosg) > max_cos ) {
 		return 0;
 	}
-
-    return 1;              
+	
+	return 1;              
 }
 
 /* Returns min(t1.n, t2.n, t3.n) */
 static int find_acl(struct tvector t1, struct tvector t2, struct tvector t3) {
 	int i = t1.n, j = t2.n, k = t3.n;
-	if (i <= j && i <= k) return i;
-	if (j <= i && j <= k) return j;
-	if (k <= i && k <= j) return k;
+	if ( i <= j && i <= k ) return i;
+	if ( j <= i && j <= k ) return j;
+	if ( k <= i && k <= j ) return k;
 }
 
 static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit, 
@@ -676,110 +688,118 @@ static int find_cell(struct tvector *tvectors, int N_tvectors, double IndexFit,
 	
 	double volume;
 	
-	// only tvectors with number of fitting reflections > acl are considered
+	/* Only tvectors with the number of fitting reflections > acl are 
+	 * considered */
 	int acl = N_reflections < 18 ? 6 : N_reflections/3;
 	
 	struct asdf_cell c = asdf_cell_new(N_reflections);
     
-    // traversing a 3d array in slices perpendicular to the main diagonal
-    int sl;
-    for (sl = 0; sl < 3 * N_tvectors - 1; sl++) {
-        
-        int i_min = sl < 2 * N_tvectors ? 0 : sl - 2 * N_tvectors;
-        int i_max = sl < N_tvectors ? sl : N_tvectors;
+	/* Traversing a 3d array in slices perpendicular to the main diagonal */
+	int sl;
+	for ( sl = 0; sl < 3 * N_tvectors - 1; sl++ ) {
         
-        for (i = i_min; i < i_max; i++) if (tvectors[i].n > acl) {
-            
-            int j_min = sl - N_tvectors - 2 * i - 1 < 0 ? i + 1 : sl - N_tvectors - i;
-            int j_max = sl - N_tvectors - i < 0 ? sl - i : N_tvectors;
-            
-
-            for (j = j_min; j < j_max; j++) if (tvectors[j].n > acl) {
-                
-                k = sl - i - j - 1;
-                
-                if (k > j && tvectors[k].n > acl &&
-                             check_cell_angles(tvectors[i].t, 
-                                               tvectors[j].t, 
-                                               tvectors[k].t, 0.99) ) {
+	int i_min = sl < 2 * N_tvectors ? 0 : sl - 2 * N_tvectors;
+	int i_max = sl < N_tvectors ? sl : N_tvectors;
+	
+	for ( i = i_min; i < i_max; i++) if (tvectors[i].n > acl ) {
+	
+	int j_min = sl - N_tvectors - 2 * i - 1 < 0 ? i + 1 : sl - N_tvectors - i;
+	int j_max = sl - N_tvectors - i < 0 ? sl - i : N_tvectors;
+		
+	for ( j = j_min; j < j_max; j++) if (tvectors[j].n > acl ) {
+	
+		k = sl - i - j - 1;
+		
+		if ( k > j && tvectors[k].n > acl &&
+		     check_cell_angles(tvectors[i].t, 
+				       tvectors[j].t, 
+				       tvectors[k].t, 0.99) ) 
+		{
+		
+			volume = calc_volume(tvectors[i].t, 
+			                     tvectors[j].t, 
+			                     tvectors[k].t);
+			
+			if ( fabs(volume) > volume_min && 
+			     fabs(volume) < volume_max )
+			{
+				
+				gsl_vector_memcpy(c.axes[0], tvectors[i].t);
+				gsl_vector_memcpy(c.axes[1], tvectors[j].t);
+				gsl_vector_memcpy(c.axes[2], tvectors[k].t);
+				
+				c.volume = volume;
+				check_refl_fitting_cell(&c, reflections, 
+				                        N_reflections, 
+				                        IndexFit);
 				
-					volume = calc_volume(tvectors[i].t, 
-					                     tvectors[j].t, 
-					                     tvectors[k].t);
+				if ( c.n < 6 ) break;
+				
+				reduce_asdf_cell(&c);
+				
+				/* if one of the cell angles > 135 or < 45 
+				 * do not continue */
+				if ( !check_cell_angles(c.axes[0], c.axes[1], 
+				     c.axes[2], 0.71) ) break;
+				
+				/* index reflections with new cell axes */
+				check_refl_fitting_cell(&c, reflections, 
+				                        N_reflections, 
+				                        IndexFit);
+				
+				acl = find_acl(tvectors[i], 
+				               tvectors[j], 
+				               tvectors[k]);
+				
+				c.acl = acl;
+				c.n_max = n_max;
+				
+				/* refine cell until the number of fitting 
+				 * reflections stops increasing */
+				n = 0;
+				cell_correct = 1;
+				while ( c.n - n && cell_correct ) {
+					n = c.n;
+					cell_correct = refine_asdf_cell(c, 
+					                          reflections, 
+					                          N_reflections, 
+					                          IndexFit);
 					
-					if (fabs(volume) > volume_min && fabs(volume) < volume_max){
-						
-						gsl_vector_memcpy(c.axes[0], tvectors[i].t);
-						gsl_vector_memcpy(c.axes[1], tvectors[j].t);
-						gsl_vector_memcpy(c.axes[2], tvectors[k].t);
-						
-						c.volume = volume;
-						check_refl_fitting_cell(&c, reflections, N_reflections, 
-											    IndexFit);
-						
-						if (c.n >= 6) {
-							reduce_asdf_cell(&c);
-                            
-                            // if one of the cell angles > 135 or < 45 do not continue
-                            if (!check_cell_angles(c.axes[0], c.axes[1], 
-                                                   c.axes[2], 0.71)) break;
-							
-							// index reflections with new cell axes
-							check_refl_fitting_cell(&c, reflections, 
-							                        N_reflections, IndexFit);
-							
-							acl = find_acl(tvectors[i], 
-						                   tvectors[j], 
-						                   tvectors[k]);
-								
-							c.acl = acl;
-                            c.n_max = n_max;
-                            
-                            // refine cell until number of fitting reflections 
-                            // stops increasing
-							n = 0;
-							cell_correct = 1;
-							while (c.n - n && cell_correct) {
-								n = c.n;
-								cell_correct = refine_asdf_cell(c, reflections, 
-								                                N_reflections, 
-								                                IndexFit);
-								                                
-								check_refl_fitting_cell(&c, reflections, 
-								                       N_reflections, IndexFit);
-							}
-							
-							if (cell_correct) {
-       							reduce_asdf_cell(&c);
-       							if (result->n < c.n) 
-       							   asdf_cell_memcpy(result, &c);      
-       							
-								acl++;                             
-                                
-                                if (acl > n_max) break;
-								if (tvectors[j].n <= acl || 
-								    tvectors[i].n <= acl) break;
-							}
-						}
+					check_refl_fitting_cell(&c, reflections, 
+					                        N_reflections, 
+					                        IndexFit);
+				}
+				
+				if ( cell_correct ) {
+					reduce_asdf_cell(&c);
+					if ( result->n < c.n ) {
+						asdf_cell_memcpy(result, &c);      
 					}
+					acl++;                             
+				
+					if (acl > n_max) break;
+					if (tvectors[j].n <= acl || 
+					    tvectors[i].n <= acl) break;
 				}
 			}
-            if (acl > n_max) break;
-			if (tvectors[i].n <= acl) break;
 		}
-        if (acl > n_max) break;
+	}
+	if ( acl > n_max ) break;
+	if ( tvectors[i].n <= acl ) break;
+	}
+	if ( acl > n_max ) break;
 	}
 	
 	asdf_cell_free(c);
 	
-	if (result->n) return 1;
+	if ( result->n ) return 1;
 	return 0;
 }
 
 static void shuffle_triplets(int **triplets, int n) {
 	int i, j;
 	int t[3]; 
-	for (i = 0; i < n - 1; i++) {
+	for ( i = 0; i < n - 1; i++ ) {
 		j = i + rand() / (RAND_MAX / (n - i) + 1);
 		memcpy(t, triplets[j], 3 * sizeof(int));
 		memcpy(triplets[j], triplets[i], 3 * sizeof(int));
@@ -793,9 +813,10 @@ static double angle_between_gsl(gsl_vector *a, gsl_vector *b) {
 	return acos(ab/gsl_blas_dnrm2(a)/gsl_blas_dnrm2(b)) * 180 / M_PI;
 }
 
-static int index_refls(gsl_vector **reflections, int N_reflections, double d_max, 
-                       double volume_min, double volume_max, double LevelFit, 
-                       double IndexFit, int i_max, struct asdf_cell *c) {
+static int index_refls(gsl_vector **reflections, int N_reflections, 
+                       double d_max, double volume_min, double volume_max,
+                       double LevelFit, double IndexFit, int i_max, 
+                       struct asdf_cell *c) {
 
 	int i, j, k, l, n;
 	
@@ -805,9 +826,9 @@ static int index_refls(gsl_vector **reflections, int N_reflections, double d_max
 
 	int **triplets = malloc(N_triplets * sizeof(int *));
 	l = 0;
-	for (i = 0; i < N_reflections; i++) {
-		for (j = i + 1; j < N_reflections; j++) {
-			for (k = j + 1; k < N_reflections; k++) {
+	for ( i = 0; i < N_reflections; i++ ) {
+		for ( j = i + 1; j < N_reflections; j++ ) {
+			for ( k = j + 1; k < N_reflections; k++ ) {
 				triplets[l] = malloc(3 * sizeof(int));
 				
 				triplets[l][0] = i;
@@ -819,7 +840,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections, double d_max
 	}
 
 	/* Triplets are processed in a random sequence if N_triplets > 10000 */
-	if (N_reflections > 40) shuffle_triplets(triplets, N_triplets);
+	if ( N_reflections > 40 ) shuffle_triplets(triplets, N_triplets);
 	
 	gsl_vector *normal = gsl_vector_alloc(3);
 		
@@ -828,41 +849,47 @@ static int index_refls(gsl_vector **reflections, int N_reflections, double d_max
 	
 	int *fits = malloc(N_reflections * sizeof(int));
 	
-	if (i_max > N_triplets) i_max = N_triplets;
+	if ( i_max > N_triplets ) i_max = N_triplets;
    
-    struct tvector *tvectors = malloc(i_max * sizeof(struct tvector));
+	struct tvector *tvectors = malloc(i_max * sizeof(struct tvector));
 	int N_tvectors = 0;
 	
 	int n_max = 0; // maximum number of reflections fitting one of tvectors 
 
-	for (i = 0; i < i_max; i++) {
-		if (calc_normal(reflections[triplets[i][0]],
-		                reflections[triplets[i][1]],
-		                reflections[triplets[i][2]],
-		                normal)) {
+	for ( i = 0; i < i_max; i++ ) {
+		if ( calc_normal(reflections[triplets[i][0]],
+		                 reflections[triplets[i][1]],
+		                 reflections[triplets[i][2]],
+		                 normal) ) 
+		{
 			
-			// calculate projections of reflections to normal
-			for (k = 0; k < N_reflections; k++) {
-				gsl_blas_ddot(normal, reflections[k], &projections[k]);
+			/* Calculate projections of reflections to normal */
+			for ( k = 0; k < N_reflections; k++ ) {
+				gsl_blas_ddot(normal, reflections[k], 
+				              &projections[k]);
 			}
 			
-			// find ds - period in 1d lattice of projections
+			/* Find ds - period in 1d lattice of projections */
 			ds = find_ds_fft(projections, N_reflections, d_max);
 			
-			// refine ds, write 1 to fits[i] if reflections[i] fits ds
-			ds = refine_ds(projections, N_reflections, ds, LevelFit, fits);
+			/* Refine ds, write 1 to fits[i] if reflections[i] 
+			 * fits ds */
+			ds = refine_ds(projections, N_reflections, ds, LevelFit, 
+			               fits);
 			
-			// n - number of reflections fitting ds
-			n = check_refl_fitting_ds(projections, N_reflections, ds, LevelFit);
+			/* n - number of reflections fitting ds */
+			n = check_refl_fitting_ds(projections, N_reflections, 
+			                          ds, LevelFit);
 			
-			// normal/ds - possible direct vector
+			/* normal/ds - possible direct vector */
 			gsl_vector_scale(normal, 1/ds);
 			
-			if (n > N_reflections/3 && n > 6) {
+			if ( n > N_reflections/3 && n > 6 ) {
                 
 				tvectors[N_tvectors] = tvector_new(N_reflections);
 				
-				gsl_vector_memcpy(tvectors[N_tvectors].t, normal);
+				gsl_vector_memcpy(tvectors[N_tvectors].t, 
+				                  normal);
 				memcpy(tvectors[N_tvectors].fits, fits, 
 				       N_reflections * sizeof(int));
 				
@@ -874,37 +901,39 @@ static int index_refls(gsl_vector **reflections, int N_reflections, double d_max
 			}
 		}
 		
-		if ((i != 0 && i % 10000 == 0) || i == i_max - 1) {
+		if ( (i != 0 && i % 10000 == 0) || i == i_max - 1 ) {
 			
-			// sort tvectors by length
+			/* Sort tvectors by length */
 			qsort(tvectors, N_tvectors, sizeof(struct tvector), 
 			      compare_tvectors);
             
-            /* three shortest independent tvectors with t.n > acl determine
-             * the final cell. acl is selected for the solution with the
-             * maximum number of fitting reflections */
-             
-			find_cell(tvectors, N_tvectors, IndexFit, volume_min, volume_max, 
-			          n_max, reflections, N_reflections, c);
+			/* Three shortest independent tvectors with t.n > acl 
+			 * determine the final cell. acl is selected for the 
+			 * solution with the maximum number of fitting 
+			 * reflections */
+		             
+			find_cell(tvectors, N_tvectors, IndexFit, volume_min, 
+			          volume_max, n_max, reflections, 
+			          N_reflections, c);
                         
-			if (c->n > 4 * n_max / 5) {
+			if ( c->n > 4 * n_max / 5 ) {
 				break;
 			}
 		}	
 	}	
 	free(fits);
 	
-	for (i = 0; i < N_tvectors; i++) {
+	for ( i = 0; i < N_tvectors; i++ ) {
 		tvector_free(tvectors[i]);
 	}
 	free(tvectors);
 	
-	for (i = 0; i < N_triplets; i++) {
+	for ( i = 0; i < N_triplets; i++ ) {
 		free(triplets[i]);
 	}
 	free(triplets);
 
-	if (c->n) return 1;
+	if ( c->n ) return 1;
 	
 	return 0;
 }
@@ -936,30 +965,30 @@ double cell_get_volume(UnitCell *cell)
 }
 
 int run_asdf(struct image *image, IndexingPrivate *ipriv) {
-    int i;
+	int i;
 	
-    double LevelFit = 1./1000;
-    double IndexFit = 1./500;
-    double d_max = 220.; // thrice the maximum expected axis length
-    double volume_min = 100.;
-    double volume_max = 1000000.;
-    
-    int i_max = 10000; // maximum number of triplets 
-    
-    struct asdf_private *dp = (struct asdf_private *)ipriv;
-    
-    if ( dp->indm & INDEXING_CHECK_CELL_AXES ) {
+	double LevelFit = 1./1000;
+	double IndexFit = 1./500;
+	double d_max = 220.; // thrice the maximum expected axis length
+	double volume_min = 100.;
+	double volume_max = 1000000.;
+	
+	int i_max = 10000; // maximum number of triplets 
+	
+	struct asdf_private *dp = (struct asdf_private *)ipriv;
+	
+	if ( dp->indm & INDEXING_CHECK_CELL_AXES ) {
 		double volume = cell_get_volume(dp->template);
 		volume_min = volume * 0.95;
 		volume_max = volume * 1.05;
 	}
 	
-    int N_reflections = image_feature_count(image->features);
-    gsl_vector *reflections[N_reflections];
-
-    for (i = 0; i < N_reflections; i++) {
+	int N_reflections = image_feature_count(image->features);
+	gsl_vector *reflections[N_reflections];
+	
+	for ( i = 0; i < N_reflections; i++ ) {
 		struct imagefeature *f;
-
+	
 		f = image_get_feature(image->features, i);
 		if ( f == NULL ) continue;
 		
@@ -969,38 +998,37 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv) {
 		gsl_vector_set(reflections[i], 2, f->rz/1e10);	
 	}
 	
-    struct asdf_cell c = asdf_cell_new(N_reflections);
-    
-    if (N_reflections == 0) return 0;
-    
-    i = index_refls(reflections, N_reflections, d_max, volume_min, volume_max, 
-                    LevelFit, IndexFit, i_max, &c);
-                    //~ 
-    
-    for (i = 0; i < N_reflections; i ++) {
+	struct asdf_cell c = asdf_cell_new(N_reflections);
+	
+	if ( N_reflections == 0 ) return 0;
+	
+	i = index_refls(reflections, N_reflections, d_max, volume_min, volume_max, 
+		        LevelFit, IndexFit, i_max, &c);
+	
+	for ( i = 0; i < N_reflections; i++ ) {
 		gsl_vector_free(reflections[i]);
 	}
-    
-    if (i) {
+	
+	if ( i ) {
 		UnitCell *uc;
 		uc = cell_new();
 		 
 		cell_set_cartesian(uc, gsl_vector_get(c.axes[0], 0) * 1e-10, 
-							   gsl_vector_get(c.axes[0], 1) * 1e-10,
-							   gsl_vector_get(c.axes[0], 2) * 1e-10,
-							   gsl_vector_get(c.axes[1], 0) * 1e-10,
-							   gsl_vector_get(c.axes[1], 1) * 1e-10,
-							   gsl_vector_get(c.axes[1], 2) * 1e-10,
-							   gsl_vector_get(c.axes[2], 0) * 1e-10,
-							   gsl_vector_get(c.axes[2], 1) * 1e-10,
-							   gsl_vector_get(c.axes[2], 2) * 1e-10);
-
+				       gsl_vector_get(c.axes[0], 1) * 1e-10,
+				       gsl_vector_get(c.axes[0], 2) * 1e-10,
+				       gsl_vector_get(c.axes[1], 0) * 1e-10,
+				       gsl_vector_get(c.axes[1], 1) * 1e-10,
+				       gsl_vector_get(c.axes[1], 2) * 1e-10,
+				       gsl_vector_get(c.axes[2], 0) * 1e-10,
+				       gsl_vector_get(c.axes[2], 1) * 1e-10,
+				       gsl_vector_get(c.axes[2], 2) * 1e-10);
+	
 		if ( check_cell(dp, image, uc) ) {
 			cell_free(uc);
 			return 1;
 		}
 		
-        cell_free(uc);
+	cell_free(uc);
 	}
 	
 	return 0;