aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/detector.c14
-rw-r--r--src/diffraction.c93
-rw-r--r--src/image.h2
3 files changed, 55 insertions, 54 deletions
diff --git a/src/detector.c b/src/detector.c
index bfc78b6e..261f9fff 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -27,22 +27,22 @@ void record_image(struct image *image)
int x, y;
double ph_per_e;
double twotheta_max, np, sa_per_pixel;
-
+
/* How many photons are scattered per electron? */
ph_per_e = PULSE_ENERGY_DENSITY * pow(THOMSON_LENGTH, 2.0)
/ image->xray_energy;
printf("%e photons are scattered per electron\n", ph_per_e);
-
+
twotheta_max = image->twotheta[0];
np = sqrt(pow(image->x_centre, 2.0) + pow(image->y_centre, 2.0));
sa_per_pixel = pow(2.0 * twotheta_max / np, 2.0);
printf("sa per pixel=%e\n", sa_per_pixel);
-
+
image->data = malloc(image->width * image->height * sizeof(uint16_t));
-
+
for ( x=0; x<image->width; x++ ) {
for ( y=0; y<image->height; y++ ) {
-
+
double counts;
double intensity;
double sa;
@@ -55,9 +55,9 @@ void record_image(struct image *image)
sa = sa_per_pixel * cos(image->twotheta[x + image->width*y]);
counts = intensity * ph_per_e * sa;
-
+
image->data[x + image->width*y] = counts;
-
+
}
}
}
diff --git a/src/diffraction.c b/src/diffraction.c
index 4f5a1d1d..9fcde28b 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -35,28 +35,28 @@ static double lattice_factor(struct threevec q, double ax, double ay, double az,
Udotq.u = (ax*q.u + ay*q.v + az*q.w)/2.0;
Udotq.v = (bx*q.u + by*q.v + bz*q.w)/2.0;
Udotq.w = (cx*q.u + cy*q.v + cz*q.w)/2.0;
-
+
if ( na > 1 ) {
f1 = sin(2.0*M_PI*(double)na*Udotq.u)
/ sin(2.0*M_PI*Udotq.u);
} else {
f1 = 1.0;
}
-
+
if ( nb > 1 ) {
f2 = sin(2.0*M_PI*(double)nb*Udotq.v)
/ sin(2.0*M_PI*Udotq.v);
} else {
f2 = 1.0;
}
-
+
if ( nc > 1 ) {
f3 = sin(2.0*M_PI*(double)nc*Udotq.w)
/ sin(2.0*M_PI*Udotq.w);
} else {
f3 = 1.0;
}
-
+
return f1 * f2 * f3;
}
@@ -134,22 +134,22 @@ static complex molecule_factor(struct molecule *mol, struct threevec q,
int i;
double F = 0.0;
double s;
-
+
s = modulus(q.u, q.v, q.w);
-
+
for ( i=0; i<mol->n_species; i++ ) {
double complex sfac;
double complex contrib = 0.0;
struct mol_species *spec;
int j;
-
+
spec = mol->species[i];
-
+
for ( j=0; j<spec->n_atoms; j++ ) {
-
+
double ph;
-
+
ph= q.u*spec->x[j] + q.v*spec->y[j] + q.w*spec->z[j];
contrib += cos(2.0*M_PI*ph) + I*sin(2.0*M_PI*ph);
}
@@ -157,7 +157,7 @@ static complex molecule_factor(struct molecule *mol, struct threevec q,
sfac = get_sfac(spec->species, en);
F += sfac * contrib * exp(-2.0 * spec->B[j] * s);
}
-
+
return F;
}
@@ -169,11 +169,11 @@ static struct molecule *load_molecule()
char line[1024];
char *rval;
int i;
-
+
mol = malloc(sizeof(struct molecule));
if ( mol == NULL ) return NULL;
mol->n_species = 0;
-
+
fh = fopen("molecule.pdb", "r");
if ( fh == NULL ) {
fprintf(stderr, "Couldn't open file\n");
@@ -181,18 +181,18 @@ static struct molecule *load_molecule()
}
do {
-
+
char el[4];
int i, r;
int done = 0;
float x, y, z, occ, B;
char *coords;
-
+
rval = fgets(line, 1023, fh);
-
+
/* Only interested in atoms */
if ( strncmp(line, "HETATM", 6) != 0 ) continue;
-
+
/* The following crimes against programming style
* were brought to you by Wizbit Enterprises, Inc. */
if ( line[76] == ' ' ) {
@@ -203,44 +203,44 @@ static struct molecule *load_molecule()
el[1] = line[77];
el[2] = '\0';
}
-
+
coords = line + 29;
r = sscanf(coords, "%f %f %f %f %f", &x, &y, &z, &occ, &B);
if ( r != 5 ) {
fprintf(stderr, "WTF?\n");
abort();
}
-
+
for ( i=0; i<mol->n_species; i++ ) {
-
+
struct mol_species *spec;
int n;
-
+
spec = mol->species[i];
-
+
if ( strcmp(spec->species, el) != 0 ) continue;
-
-
+
+
n = mol->species[i]->n_atoms;
-
+
spec->x[n] = x;
spec->y[n] = y;
spec->z[n] = z;
spec->occ[n] = occ;
spec->B[n] = B;
mol->species[i]->n_atoms++;
-
+
done = 1;
-
+
}
-
+
if ( !done ) {
-
+
/* Need to create record for this species */
struct mol_species *spec;
-
+
spec = malloc(sizeof(struct mol_species));
-
+
memcpy(spec->species, el, 4);
spec->x[0] = x;
spec->y[0] = y;
@@ -248,23 +248,23 @@ static struct molecule *load_molecule()
spec->occ[0] = occ;
spec->B[0] = B;
spec->n_atoms = 1;
-
+
mol->species[mol->n_species] = spec;
mol->n_species++;
-
+
}
-
-
+
+
} while ( rval != NULL );
-
+
fclose(fh);
-
+
printf("There are %i species\n", mol->n_species);
for ( i=0; i<mol->n_species; i++ ) {
printf("'%s': %i\n", mol->species[i]->species,
mol->species[i]->n_atoms);
}
-
+
return mol;
}
@@ -275,33 +275,34 @@ void get_diffraction(struct image *image, UnitCell *cell)
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
-
+
/* Generate the array of reciprocal space vectors in image->qvecs */
get_ewald(image);
image->molecule = load_molecule();
if ( image->molecule == NULL ) return;
-
+
cell_get_cartesian(cell, &ax, &ay, &az,
&bx, &by, &bz,
&cx, &cy, &cz);
-
+
image->sfacs = malloc(image->width * image->height * sizeof(complex));
-
+
for ( x=0; x<image->width; x++ ) {
for ( y=0; y<image->height; y++ ) {
-
+
double f_lattice;
complex f_molecule;
struct threevec q;
-
+
q = image->qvecs[x + image->width*y];
-
+
f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz);
f_molecule = molecule_factor(image->molecule, q,
image->xray_energy);
-
+
image->sfacs[x + image->width*y] = f_lattice * f_molecule;
}
+ printf("x=%i\n", x);
}
}
diff --git a/src/image.h b/src/image.h
index 7c64c6d6..ec5a20d5 100644
--- a/src/image.h
+++ b/src/image.h
@@ -61,7 +61,7 @@ struct mol_species
{
char species[4]; /* Species name */
int n_atoms; /* Number of atoms of this species */
-
+
float x[32*1024];
float y[32*1024];
float z[32*1024];
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 08:32:14 +0000