pattern_sim: Implement phased gradients

13 files changed, 162 insertions, 34 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index e7752078..f58c5046 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -100,13 +100,13 @@ int main(int argc, char *argv[])
 
 	cell = load_cell_from_pdb("molecule.pdb");
 	c1 = new_list_count();
-	ref1 = read_reflections(afile, c1);
+	ref1 = read_reflections(afile, c1, NULL);
 	if ( ref1 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", afile);
 		return 1;
 	}
 	c2 = new_list_count();
-	ref2 = read_reflections(bfile, c2);
+	ref2 = read_reflections(bfile, c2, NULL);
 	if ( ref2 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", bfile);
 		return 1;
@@ -146,7 +146,7 @@ int main(int argc, char *argv[])
 	STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R*100.0, scale);
 
 	if ( outfile != NULL ) {
-		write_reflections(outfile, NULL, out, 1, cell, 1);
+		write_reflections(outfile, NULL, out, NULL, 1, cell, 1);
 	}
 
 	return 0;
diff --git a/src/diffraction.c b/src/diffraction.c
index 9bd7c4e5..4adf5932 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -148,9 +148,107 @@ static double interpolate_intensity(const double *ref,
 }
 
 
+static double complex interpolate_phased_linear(const double *ref,
+                                                const unsigned int *counts,
+                                                const double *phases,
+                                                float hd,
+                                                signed int k, signed int l)
+{
+	signed int h;
+	double val1, val2;
+	float f;
+	unsigned int c1, c2;
+	double ph1, ph2;
+	double re1, re2, im1, im2;
+	double re, im;
+
+	h = (signed int)hd;
+	if ( hd < 0.0 ) h -= 1;
+	f = hd - (float)h;
+	assert(f >= 0.0);
+
+	val1 = lookup_intensity(ref, h, k, l);
+	val2 = lookup_intensity(ref, h+1, k, l);
+	c1 = lookup_count(counts, h, k, l);
+	c2 = lookup_count(counts, h+1, k, l);
+	ph1 = lookup_phase(phases, h, k, l);
+	ph2 = lookup_phase(phases, h+1, k, l);
+
+	if ( c1 == 0 ) {
+		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+		      h, k, l);
+		return 1.0e20;
+	}
+
+	if ( c2 == 0 ) {
+		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+		      h+1, k, l);
+		return 1.0e20;
+	}
+
+	val1 = val1 / (double)c1;
+	val2 = val2 / (double)c2;
+
+	/* Calculate real and imaginary parts */
+	re1 = val1 * cos(ph1);
+	im1 = val1 * sin(ph1);
+	re2 = val2 * cos(ph2);
+	im2 = val2 * sin(ph2);
+
+	re = (1.0-f)*re1 + f*re2;
+	im = (1.0-f)*im1 + f*im2;
+
+	return re + im*I;
+}
+
+
+static double complex interpolate_phased_bilinear(const double *ref,
+                                                  const unsigned int *counts,
+                                                  const double *phases,
+                                                  float hd, float kd,
+                                                  signed int l)
+{
+	signed int k;
+	double complex val1, val2;
+	float f;
+
+	k = (signed int)kd;
+	if ( kd < 0.0 ) k -= 1;
+	f = kd - (float)k;
+	assert(f >= 0.0);
+
+	val1 = interpolate_phased_linear(ref, counts, phases, hd, k, l);
+	val2 = interpolate_phased_linear(ref, counts, phases, hd, k+1, l);
+
+	return (1.0-f)*val1 + f*val2;
+}
+
+
+static double interpolate_phased_intensity(const double *ref,
+                                           const unsigned int *counts,
+                                           const double *phases,
+                                           float hd, float kd, float ld)
+{
+	signed int l;
+	double complex val1, val2;
+	float f;
+
+	l = (signed int)ld;
+	if ( ld < 0.0 ) l -= 1;
+	f = ld - (float)l;
+	assert(f >= 0.0);
+
+	val1 = interpolate_phased_bilinear(ref, counts, phases, hd, kd, l);
+	val2 = interpolate_phased_bilinear(ref, counts, phases, hd, kd, l+1);
+
+	return cabs((1.0-f)*val1 + f*val2);
+}
+
+
 /* Look up the structure factor for the nearest Bragg condition */
 static double molecule_factor(const double *intensities,
-                              const unsigned int *counts, struct rvec q,
+                              const unsigned int *counts, const double *phases,
+                              struct rvec q,
                               double ax, double ay, double az,
                               double bx, double by, double bz,
                               double cx, double cy, double cz,
@@ -180,6 +278,9 @@ static double molecule_factor(const double *intensities,
 		r = interpolate_intensity(intensities, counts, hd, kd, ld);
 		break;
 	case GRADIENT_PHASED :
+		r = interpolate_phased_intensity(intensities, counts, phases,
+		                                 hd, kd, ld);
+		break;
 	default:
 		ERROR("This gradient method not implemented yet.\n");
 		exit(1);
@@ -290,7 +391,8 @@ double get_tt(struct image *image, unsigned int xs, unsigned int ys)
 
 void get_diffraction(struct image *image, int na, int nb, int nc,
                      const double *intensities, const unsigned int *counts,
-                     UnitCell *cell, int do_water, GradientMethod m)
+                     const double *phases, UnitCell *cell, int do_water,
+                     GradientMethod m)
 {
 	unsigned int xs, ys;
 	double ax, ay, az;
@@ -345,7 +447,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 				I_molecule = 1.0e10;
 			} else {
 				I_molecule = molecule_factor(intensities,
-				                             counts, q,
+				                             counts, phases, q,
 			                                     ax,ay,az,
 			                                     bx,by,bz,cx,cy,cz,
 				                             m);
diff --git a/src/diffraction.h b/src/diffraction.h
index ed583feb..c7afaa3e 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -27,8 +27,8 @@ typedef enum {
 
 extern void get_diffraction(struct image *image, int na, int nb, int nc,
                             const double *intensities,
-                            const unsigned int *counts, UnitCell *cell,
-                            int do_water, GradientMethod m);
+                            const unsigned int *counts, const double *phases,
+                            UnitCell *cell, int do_water, GradientMethod m);
 extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
                          unsigned int sampling, float *ttp, float k);
 
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 735ae52a..e6fe4cee 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -106,6 +106,7 @@ int main(int argc, char *argv[])
 {
 	int c;
 	double *ideal_ref;
+	double *phases;
 	struct molecule *mol;
 	char *template = NULL;
 	int config_noisify = 0;
@@ -177,10 +178,12 @@ int main(int argc, char *argv[])
 
 	mol = load_molecule(filename);
 	cts = new_list_count();
+	phases = new_list_intensity(); /* "intensity" type used for phases */
 	if ( input == NULL ) {
-		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9), cts);
+		ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
+		                            cts, phases);
 	} else {
-		ideal_ref = read_reflections(input, cts);
+		ideal_ref = read_reflections(input, cts, phases);
 		free(input);
 	}
 
