aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am3
-rw-r--r--src/diffraction.c89
-rw-r--r--src/ewald.c57
-rw-r--r--src/ewald.h23
-rw-r--r--src/hdf5-file.c8
-rw-r--r--src/hdf5-file.h2
-rw-r--r--src/image.c1
-rw-r--r--src/image.h12
-rw-r--r--src/main.c31
-rw-r--r--src/utils.c34
-rw-r--r--src/utils.h5
11 files changed, 205 insertions, 60 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index eca48660..f5ce4e6c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -2,5 +2,6 @@ bin_PROGRAMS = pattern_sim
AM_CFLAGS = -Wall -g @CFLAGS@
-pattern_sim_SOURCES = main.c diffraction.c utils.c image.c cell.c hdf5-file.c
+pattern_sim_SOURCES = main.c diffraction.c utils.c image.c cell.c hdf5-file.c \
+ ewald.c
pattern_sim_LDADD = @LIBS@
diff --git a/src/diffraction.c b/src/diffraction.c
index 93c6a22c..b2da4622 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -17,43 +17,64 @@
#include "image.h"
#include "utils.h"
#include "cell.h"
-
-
-static void mapping_rotate(double x, double y, double z,
- double *ddx, double *ddy, double *ddz,
- double omega, double tilt)
-{
- double nx, ny, nz;
- double x_temp, y_temp, z_temp;
-
- /* First: rotate image clockwise until tilt axis is aligned
- * horizontally. */
- nx = x*cos(omega) + y*sin(omega);
- ny = -x*sin(omega) + y*cos(omega);
- nz = z;
-
- /* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the
- * "wrong" way. This is because the crystal is rotated in the
- * experiment, not the Ewald sphere. */
- x_temp = nx; y_temp = ny; z_temp = nz;
- nx = x_temp;
- ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
- nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
-
- /* Finally, reverse the omega rotation to restore the location of the
- * image in 3D space */
- x_temp = nx; y_temp = ny; z_temp = nz;
- nx = x_temp*cos(-omega) + y_temp*sin(-omega);
- ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
- nz = z_temp;
-
- *ddx = nx;
- *ddy = ny;
- *ddz = nz;
-}
+#include "ewald.h"
void get_diffraction(struct image *image, UnitCell *cell)
{
+ int x, y;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ int na = 64;
+ int nb = 64;
+ int nc = 64;
+
+ /* Generate the array of reciprocal space vectors in image->qvecs */
+ get_ewald(image);
+
+ cell_get_cartesian(cell, &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz);
+
+ image->sfacs = malloc(image->width * image->height * sizeof(double));
+
+ for ( x=0; x<image->width; x++ ) {
+ for ( y=0; y<image->height; y++ ) {
+
+ struct threevec q;
+ struct threevec Udotq;
+ double f1, f2, f3;
+
+ q = image->qvecs[x + image->width*y];
+
+ 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;
+ }
+
+ image->sfacs[x + image->width*y] = f1 * f2 * f3;
+ }
+ }
}
diff --git a/src/ewald.c b/src/ewald.c
new file mode 100644
index 00000000..9d8a0450
--- /dev/null
+++ b/src/ewald.c
@@ -0,0 +1,57 @@
+/*
+ * ewald.c
+ *
+ * Calculate q-vector arrays
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * pattern_sim - Simulate diffraction patterns from small crystals
+ *
+ */
+
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "image.h"
+#include "utils.h"
+#include "cell.h"
+
+
+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));
+
+ 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;
+
+ 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);
+ twotheta = atan2(r, image->camera_len);
+
+ qx = k * sin(twothetax);
+ qy = k * sin(twothetay);
+ qz = k - k * cos(twotheta);
+
+ image->qvecs[x + image->width*y].u = qx;
+ image->qvecs[x + image->width*y].v = qy;
+ image->qvecs[x + image->width*y].w = qz;
+
+ }
+ }
+}
diff --git a/src/ewald.h b/src/ewald.h
new file mode 100644
index 00000000..ffa6e162
--- /dev/null
+++ b/src/ewald.h
@@ -0,0 +1,23 @@
+/*
+ * ewald.h
+ *
+ * Calculate q-vector arrays
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * pattern_sim - Simulate diffraction patterns from small crystals
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef EWALD_H
+#define EWALD_H
+
+#include "image.h"
+
+extern void get_ewald(struct image *image);
+
+#endif /* EWALD_H */
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index ca1d5450..6eb84aeb 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -22,7 +22,7 @@
#include "image.h"
-int hdf5_write(const char *filename, const uint16_t *data,
+int hdf5_write(const char *filename, const double *data,
int width, int height)
{
hid_t fh, gh, sh, dh; /* File, group, dataspace and data handles */
@@ -49,7 +49,7 @@ int hdf5_write(const char *filename, const uint16_t *data,
max_size[1] = height;
sh = H5Screate_simple(2, size, max_size);
- dh = H5Dcreate(gh, "data", H5T_NATIVE_UINT16, sh,
+ dh = H5Dcreate(gh, "data", H5T_NATIVE_FLOAT, sh,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if ( dh < 0 ) {
fprintf(stderr, "Couldn't create dataset\n");
@@ -63,7 +63,7 @@ int hdf5_write(const char *filename, const uint16_t *data,
(int)size[1], (int)size[0],
(int)max_size[1], (int)max_size[0]);
- r = H5Dwrite(dh, H5T_NATIVE_UINT16, H5S_ALL,
+ r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL,
H5S_ALL, H5P_DEFAULT, data);
if ( r < 0 ) {
fprintf(stderr, "Couldn't write data\n");
@@ -72,6 +72,8 @@ int hdf5_write(const char *filename, const uint16_t *data,
return 1;
}
+ H5Gclose(gh);
+ H5Dclose(dh);
H5Fclose(fh);
return 0;
diff --git a/src/hdf5-file.h b/src/hdf5-file.h
index a263368a..7f74efa3 100644
--- a/src/hdf5-file.h
+++ b/src/hdf5-file.h
@@ -18,7 +18,7 @@
#include <stdint.h>
-extern int hdf5_write(const char *filename, const uint16_t *data,
+extern int hdf5_write(const char *filename, const double *data,
int width, int height);
extern int hdf5_read(struct image *image, const char *filename);
diff --git a/src/image.c b/src/image.c
index 874d1f4a..04d0f35d 100644
--- a/src/image.c
+++ b/src/image.c
@@ -14,6 +14,7 @@
#include <stdlib.h>
#include <assert.h>
#include <math.h>
+#include <stdio.h>
#include "image.h"
#include "utils.h"
diff --git a/src/image.h b/src/image.h
index c596abfe..35bd5da9 100644
--- a/src/image.h
+++ b/src/image.h
@@ -46,10 +46,22 @@ struct imagefeature {
/* An opaque type representing a list of image features */
typedef struct _imagefeaturelist ImageFeatureList;
+
+/* A 3D vector in reciprocal space */
+struct threevec
+{
+ double u;
+ double v;
+ double w;
+};
+
+
/* Structure describing an image */
struct image {
uint16_t *data;
+ double *sfacs;
+ struct threevec *qvecs;
/* Radians. Defines where the pattern lies in reciprocal space */
double tilt;
diff --git a/src/main.c b/src/main.c
index 8ab11c85..ee02d340 100644
--- a/src/main.c
+++ b/src/main.c
@@ -19,7 +19,7 @@
#include <unistd.h>
#include "image.h"
-#include "relrod.h"
+#include "diffraction.h"
#include "cell.h"
#include "utils.h"
#include "hdf5-file.h"
@@ -39,11 +39,9 @@ static void main_show_help(const char *s)
int main(int argc, char *argv[])
{
- int c, i;
+ int c;
UnitCell *cell;
struct image image;
- int nrefl;
- float t;
while ((c = getopt(argc, argv, "h")) != -1) {
@@ -70,28 +68,21 @@ int main(int argc, char *argv[])
image.width = 512;
image.height = 512;
image.omega = deg2rad(40.0);
+ image.tilt = deg2rad(0.0);
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.lambda = 0.2e-9; /* LCLS wavelength */
- image.data = malloc(512*512*2);
-
- for ( t=0.0; t<180.0; t+=10.0 ) {
-
- char filename[32];
-
- memset(image.data, 0, 512*512*2);
- image.tilt = deg2rad(t);
-
- get_diffraction(&image, cell);
-
- /* Write the output file */
- snprintf(filename, 32, "simulated-%.0f.h5", t);
- hdf5_write(filename, image.data, image.width, image.height);
-
- }
+ image.qvecs = NULL;
+ image.sfacs = NULL;
+ image.data = NULL;
+
+ get_diffraction(&image, cell);
+
+ /* Write the output file */
+ hdf5_write("results/sim.h5", image.sfacs, image.width, image.height);
return 0;
}
diff --git a/src/utils.c b/src/utils.c
index 38855671..fa239f90 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -122,3 +122,37 @@ int sign(double a)
if ( a > 0 ) return +1;
return 0;
}
+
+
+void mapping_rotate(double x, double y, double z,
+ double *ddx, double *ddy, double *ddz,
+ double omega, double tilt)
+{
+ double nx, ny, nz;
+ double x_temp, y_temp, z_temp;
+
+ /* First: rotate image clockwise until tilt axis is aligned
+ * horizontally. */
+ nx = x*cos(omega) + y*sin(omega);
+ ny = -x*sin(omega) + y*cos(omega);
+ nz = z;
+
+ /* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the
+ * "wrong" way. This is because the crystal is rotated in the
+ * experiment, not the Ewald sphere. */
+ x_temp = nx; y_temp = ny; z_temp = nz;
+ nx = x_temp;
+ ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
+ nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
+
+ /* Finally, reverse the omega rotation to restore the location of the
+ * image in 3D space */
+ x_temp = nx; y_temp = ny; z_temp = nz;
+ nx = x_temp*cos(-omega) + y_temp*sin(-omega);
+ ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
+ nz = z_temp;
+
+ *ddx = nx;
+ *ddy = ny;
+ *ddz = nz;
+}
diff --git a/src/utils.h b/src/utils.h
index a60abe7e..f1b72af8 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -33,7 +33,10 @@ extern double distance3d(double x1, double y1, double z1,
extern size_t skipspace(const char *s);
extern void chomp(char *s);
extern int sign(double a);
-
+extern void mapping_rotate(double x, double y, double z,
+ double *ddx, double *ddy, double *ddz,
+ double omega, double tilt);
+
#define rad2deg(a) ((a)*180/M_PI)
#define deg2rad(a) ((a)*M_PI/180)
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 02:57:43 +0000