diff options
author | Thomas White <taw@physics.org> | 2010-12-03 15:04:02 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2012-02-22 15:27:07 +0100 |
commit | 27d14f0b1391bbc6d667eb9c82e78c3b02052e5a (patch) | |
tree | 136e6a3ac7ec8eaf4f36a9a88913838c090125c3 /src/diffraction.c | |
parent | 796feb582e9dff7511c411f0e97dcdf382a6f85d (diff) |
Use symmetry when simulating (on the CPU only)
Diffstat (limited to 'src/diffraction.c')
-rw-r--r-- | src/diffraction.c | 122 |
1 files changed, 96 insertions, 26 deletions
diff --git a/src/diffraction.c b/src/diffraction.c index 93ac7f8e..aa06c44b 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -23,6 +23,7 @@ #include "diffraction.h" #include "sfac.h" #include "beam-parameters.h" +#include "symmetry.h" #define SAMPLING (4) @@ -69,8 +70,63 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az, } -static double interpolate_linear(const double *ref, - float hd, signed int k, signed int l) +static double sym_lookup_intensity(const double *intensities, + const unsigned char *flags, const char *sym, + signed int h, signed int k, signed int l) +{ + int i; + double ret = 0.0; + + for ( i=0; i<num_general_equivs(sym); i++ ) { + + signed int he; + signed int ke; + signed int le; + double f, val; + + get_general_equiv(h, k, l, &he, &ke, &le, sym, i); + + f = (double)lookup_flag(flags, he, ke, le); + val = lookup_intensity(intensities, he, ke, le); + + ret += f*val; + + } + + return ret; +} + + +static double sym_lookup_phase(const double *phases, + const unsigned char *flags, const char *sym, + signed int h, signed int k, signed int l) +{ + int i; + double ret = 0.0; + + for ( i=0; i<num_general_equivs(sym); i++ ) { + + signed int he; + signed int ke; + signed int le; + double f, val; + + get_general_equiv(h, k, l, &he, &ke, &le, sym, i); + + f = (double)lookup_flag(flags, he, ke, le); + val = lookup_phase(phases, he, ke, le); + + ret += f*val; + + } + + return ret; +} + + +static double interpolate_linear(const double *ref, const unsigned char *flags, + const char *sym, float hd, + signed int k, signed int l) { signed int h; double val1, val2; @@ -81,8 +137,8 @@ static double interpolate_linear(const double *ref, f = hd - (float)h; assert(f >= 0.0); - val1 = lookup_intensity(ref, h, k, l); - val2 = lookup_intensity(ref, h+1, k, l); + val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); + val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); val1 = val1; val2 = val2; @@ -92,6 +148,7 @@ static double interpolate_linear(const double *ref, static double interpolate_bilinear(const double *ref, + const unsigned char *flags, const char *sym, float hd, float kd, signed int l) { signed int k; @@ -103,14 +160,15 @@ static double interpolate_bilinear(const double *ref, f = kd - (float)k; assert(f >= 0.0); - val1 = interpolate_linear(ref, hd, k, l); - val2 = interpolate_linear(ref, hd, k+1, l); + val1 = interpolate_linear(ref, flags, sym, hd, k, l); + val2 = interpolate_linear(ref, flags, sym, hd, k+1, l); return (1.0-f)*val1 + f*val2; } static double interpolate_intensity(const double *ref, + const unsigned char *flags, const char *sym, float hd, float kd, float ld) { signed int l; @@ -122,8 +180,8 @@ static double interpolate_intensity(const double *ref, f = ld - (float)l; assert(f >= 0.0); - val1 = interpolate_bilinear(ref, hd, kd, l); - val2 = interpolate_bilinear(ref, hd, kd, l+1); + val1 = interpolate_bilinear(ref, flags, sym, hd, kd, l); + val2 = interpolate_bilinear(ref, flags, sym, hd, kd, l+1); return (1.0-f)*val1 + f*val2; } @@ -131,6 +189,8 @@ static double interpolate_intensity(const double *ref, static double complex interpolate_phased_linear(const double *ref, const double *phases, + const unsigned char *flags, + const char *sym, float hd, signed int k, signed int l) { @@ -146,10 +206,10 @@ static double complex interpolate_phased_linear(const double *ref, f = hd - (float)h; assert(f >= 0.0); - val1 = lookup_intensity(ref, h, k, l); - val2 = lookup_intensity(ref, h+1, k, l); - ph1 = lookup_phase(phases, h, k, l); - ph2 = lookup_phase(phases, h+1, k, l); + val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); + val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); + ph1 = sym_lookup_phase(phases, flags, sym, h, k, l); + ph2 = sym_lookup_phase(phases, flags, sym, h+1, k, l); val1 = val1; val2 = val2; @@ -169,6 +229,8 @@ static double complex interpolate_phased_linear(const double *ref, static double complex interpolate_phased_bilinear(const double *ref, const double *phases, + const unsigned char *flags, + const char *sym, float hd, float kd, signed int l) { @@ -181,8 +243,8 @@ static double complex interpolate_phased_bilinear(const double *ref, f = kd - (float)k; assert(f >= 0.0); - val1 = interpolate_phased_linear(ref, phases, hd, k, l); - val2 = interpolate_phased_linear(ref, phases, hd, k+1, l); + val1 = interpolate_phased_linear(ref, phases, flags, sym, hd, k, l); + val2 = interpolate_phased_linear(ref, phases, flags, sym, hd, k+1, l); return (1.0-f)*val1 + f*val2; } @@ -190,6 +252,8 @@ static double complex interpolate_phased_bilinear(const double *ref, static double interpolate_phased_intensity(const double *ref, const double *phases, + const unsigned char *flags, + const char *sym, float hd, float kd, float ld) { signed int l; @@ -201,20 +265,22 @@ static double interpolate_phased_intensity(const double *ref, f = ld - (float)l; assert(f >= 0.0); - val1 = interpolate_phased_bilinear(ref, phases, hd, kd, l); - val2 = interpolate_phased_bilinear(ref, phases, hd, kd, l+1); + val1 = interpolate_phased_bilinear(ref, phases, flags, sym, + hd, kd, l); + val2 = interpolate_phased_bilinear(ref, phases, flags, sym, + 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 double *phases, - struct rvec q, +static double molecule_factor(const double *intensities, const double *phases, + const unsigned char *flags, struct rvec q, double ax, double ay, double az, double bx, double by, double bz, double cx, double cy, double cz, - GradientMethod m) + GradientMethod m, const char *sym) { float hd, kd, ld; signed int h, k, l; @@ -224,6 +290,9 @@ static double molecule_factor(const double *intensities,const double *phases, kd = q.u * bx + q.v * by + q.w * bz; ld = q.u * cx + q.v * cy + q.w * cz; + /* No flags -> flat intensity distribution */ + if ( flags == NULL ) return 1.0e5; + switch ( m ) { case GRADIENT_MOSAIC : h = (signed int)rint(hd); @@ -232,14 +301,14 @@ static double molecule_factor(const double *intensities,const double *phases, if ( abs(h) > INDMAX ) r = 0.0; else if ( abs(k) > INDMAX ) r = 0.0; else if ( abs(l) > INDMAX ) r = 0.0; - else r = lookup_intensity(intensities, h, k, l); + else r = sym_lookup_intensity(intensities, flags, sym, h, k, l); break; case GRADIENT_INTERPOLATE : - r = interpolate_intensity(intensities, hd, kd, ld); + r = interpolate_intensity(intensities, flags, sym, hd, kd, ld); break; case GRADIENT_PHASED : - r = interpolate_phased_intensity(intensities, phases, - hd, kd, ld); + r = interpolate_phased_intensity(intensities, phases, flags, + sym, hd, kd, ld); break; default: ERROR("This gradient method not implemented yet.\n"); @@ -298,7 +367,8 @@ double water_diffraction(struct rvec q, double en, void get_diffraction(struct image *image, int na, int nb, int nc, const double *intensities, const double *phases, - UnitCell *cell, int do_water, GradientMethod m) + const unsigned char *flags, UnitCell *cell, int do_water, + GradientMethod m, const char *sym) { unsigned int xs, ys; double ax, ay, az; @@ -353,10 +423,10 @@ void get_diffraction(struct image *image, int na, int nb, int nc, I_molecule = 1.0e10; } else { I_molecule = molecule_factor(intensities, - phases, q, + phases, flags, q, ax,ay,az, bx,by,bz,cx,cy,cz, - m); + m, sym); } I_lattice = pow(f_lattice, 2.0); |