@@ -243,7 +246,8 @@ int main(int argc, char *argv[])
 
 	}
 
-	write_reflections(output, counts, ideal_ref, config_za, mol->cell, 1);
+	write_reflections(output, counts, ideal_ref, phases,
+	                  config_za, mol->cell, 1);
 
 	return 0;
 }
diff --git a/src/indexamajig.c b/src/indexamajig.c
index c365048a..e9f936b1 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -216,7 +216,8 @@ static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
 		get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
 	} else {
 		get_diffraction(simage, 24, 24, 40,
-		                intensities, counts, cell, 0, GRADIENT_MOSAIC);
+		                intensities, counts, NULL, cell, 0,
+		                GRADIENT_MOSAIC);
 	}
 	record_image(simage, 0);
 
@@ -508,7 +509,7 @@ int main(int argc, char *argv[])
 
 	if ( intfile != NULL ) {
 		counts = new_list_count();
-		intensities = read_reflections(intfile, counts);
+		intensities = read_reflections(intfile, counts, NULL);
 	} else {
 		intensities = NULL;
 		counts = NULL;
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 8c7a7e4b..c89da879 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -167,6 +167,7 @@ int main(int argc, char *argv[])
 	double *powder;
 	char *intfile = NULL;
 	double *intensities;
+	double *phases;
 	int config_simdetails = 0;
 	int config_nearbragg = 0;
 	int config_randomquat = 0;
@@ -286,9 +287,11 @@ int main(int argc, char *argv[])
 		STATUS("I'll simulate a flat intensity distribution.\n");
 		intensities = NULL;
 		counts = NULL;
+		phases = NULL;
 	} else {
 		counts = new_list_count();
-		intensities = read_reflections(intfile, counts);
+		phases = new_list_phase();
+		intensities = read_reflections(intfile, counts, phases);
 		free(intfile);
 	}
 
