diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/diffraction.c | 9 | ||||
| -rw-r--r-- | src/ewald.c | 23 | ||||
| -rw-r--r-- | src/main.c | 40 |
3 files changed, 42 insertions, 30 deletions
diff --git a/src/diffraction.c b/src/diffraction.c index 87f2a5f7..6bb48332 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -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, diff --git a/src/ewald.c b/src/ewald.c index 7a620c42..d1f47da7 100644 --- a/src/ewald.c +++ b/src/ewald.c @@ -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; } |