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

Loads of lattice stuff

parent 2008d998
......@@ -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@
......@@ -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;
}
}
}
/*
* 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;
}
}
}
/*
* 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 */
......@@ -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;
......
......@@ -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);
......
......@@ -14,6 +14,7 @@
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "image.h"
#include "utils.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;
......
......@@ -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;
}
......@@ -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;
}
......@@ -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)
......
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