aboutsummaryrefslogtreecommitdiff
path: root/src/diffraction.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2009-11-17 15:43:05 +0100
committerThomas White <taw@physics.org>2009-11-17 15:43:05 +0100
commit50ffbc353574996546e51a61ec4566f461769fc3 (patch)
treee5df514215814f4f6b1e0f92195b9e209e35fe90 /src/diffraction.c
parent04c4c604f03d7c9729684b7b226a9ef695ece21d (diff)
Move scattering factor stuff to new file
Diffstat (limited to 'src/diffraction.c')
-rw-r--r--src/diffraction.c142
1 files changed, 1 insertions, 141 deletions
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,