1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
/*
* diffraction.c
*
* Calculate diffraction patterns by Fourier methods
*
* (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"
#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;
}
}
}
|