@@ -379,7 +382,7 @@ int main(int argc, char *argv[])
 			get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
 		} else {
 			get_diffraction(&image, na, nb, nc, intensities, counts,
-			                cell, !config_nowater, grad);
+			                phases, cell, !config_nowater, grad);
 		}
 		if ( image.data == NULL ) {
 			ERROR("Diffraction calculation failed.\n");
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 7c42cd3e..306ee49f 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -176,7 +176,7 @@ static void process_reflections(double *ref, unsigned int *counts,
 	if ( do_zoneaxis ) {
 		char name[64];
 		snprintf(name, 63, "ZA-%u.dat", n_patterns);
-		write_reflections(name, counts, ref, 1, cell, 1);
+		write_reflections(name, counts, ref, NULL, 1, cell, 1);
 	}
 
 	fh = fopen("results/convergence.dat", "a");
@@ -323,7 +323,7 @@ int main(int argc, char *argv[])
 	if ( intfile != NULL ) {
 		truecounts = new_list_count();
 		STATUS("Comparing against '%s'\n", intfile);
-		trueref = read_reflections(intfile, truecounts);
+		trueref = read_reflections(intfile, truecounts, NULL);
 		free(intfile);
 	} else {
 		trueref = NULL;
@@ -457,13 +457,14 @@ int main(int argc, char *argv[])
 	}
 
 	if ( output != NULL ) {
-		write_reflections(output, model_counts, model, 0, cell, 1);
+		write_reflections(output, model_counts, model, NULL,
+		                  0, cell, 1);
 	}
 
 	if ( config_zoneaxis ) {
 		char name[64];
 		snprintf(name, 63, "ZA-%u.dat", n_patterns);
-		write_reflections(name, model_counts, model, 1, cell, 10);
+		write_reflections(name, model_counts, model, NULL, 1, cell, 10);
 	}
 
 	STATUS("There were %u patterns.\n", n_patterns);
diff --git a/src/reflections.c b/src/reflections.c
index c730977e..4b48c38f 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -22,8 +22,8 @@
 
 
 void write_reflections(const char *filename, unsigned int *counts,
-                       double *ref, int zone_axis, UnitCell *cell,
-                       unsigned int min_counts)
+                       double *ref, double *phases, int zone_axis,
+                       UnitCell *cell, unsigned int min_counts)
 {
 	FILE *fh;
 	signed int h, k, l;
@@ -50,7 +50,7 @@ void write_reflections(const char *filename, unsigned int *counts,
 		fprintf(fh, "angle %5.3f deg\n", rad2deg(alpha));
 		fprintf(fh, "scale 10\n");
 	} else {
-		fprintf(fh, " h   k   l    I    sigma(I)   1/d / nm^-1\n");
+		fprintf(fh, " h   k   l    I    phase(I)    sigma(I)   1/d / nm^-1\n");
 	}
 
 	for ( h=-INDMAX; h<INDMAX; h++ ) {
@@ -59,6 +59,7 @@ void write_reflections(const char *filename, unsigned int *counts,
 
 		int N;
 		double intensity, s;
+		char ph[32];
 
 		if ( counts ) {
 			N = lookup_count(counts, h, k, l);
@@ -68,7 +69,14 @@ void write_reflections(const char *filename, unsigned int *counts,
 		}
 		if ( zone_axis && (l != 0) ) continue;
 
-		intensity = lookup_intensity(ref, h, k, l) / N;
+		intensity = lookup_phase(ref, h, k, l) / N;
+		if ( phases != NULL ) {
+			double p;
+			p = lookup_intensity(phases, h, k, l);
+			snprintf(ph, 31, "%f", p);
+		} else {
+			strncpy(ph, "-", 31);
+		}
 
 		if ( cell != NULL ) {
 			s = 2.0*resolution(cell, h, k, l);
@@ -77,7 +85,7 @@ void write_reflections(const char *filename, unsigned int *counts,
 		}
 
 		/* h, k, l, I, sigma(I), s */
-		fprintf(fh, "%3i %3i %3i %f %f %f\n", h, k, l, intensity,
+		fprintf(fh, "%3i %3i %3i %f %s %f %f\n", h, k, l, intensity, ph,
 		                                      0.0, s/1.0e9);
 
 	}
@@ -87,7 +95,8 @@ void write_reflections(const char *filename, unsigned int *counts,
 }
 
 
-double *read_reflections(const char *filename, unsigned int *counts)
+double *read_reflections(const char *filename, unsigned int *counts,
+                         double *phases)
 {
 	double *ref;
 	FILE *fh;
@@ -105,14 +114,15 @@ double *read_reflections(const char *filename, unsigned int *counts)
 
 		char line[1024];
 		signed int h, k, l;
-		float intensity;
+		float intensity, ph;
 		int r;
 
 		rval = fgets(line, 1023, fh);
-		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
-		if ( r != 4 ) continue;
+		r = sscanf(line, "%i %i %i %f %f", &h, &k, &l, &intensity, &ph);
+		if ( r != 5 ) continue;
 
 		set_intensity(ref, h, k, l, intensity);
+		if ( phases != NULL ) set_phase(phases, h, k, l, ph);
 		if ( counts != NULL ) set_count(counts, h, k, l, 1);
 
 	} while ( rval != NULL );
diff --git a/src/reflections.h b/src/reflections.h
index 052ddf9f..060a32c8 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -21,10 +21,11 @@
 
 
 extern void write_reflections(const char *filename, unsigned int *counts,
-                              double *ref, int zone_axis, UnitCell *cell,
-                              unsigned int min_counts);
+                              double *ref, double *phases, int zone_axis,
+                              UnitCell *cell, unsigned int min_counts);
 
-extern double *read_reflections(const char *filename, unsigned int *counts);
+extern double *read_reflections(const char *filename, unsigned int *counts,
+                                double *phases);
 
 extern double *ideal_intensities(double complex *sfac);
 
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 160d08dd..c4555242 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -99,7 +99,7 @@ int main(int argc, char *argv[])
 
 	cell = load_cell_from_pdb("molecule.pdb");
 	c1 = new_list_count();
-	ref1 = read_reflections(infile, c1);
+	ref1 = read_reflections(infile, c1, NULL);
 	if ( ref1 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", infile);
 		return 1;
diff --git a/src/sfac.c b/src/sfac.c
index aa7b3bf6..a8b94051 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -491,7 +491,7 @@ void free_molecule(struct molecule *mol)
 
 
 double *get_reflections(struct molecule *mol, double en, double res,
-                        unsigned int *counts)
+                        unsigned int *counts, double *phases)
 {
 	double *reflections;
 	double asx, asy, asz;
@@ -555,6 +555,7 @@ double *get_reflections(struct molecule *mol, double en, double res,
 		}
 
 		set_intensity(reflections, h, k, l, pow(cabs(F), 2.0));
+		if ( phases != NULL ) set_phase(phases, h, k, l, carg(F));
 		set_count(counts, h, k, l, 1);
 
 	}
diff --git a/src/sfac.h b/src/sfac.h
index 54104700..ea30c83e 100644
--- a/src/sfac.h
+++ b/src/sfac.h
@@ -61,7 +61,7 @@ extern double complex get_sfac(const char *n, double s, double en);
 extern struct molecule *load_molecule(const char *filename);
 extern void free_molecule(struct molecule *mol);
 extern double *get_reflections(struct molecule *mol, double en, double res,
-                               unsigned int *counts);
+                               unsigned int *counts, double *phases);
 
 
 #endif	/* SFAC_H */
diff --git a/src/utils.h b/src/utils.h
index 441fa51d..b01725dd 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -160,6 +160,11 @@ static inline double angle_between(double x1, double y1, double z1,
 #define TYPE double
 #include "list_tmp.h"
 
+/* CAs above, but for phase values */
+#define LABEL(x) x##_phase
+#define TYPE double
+#include "list_tmp.h"
+
 /* As above, but for complex structure factors */
 #define LABEL(x) x##_sfac
 #define TYPE double complex