aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/compare_hkl.c6
-rw-r--r--src/diffraction.c108
-rw-r--r--src/diffraction.h4
-rw-r--r--src/get_hkl.c10
-rw-r--r--src/indexamajig.c5
-rw-r--r--src/pattern_sim.c7
-rw-r--r--src/process_hkl.c9
-rw-r--r--src/reflections.c28
-rw-r--r--src/reflections.h7
-rw-r--r--src/render_hkl.c2
-rw-r--r--src/sfac.c3
-rw-r--r--src/sfac.h2
-rw-r--r--src/utils.h5
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
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-02 07:53:08 +0000