aboutsummaryrefslogtreecommitdiff
path: root/src/cell.c
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2009-10-14 17:35:18 +0200
committerThomas White <taw@bitwiz.org.uk>2009-10-14 17:35:18 +0200
commit4e1e2ef4472e28a146e6c83d053f1fc4d2419f88 (patch)
treeb1702fc01bb3483bd566b0e9b38e84e376216462 /src/cell.c
parent58d51685fbaad1cdfce7fcbd15ab70456f9ef328 (diff)
Calculate reflections using reciprocal cell
Diffstat (limited to 'src/cell.c')
-rw-r--r--src/cell.c67
1 files changed, 59 insertions, 8 deletions
diff --git a/src/cell.c b/src/cell.c
index 6fa91cd3..81638016 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -16,6 +16,9 @@
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
#include "cell.h"
#include "utils.h"
@@ -74,9 +77,9 @@ static void cell_update_crystallographic(UnitCell *cell)
cell->gamma = angle_between(cell->ax, cell->ay, cell->az,
cell->bx, cell->by, cell->bz);
- printf("a=%f nm\n", cell->a);
- printf("b=%f nm\n", cell->b);
- printf("c=%f nm\n", cell->c);
+ printf("a=%f nm\n", cell->a/1e9);
+ printf("b=%f nm\n", cell->b/1e9);
+ printf("c=%f nm\n", cell->c/1e9);
printf("alpha = %f deg\n", rad2deg(cell->alpha));
printf(" beta = %f deg\n", rad2deg(cell->beta));
printf("gamma = %f deg\n", rad2deg(cell->gamma));
@@ -104,7 +107,7 @@ UnitCell *cell_new()
void cell_set_parameters(UnitCell *cell, double a, double b, double c,
- double alpha, double beta, double gamma)
+ double alpha, double beta, double gamma)
{
if ( !cell ) return;
@@ -135,7 +138,7 @@ void cell_set_cartesian(UnitCell *cell,
UnitCell *cell_new_from_parameters(double a, double b, double c,
- double alpha, double beta, double gamma)
+ double alpha, double beta, double gamma)
{
UnitCell *cell;
@@ -149,9 +152,9 @@ UnitCell *cell_new_from_parameters(double a, double b, double c,
void cell_get_cartesian(UnitCell *cell,
- double *ax, double *ay, double *az,
- double *bx, double *by, double *bz,
- double *cx, double *cy, double *cz)
+ double *ax, double *ay, double *az,
+ double *bx, double *by, double *bz,
+ double *cx, double *cy, double *cz)
{
if ( !cell ) return;
@@ -159,3 +162,51 @@ void cell_get_cartesian(UnitCell *cell,
*bx = cell->bx; *by = cell->by; *bz = cell->bz;
*cx = cell->cx; *cy = cell->cy; *cz = cell->cz;
}
+
+
+void cell_get_reciprocal(UnitCell *cell,
+ double *asx, double *asy, double *asz,
+ double *bsx, double *bsy, double *bsz,
+ double *csx, double *csy, double *csz)
+{
+ int s;
+ gsl_matrix *m;
+ gsl_matrix *inv;
+ gsl_permutation *perm;
+
+ m = gsl_matrix_alloc(3, 3);
+ gsl_matrix_set(m, 0, 0, cell->ax);
+ gsl_matrix_set(m, 0, 1, cell->bx);
+ gsl_matrix_set(m, 0, 2, cell->cx);
+ gsl_matrix_set(m, 1, 0, cell->ay);
+ gsl_matrix_set(m, 1, 1, cell->by);
+ gsl_matrix_set(m, 1, 2, cell->cy);
+ gsl_matrix_set(m, 2, 0, cell->az);
+ gsl_matrix_set(m, 2, 1, cell->bz);
+ gsl_matrix_set(m, 2, 2, cell->cz);
+
+ /* Invert */
+ perm = gsl_permutation_alloc(m->size1);
+ inv = gsl_matrix_alloc(m->size1, m->size2);
+ gsl_linalg_LU_decomp(m, perm, &s);
+ gsl_linalg_LU_invert(m, perm, inv);
+ gsl_permutation_free(perm);
+ gsl_matrix_free(m);
+
+ /* Transpose */
+ gsl_matrix_transpose(inv);
+
+ *asx = gsl_matrix_get(inv, 0, 0);
+ *bsx = gsl_matrix_get(inv, 0, 1);
+ *csx = gsl_matrix_get(inv, 0, 2);
+ *asy = gsl_matrix_get(inv, 1, 0);
+ *bsy = gsl_matrix_get(inv, 1, 1);
+ *csy = gsl_matrix_get(inv, 1, 2);
+ *asz = gsl_matrix_get(inv, 2, 0);
+ *bsz = gsl_matrix_get(inv, 2, 1);
+ *csz = gsl_matrix_get(inv, 2, 2);
+
+ printf("a* = %f %f %f\n", *asx/1e9, *asy/1e9, *asz/1e9);
+ printf("b* = %f %f %f\n", *bsx/1e9, *bsy/1e9, *bsz/1e9);
+ printf("c* = %f %f %f\n", *csx/1e9, *csy/1e9, *csz/1e9);
+}