Commit 8ead809d authored by Thomas White's avatar Thomas White
Browse files

pattern_sim: Implement phased gradients

parent 509f08dc
......@@ -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;
......
......@@ -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);
......
......@@ -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);
......
......@@ -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;
}
......@@ -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;
......
......@@ -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");
......
......@@ -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);
......
......@@ -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 );
......
......@@ -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);
......
......@@ -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;
......
......@@ -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);
}
......
......@@ -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 */
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment