From 50ffbc353574996546e51a61ec4566f461769fc3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 17 Nov 2009 15:43:05 +0100 Subject: Move scattering factor stuff to new file --- Makefile.am | 3 +- src/Makefile.am | 2 +- src/diffraction.c | 142 +--------------------------------------------- src/sfac.c | 164 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/sfac.h | 23 ++++++++ 5 files changed, 191 insertions(+), 143 deletions(-) create mode 100644 src/sfac.c create mode 100644 src/sfac.h diff --git a/Makefile.am b/Makefile.am index 9087ee19..955cc8f0 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,3 +1,4 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h src/relrod.h \ - src/utils.h src/diffraction.h src/detector.h src/ewald.h + src/utils.h src/diffraction.h src/detector.h src/ewald.h \ + src/sfac.h SUBDIRS = src diff --git a/src/Makefile.am b/src/Makefile.am index c224a2b0..1220b7f2 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -3,5 +3,5 @@ bin_PROGRAMS = pattern_sim AM_CFLAGS = -Wall -g @CFLAGS@ pattern_sim_SOURCES = main.c diffraction.c utils.c image.c cell.c hdf5-file.c \ - ewald.c detector.c + ewald.c detector.c sfac.c pattern_sim_LDADD = @LIBS@ diff --git a/src/diffraction.c b/src/diffraction.c index e02d80f4..b5c07a71 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -21,6 +21,7 @@ #include "cell.h" #include "ewald.h" #include "diffraction.h" +#include "sfac.h" static double lattice_factor(struct threevec q, double ax, double ay, double az, @@ -59,147 +60,6 @@ static double lattice_factor(struct threevec q, double ax, double ay, double az, } -/* Look up f1 and f2 for this atom at this energy (in J/photon) */ -static double complex get_f1f2(const char *n, double en) -{ - FILE *fh; - char filename[64]; - char line[1024]; - char *rval; - float last_E, last_f1, last_f2; - - snprintf(filename, 63, "scattering-factors/%s.nff", n); - fh = fopen(filename, "r"); - if ( fh == NULL ) { - fprintf(stderr, "Couldn't open file '%s'\n", filename); - return 0.0; - } - - en = J_to_eV(en); - - /* Discard first line */ - fgets(line, 1023, fh); - - last_E = 0.0; - last_f1 = 0.0; - last_f2 = 0.0; - do { - - int r; - float E, f1, f2; - - rval = fgets(line, 1023, fh); - - r = sscanf(line, "%f %f %f", &E, &f1, &f2); - if ( r != 3 ) { - fprintf(stderr, "WTF?\n"); - abort(); - } - - /* Find the first energy greater than the required value */ - if ( E < en ) { - /* Store old values ready for interpolation*/ - last_E = E; - last_f1 = f1; - last_f2 = f2; - } else { - - /* Perform (linear) interpolation */ - float f; - float actual_f1, actual_f2; - - f = (en - last_E) / (E - last_E); - - actual_f1 = last_f1 + f * (f1 - last_f1); - actual_f2 = last_f2 + f * (f2 - last_f2); - - fclose(fh); - return actual_f1 + I*actual_f2; - - } - - } while ( rval != NULL ); - - fclose(fh); - - fprintf(stderr, "Couldn't find scattering factors for '%s' at %f eV!\n", - n, en); - - return 0.0; -} - - -/* s = sin(theta)/lambda */ -static double get_waas_kirf(const char *n, double s) -{ - FILE *fh; - char *rval; - double f; - float a1, a2, a3, a4, a5, c, b1, b2, b3, b4, b5; - double s2; - - fh = fopen("scattering-factors/f0_WaasKirf.dat", "r"); - if ( fh == NULL ) { - fprintf(stderr, "Couldn't open f0_WaasKirf.dat\n"); - return 0.0; - } - - do { - - int r; - char line[1024]; - char sp[1024]; - int Z; - - rval = fgets(line, 1023, fh); - - if ( (line[0] != '#') || (line[1] != 'S') ) continue; - - r = sscanf(line, "#S %i %s", &Z, sp); - if ( (r != 2) || (strcmp(sp, n) != 0) ) continue; - - /* Skip two lines */ - fgets(line, 1023, fh); - fgets(line, 1023, fh); - - /* Read scattering coefficients */ - rval = fgets(line, 1023, fh); - r = sscanf(line, " %f %f %f %f %f %f %f %f %f %f %f", - &a1, &a2, &a3, &a4, &a5, &c, &b1, &b2, &b3, &b4, &b5); - if ( r != 11 ) { - fprintf(stderr, "Couldn't read scattering factors\n"); - return 0.0; - } - - break; - - } while ( rval != NULL ); - - fclose(fh); - - s2 = pow(s/1e10, 2.0); /* s2 is s squared in Angstroms squared */ - f = c + a1*exp(-b1*s2) + a2*exp(-b2*s2) + a3*exp(-b3*s2) - + a4*exp(-b4*s2) + a5*exp(-b5*s2); - - return f; -} - - -/* Get complex scattering factors for element 'n' at energy 'en' (J/photon), - * at resolution 's' = sin(theta)/lambda (in m^-1) */ -static double complex get_sfac(const char *n, double s, double en) -{ - double complex f1f2; - double fq, fq0; - - f1f2 = get_f1f2(n, en); - fq = get_waas_kirf(n, s); - fq0 = get_waas_kirf(n, 0.0); - - return fq - fq0 + f1f2; -} - - /* Return structure factor for molecule 'mol' at energy en' (J/photon) at * scattering vector 'q' */ static double complex molecule_factor(struct molecule *mol, struct threevec q, diff --git a/src/sfac.c b/src/sfac.c new file mode 100644 index 00000000..4c3dabd3 --- /dev/null +++ b/src/sfac.c @@ -0,0 +1,164 @@ +/* + * sfac.c + * + * Scattering factors + * + * (c) 2007-2009 Thomas White + * + * pattern_sim - Simulate diffraction patterns from small crystals + * + */ + + +#include +#include +#include +#include +#include + +#include "utils.h" +#include "sfac.h" + + +/* Look up f1 and f2 for this atom at this energy (in J/photon) */ +static double complex get_f1f2(const char *n, double en) +{ + FILE *fh; + char filename[64]; + char line[1024]; + char *rval; + double last_E, last_f1, last_f2; + + snprintf(filename, 63, "scattering-factors/%s.nff", n); + fh = fopen(filename, "r"); + if ( fh == NULL ) { + fprintf(stderr, "Couldn't open file '%s'\n", filename); + return 0.0; + } + + en = J_to_eV(en); + + /* Discard first line */ + fgets(line, 1023, fh); + + last_E = 0.0; + last_f1 = 0.0; + last_f2 = 0.0; + do { + + int r; + double E, f1, f2; + float E_f, f1_f, f2_f; + + rval = fgets(line, 1023, fh); + + r = sscanf(line, "%f %f %f", &E_f, &f1_f, &f2_f); + if ( r != 3 ) { + fprintf(stderr, "WTF?\n"); + abort(); + } + /* Promote to double precision */ + E = E_f; f1 = f1_f; f2 = f2_f; + + /* Find the first energy greater than the required value */ + if ( E < en ) { + /* Store old values ready for interpolation*/ + last_E = E; + last_f1 = f1; + last_f2 = f2; + } else { + + /* Perform (linear) interpolation */ + double f; + double actual_f1, actual_f2; + + f = (en - last_E) / (E - last_E); + + actual_f1 = last_f1 + f * (f1 - last_f1); + actual_f2 = last_f2 + f * (f2 - last_f2); + + fclose(fh); + return actual_f1 + I*actual_f2; + + } + + } while ( rval != NULL ); + + fclose(fh); + + fprintf(stderr, "Couldn't find scattering factors for '%s' at %f eV!\n", + n, en); + + return 0.0; +} + + +/* s = sin(theta)/lambda */ +static double get_waas_kirf(const char *n, double s) +{ + FILE *fh; + char *rval; + double f; + float a1, a2, a3, a4, a5, c, b1, b2, b3, b4, b5; + double s2; + + fh = fopen("scattering-factors/f0_WaasKirf.dat", "r"); + if ( fh == NULL ) { + fprintf(stderr, "Couldn't open f0_WaasKirf.dat\n"); + return 0.0; + } + + do { + + int r; + char line[1024]; + char sp[1024]; + int Z; + + rval = fgets(line, 1023, fh); + + if ( (line[0] != '#') || (line[1] != 'S') ) continue; + + r = sscanf(line, "#S %i %s", &Z, sp); + if ( (r != 2) || (strcmp(sp, n) != 0) ) continue; + + /* Skip two lines */ + fgets(line, 1023, fh); + fgets(line, 1023, fh); + + /* Read scattering coefficients */ + rval = fgets(line, 1023, fh); + r = sscanf(line, " %f %f %f %f %f %f %f %f %f %f %f", + &a1, &a2, &a3, &a4, &a5, &c, &b1, &b2, &b3, &b4, &b5); + if ( r != 11 ) { + fprintf(stderr, "Couldn't read scattering factors\n"); + return 0.0; + } + + break; + + } while ( rval != NULL ); + + fclose(fh); + + s2 = pow(s/1e10, 2.0); /* s2 is s squared in Angstroms squared */ + f = c + a1*exp(-b1*s2) + a2*exp(-b2*s2) + a3*exp(-b3*s2) + + a4*exp(-b4*s2) + a5*exp(-b5*s2); + + return f; +} + + +/* Get complex scattering factors for element 'n' at energy 'en' (J/photon), + * at resolution 's' = sin(theta)/lambda (in m^-1) */ +double complex get_sfac(const char *n, double s, double en) +{ + double complex f1f2; + double fq, fq0; + + f1f2 = get_f1f2(n, en); + fq = get_waas_kirf(n, s); + fq0 = get_waas_kirf(n, 0.0); + + return fq - fq0 + f1f2; +} diff --git a/src/sfac.h b/src/sfac.h new file mode 100644 index 00000000..afa84f84 --- /dev/null +++ b/src/sfac.h @@ -0,0 +1,23 @@ +/* + * sfac.h + * + * Scattering factors + * + * (c) 2007-2009 Thomas White + * + * pattern_sim - Simulate diffraction patterns from small crystals + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef SFAC_H +#define SFAC_H + +#include + +extern double complex get_sfac(const char *n, double s, double en); + +#endif /* SFAC_H */ -- cgit v1.2.3