aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/cell.c30
-rw-r--r--src/cell.h3
-rw-r--r--src/diffraction.c60
-rw-r--r--src/sfac.c62
-rw-r--r--src/sfac.h1
5 files changed, 98 insertions, 58 deletions
diff --git a/src/cell.c b/src/cell.c
index 5e8b1ce6..6fc29e56 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -206,3 +206,33 @@ void cell_get_reciprocal(UnitCell *cell,
*bsz = gsl_matrix_get(inv, 2, 1);
*csz = gsl_matrix_get(inv, 2, 2);
}
+
+
+double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
+{
+ const double a = cell->a;
+ const double b = cell->b;
+ const double c = cell->c;
+ const double alpha = cell->alpha;
+ const double beta = cell->beta;
+ const double gamma = cell->gamma;
+
+ const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha)
+ - cos(beta)*cos(beta)
+ - cos(gamma)*cos(gamma)
+ + 2*cos(alpha)*cos(beta)*cos(gamma) );
+
+ const double S11 = b*b*c*c*sin(alpha)*sin(alpha);
+ const double S22 = a*a*c*c*sin(beta)*sin(beta);
+ const double S33 = a*a*b*b*sin(gamma)*sin(gamma);
+ const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma));
+ const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha));
+ const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta));
+
+ const double brackets = S11*h*h + S22*k*k + S33*l*l
+ + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l;
+ const double oneoverdsq = brackets / Vsq;
+ const double oneoverd = sqrt(oneoverdsq);
+
+ return oneoverd / 2;
+}
diff --git a/src/cell.h b/src/cell.h
index b2203561..860cce05 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -62,4 +62,7 @@ extern void cell_get_reciprocal(UnitCell *cell,
double *bsx, double *bsy, double *bsz,
double *csx, double *csy, double *csz);
+extern double resolution(UnitCell *cell,
+ signed int h, signed int k, signed int l);
+
#endif /* CELL_H */
diff --git a/src/diffraction.c b/src/diffraction.c
index e112e7ce..9a9e31b6 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -84,6 +84,7 @@ static double complex molecule_factor(struct molecule *mol, struct threevec q,
{
double hd, kd, ld;
signed int h, k, l;
+ double complex r;
hd = q.u * ax + q.v * ay + q.w * az;
kd = q.u * bx + q.v * by + q.w * bz;
@@ -92,7 +93,9 @@ static double complex molecule_factor(struct molecule *mol, struct threevec q,
k = (signed int)nearbyint(kd);
l = (signed int)nearbyint(ld);
- return get_integral(mol->reflections, h, k, l);
+ r = get_integral(mol->reflections, h, k, l);
+
+ return r;
}
@@ -102,51 +105,6 @@ static double complex water_factor(struct threevec q, double en)
}
-static void get_reflections_cached(struct molecule *mol, double en)
-{
- char s[1024];
- FILE *fh;
- size_t r;
-
- snprintf(s, 1023, "reflections-%ieV.cache", (int)J_to_eV(en));
- fh = fopen(s, "rb");
- if ( fh == NULL ) {
- printf("No cache file found (looked for %s)\n", s);
- goto calc;
- }
-
- mol->reflections = reflist_new();
- r = fread(mol->reflections, sizeof(double complex), IDIM*IDIM*IDIM, fh);
- if ( r < IDIM*IDIM*IDIM ) {
- printf("Found cache file (%s), but failed to read.\n", s);
- goto calc;
- }
-
- printf("Read structure factors (at Bragg positions) from %s\n", s);
- return;
-
-calc:
- printf("Calculating structure factors at Bragg positions...\n");
- mol->reflections = get_reflections(mol, en);
- fh = fopen(s, "wb");
- if ( fh == NULL ) {
- printf("Failed to write cache file (%s)\n", s);
- return;
- }
-
- r = fwrite(mol->reflections, sizeof(double complex),
- IDIM*IDIM*IDIM, fh);
- if ( r < IDIM*IDIM*IDIM ) {
- printf("Failed to write cache file (%s)\n", s);
- return;
- }
- fclose(fh);
-
- printf("Successfully saved structure factors at Bragg positions to"
- " file %s\n", s);
-}
-
-
void get_diffraction(struct image *image)
{
int x, y;
@@ -189,14 +147,8 @@ void get_diffraction(struct image *image)
f_molecule = molecule_factor(image->molecule, q,
ax,ay,az,bx,by,bz,cx,cy,cz);
- /* Nasty approximation follows */
- if ( y == image->height/2 ) {
- f_water = water_factor(q, image->xray_energy);
- } else {
- f_water = 0.0;
- }
-
- val = (f_molecule) + f_water;
+ f_water = water_factor(q, image->xray_energy);
+ val = (f_molecule * f_lattice) + f_water;
image->sfacs[x + image->width*y] = val;
}
diff --git a/src/sfac.c b/src/sfac.c
index 7257bf01..7bf30f82 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -175,7 +175,7 @@ static double get_waas_kirf(const char *n, double s)
fclose(fh);
- s2 = pow(s/1e10, 2.0); /* s2 is s squared in Angstroms squared */
+ s2 = pow(s/1e10, 2.0); /* s2 is s squared in Angstroms^-2 */
f = c + a1*exp(-b1*s2) + a2*exp(-b2*s2) + a3*exp(-b3*s2)
+ a4*exp(-b4*s2) + a5*exp(-b5*s2);
@@ -390,6 +390,8 @@ double complex *get_reflections(struct molecule *mol, double en)
double bsx, bsy, bsz;
double csx, csy, csz;
signed int h, k, l;
+ //double tscat = 0.0;
+ //double F00;
cell_get_reciprocal(mol->cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
@@ -401,10 +403,12 @@ double complex *get_reflections(struct molecule *mol, double en)
for ( k=-INDMAX; k<=INDMAX; k++ ) {
for ( l=-INDMAX; l<=INDMAX; l++ ) {
- double complex F;
+ double complex F = 0.0;
int i;
double s;
+ s = resolution(mol->cell, h, k, l);
+
/* Atoms are grouped by species for faster calculation */
for ( i=0; i<mol->n_species; i++ ) {
@@ -426,8 +430,7 @@ double complex *get_reflections(struct molecule *mol, double en)
ph = u*spec->x[j] + v*spec->y[j] + w*spec->z[j];
/* Conversion from revolutions to radians */
- contrib += cos(-2.0*M_PI*ph)
- + I*sin(-2.0*M_PI*ph);
+ contrib += cexp(-2.0*M_PI*I*ph);
}
@@ -437,12 +440,63 @@ double complex *get_reflections(struct molecule *mol, double en)
}
integrate_reflection(reflections, h, k, l, F);
+ //if ( (h!=0) || (k!=0) || (l!=0) ) {
+ // tscat += cabs(F);
+ //} else {
+ // F00 = cabs(F);
+ //}
}
}
progress_bar((h+INDMAX+1), 2*INDMAX);
}
printf("\n");
+ //printf("Total scattered = %f, F000 = %f\n", tscat, F00);
return reflections;
}
+
+
+void get_reflections_cached(struct molecule *mol, double en)
+{
+ char s[1024];
+ FILE *fh;
+ size_t r;
+
+ snprintf(s, 1023, "reflections-%ieV.cache", (int)J_to_eV(en));
+ fh = fopen(s, "rb");
+ if ( fh == NULL ) {
+ printf("No cache file found (looked for %s)\n", s);
+ goto calc;
+ }
+
+ mol->reflections = reflist_new();
+ r = fread(mol->reflections, sizeof(double complex), IDIM*IDIM*IDIM, fh);
+ if ( r < IDIM*IDIM*IDIM ) {
+ printf("Found cache file (%s), but failed to read.\n", s);
+ goto calc;
+ }
+
+ printf("Read structure factors (at Bragg positions) from %s\n", s);
+ return;
+
+calc:
+ printf("Calculating structure factors at Bragg positions...\n");
+ mol->reflections = get_reflections(mol, en);
+ fh = fopen(s, "wb");
+ if ( fh == NULL ) {
+ printf("Failed to write cache file (%s)\n", s);
+ return;
+ }
+
+ r = fwrite(mol->reflections, sizeof(double complex),
+ IDIM*IDIM*IDIM, fh);
+ if ( r < IDIM*IDIM*IDIM ) {
+ printf("Failed to write cache file (%s)\n", s);
+ return;
+ }
+ fclose(fh);
+
+ printf("Successfully saved structure factors at Bragg positions to"
+ " file %s\n", s);
+}
diff --git a/src/sfac.h b/src/sfac.h
index fc4b2bd9..131a05e8 100644
--- a/src/sfac.h
+++ b/src/sfac.h
@@ -55,5 +55,6 @@ struct molecule
extern double complex get_sfac(const char *n, double s, double en);
extern struct molecule *load_molecule(void);
extern double complex *get_reflections(struct molecule *mol, double en);
+extern void get_reflections_cached(struct molecule *mol, double en);
#endif /* SFAC_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 21:13:37 +0000