Commit fb42533c authored by Thomas White's avatar Thomas White
Browse files

Tidy up and fix

parent c72fcf4b
......@@ -129,7 +129,7 @@ static double complex water_factor(struct threevec q, double en)
/* Density of water molecules */
molecules_per_m3 = WATER_DENSITY * (AVOGADRO / WATER_MOLAR_MASS);
/* Number of water molecules per slice */
molecules_per_m = (2*rb*2*rc) * molecules_per_m3 / N_WATER_SLICES;
......@@ -168,8 +168,11 @@ void get_diffraction(struct image *image, UnitCell *cell)
/* Generate the array of reciprocal space vectors in image->qvecs */
get_ewald(image);
image->molecule = load_molecule();
if ( image->molecule == NULL ) return;
if ( image->molecule == NULL ) {
image->molecule = load_molecule();
if ( image->molecule == NULL ) return;
}
cell_get_cartesian(cell, &ax, &ay, &az,
&bx, &by, &bz,
......
......@@ -55,32 +55,32 @@ void get_ewald(struct image *image)
{
int x, y;
double k; /* Wavenumber */
k = 1/image->lambda;
image->qvecs = malloc(image->width * image->height
* sizeof(struct threevec));
image->twotheta = malloc(image->width * image->height
* sizeof(double));
for ( x=0; x<image->width; x++ ) {
for ( y=0; y<image->height; y++ ) {
double rx, ry, r;
double twothetax, twothetay, twotheta;
double qx, qy, qz;
struct threevec q;
/* Calculate q vectors for Ewald sphere */
rx = ((double)x - image->x_centre) / image->resolution;
ry = ((double)y - image->y_centre) / image->resolution;
r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
twothetax = atan2(rx, image->camera_len);
twothetay = atan2(ry, image->camera_len);
twothetay = atan2(ry, image->camera_len);
twotheta = atan2(r, image->camera_len);
qx = k * sin(twothetax);
qy = k * sin(twothetay);
qz = k - k * cos(twotheta);
......@@ -89,11 +89,8 @@ void get_ewald(struct image *image)
image->qvecs[x + image->width*y] = quat_rot(q,
image->orientation);
image->qvecs[x + image->width*y].u = qx;
image->qvecs[x + image->width*y].v = qy;
image->qvecs[x + image->width*y].w = qz;
image->twotheta[x + image->width*y] = twotheta;
}
}
}
......@@ -67,6 +67,21 @@ int main(int argc, char *argv[])
deg2rad(90.0),
deg2rad(120.0));
/* Define image parameters */
image.width = 512;
image.height = 512;
image.fmode = FORMULATION_CLEN;
image.x_centre = 255.5;
image.y_centre = 255.5;
image.camera_len = 0.2; /* 20 cm */
image.resolution = 5120; /* 512 pixels in 10 cm */
image.xray_energy = eV_to_J(2.0e3); /* 2 keV energy */
image.lambda = ph_en_to_lambda(image.xray_energy); /* Wavelength */
image.molecule = NULL;
/* Splurge a few useful numbers */
printf("Wavelength is %f nm\n", image.lambda/1.0e-9);
again:
/* Read quaternion from stdin */
......@@ -102,23 +117,12 @@ again:
} while ( !done );
/* Define image parameters */
image.width = 512;
image.height = 512;
image.fmode = FORMULATION_CLEN;
image.x_centre = 255.5;
image.y_centre = 255.5;
image.camera_len = 0.2; /* 20 cm */
image.resolution = 5120; /* 512 pixels in 10 cm */
image.xray_energy = eV_to_J(2.0e3); /* 2 keV energy */
image.lambda = ph_en_to_lambda(image.xray_energy); /* Wavelength */
/* Ensure no residual information */
image.qvecs = NULL;
image.sfacs = NULL;
image.data = NULL;
image.twotheta = NULL;
/* Splurge a few useful numbers */
printf("Wavelength is %f nm\n", image.lambda/1.0e-9);
image.hdr = NULL;
get_diffraction(&image, cell);
record_image(&image);
......@@ -128,6 +132,14 @@ again:
/* Write the output file */
hdf5_write(filename, image.data, image.width, image.height);
printf("Mock write: %i\n", number-1);
/* Clean up */
free(image.data);
free(image.qvecs);
free(image.hdr);
free(image.sfacs);
free(image.twotheta);
/* Do it all again */
goto again;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment