From e2fe7137e7994b11bfce1f3a74af8f5ec9025b48 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Tue, 11 Jul 2023 11:01:23 +0200
Subject: Redefine ray vector without trigonometry

This makes the gradients very simple!
---
 libcrystfel/src/geometry.c       | 12 +++------
 libcrystfel/src/predict-refine.c | 57 +++++++++-------------------------------
 2 files changed, 15 insertions(+), 54 deletions(-)

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 9dc471ee..674fe4d0 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -53,17 +53,11 @@ static int locate_peak_on_panel(double x, double y, double z, double k,
                                 double det_shift_x, double det_shift_y,
                                 double *pfs, double *pss)
 {
-	double ctt, tta, phi;
 	gsl_vector *v;
 	gsl_vector *t;
 	gsl_matrix *M;
 	double fs, ss, one_over_mu;
 
-	/* Calculate 2theta (scattering angle) and azimuth (phi) */
-	tta = atan2(sqrt(x*x+y*y), k+z);
-	ctt = cos(tta);
-	phi = atan2(y, x);
-
 	/* Set up matrix equation */
 	M = gsl_matrix_alloc(3, 3);
 	v = gsl_vector_alloc(3);
@@ -73,9 +67,9 @@ static int locate_peak_on_panel(double x, double y, double z, double k,
 		return 0;
 	}
 
-	gsl_vector_set(t, 0, sin(tta)*cos(phi));
-	gsl_vector_set(t, 1, sin(tta)*sin(phi));
-	gsl_vector_set(t, 2, ctt);
+	gsl_vector_set(t, 0, x);
+	gsl_vector_set(t, 1, y);
+	gsl_vector_set(t, 2, k+z);
 
 	gsl_matrix_set(M, 0, 0, p->cnx+(det_shift_x/p->pixel_pitch));
 	gsl_matrix_set(M, 0, 1, p->fsx);
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index d8dbbca0..b87173b7 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -148,76 +148,49 @@ int fs_ss_gradient_physics(int param, Reflection *refl, UnitCell *cell,
                          float *fsg, float *ssg)
 {
 	signed int h, k, l;
-	double xl, yl, zl, kpred;
-	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
 	gsl_vector *dRdp;
 	gsl_vector *v;
 
 	get_indices(refl, &h, &k, &l);
-	kpred = get_kpred(refl);
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-	xl = h*asx + k*bsx + l*csx;
-	yl = h*asy + k*bsy + l*csy;
-	zl = h*asz + k*bsz + l*csz;
 
 	dRdp = gsl_vector_calloc(3);
 
 	switch ( param ) {
 
 		case GPARAM_ASX :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 0, h);
 		break;
 
 		case GPARAM_BSX :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 0, k);
 		break;
 
 		case GPARAM_CSX :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 0, l);
 		break;
 
 		case GPARAM_ASY :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 1, h);
 		break;
 
 		case GPARAM_BSY :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 1, k);
 		break;
 
 		case GPARAM_CSY :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 1, l);
 		break;
 
 		case GPARAM_ASZ :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 2, h);
 		break;
 
 		case GPARAM_BSZ :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 2, k);
 		break;
 
 		case GPARAM_CSZ :
-		gsl_vector_set(dRdp, 0, 0.0);
-		gsl_vector_set(dRdp, 1, 0.0);
-		gsl_vector_set(dRdp, 2, 0.0);
+		gsl_vector_set(dRdp, 2, l);
 		break;
 
 		default :
@@ -321,7 +294,6 @@ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell,
 	signed int h, k, l;
 	double xl, yl, zl, kpred;
 	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
-	double tta, ctt, phi;
 	gsl_vector *t;
 	gsl_vector *v;
 	gsl_matrix *M;
@@ -337,11 +309,6 @@ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell,
 	yl = h*asy + k*bsy + l*csy;
 	zl = h*asz + k*bsz + l*csz;
 
-	/* Calculate 2theta (scattering angle) and azimuth (phi) */
-	tta = atan2(sqrt(xl*xl+yl*yl), kpred+zl);
-	ctt = cos(tta);
-	phi = atan2(yl, xl);
-
 	/* Set up matrix equation */
 	M = gsl_matrix_alloc(3, 3);
 	v = gsl_vector_alloc(3);
@@ -351,9 +318,9 @@ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell,
 		return 1;
 	}
 
-	gsl_vector_set(t, 0, sin(tta)*cos(phi));
-	gsl_vector_set(t, 1, sin(tta)*sin(phi));
-	gsl_vector_set(t, 2, ctt);
+	gsl_vector_set(t, 0, xl);
+	gsl_vector_set(t, 1, yl);
+	gsl_vector_set(t, 2, kpred+zl);
 
 	gsl_matrix_set(M, 0, 0, p->cnx);
 	gsl_matrix_set(M, 0, 1, p->fsx);
-- 
cgit v1.2.3