aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2013-02-06 19:15:06 +0100
committerThomas White <taw@bitwiz.org.uk>2013-02-06 19:15:06 +0100
commit53382f60c7f8955c016725885c04f510f7fe96ee (patch)
tree6dd7fcade933eb3f8d69de1d78fc2be275de5098
parent761da16096472d157b2dd9c871e81ab17cdef12d (diff)
Stuff for partial_sim
-rw-r--r--src/partial_sim.c66
1 files changed, 40 insertions, 26 deletions
diff --git a/src/partial_sim.c b/src/partial_sim.c
index facf14ef..4e96fbbc 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -54,11 +54,12 @@
#define NBINS 50
-static void mess_up_cell(UnitCell *cell, double cnoise)
+static void mess_up_cell(Crystal *cr, double cnoise)
{
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
+ UnitCell *cell = crystal_get_cell(cr);
//STATUS("Real:\n");
//cell_print(cell);
@@ -82,17 +83,18 @@ static void mess_up_cell(UnitCell *cell, double cnoise)
/* For each reflection in "partial", fill in what the intensity would be
* according to "full" */
-static void calculate_partials(RefList *partial, double osf,
+static void calculate_partials(Crystal *cr,
RefList *full, const SymOpList *sym,
int random_intensities,
pthread_mutex_t *full_lock,
unsigned long int *n_ref, double *p_hist,
- double *p_max, double max_q, UnitCell *cell)
+ double *p_max, double max_q)
{
Reflection *refl;
RefListIterator *iter;
+ double res;
- for ( refl = first_refl(partial, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -137,16 +139,17 @@ static void calculate_partials(RefList *partial, double osf,
}
}
- Ip = osf * p * If;
+ Ip = crystal_get_osf(cr) * p * If;
- bin = NBINS*2.0*resolution(cell, h, k, l) / max_q;
+ res = resolution(crystal_get_cell(cr), h, k, l);
+ bin = NBINS*2.0*res/max_q;
if ( (bin < NBINS) && (bin>=0) ) {
p_hist[bin] += p;
n_ref[bin]++;
if ( p > p_max[bin] ) p_max[bin] = p;
} else {
STATUS("Reflection out of histogram range: %e %i %f\n",
- resolution(cell, h, k, l), bin, p);
+ res, bin, p);
}
Ip = gaussian_noise(Ip, 100.0);
@@ -206,13 +209,14 @@ struct queue_args
unsigned long int n_ref[NBINS];
double p_max[NBINS];
- FILE *stream;
+ Stream *stream;
};
struct worker_args
{
struct queue_args *qargs;
+ Crystal *crystal;
struct image image;
/* Histogram for this image */
@@ -243,21 +247,34 @@ static void *create_job(void *vqargs)
static void run_job(void *vwargs, int cookie)
{
- double osf;
struct quaternion orientation;
struct worker_args *wargs = vwargs;
struct queue_args *qargs = wargs->qargs;
int i;
+ Crystal *cr;
+ UnitCell *cell;
+ RefList *reflections;
+
+ cr = crystal_new();
+ if ( cr == NULL ) {
+ ERROR("Failed to create crystal.\n");
+ return;
+ }
+ wargs->crystal = cr;
+ crystal_set_image(cr, &wargs->image);
- osf = gaussian_noise(1.0, 0.3);
+ crystal_set_osf(cr, gaussian_noise(1.0, 0.3));
+ crystal_set_profile_radius(cr, wargs->image.beam->profile_radius);
/* Set up a random orientation */
orientation = random_quaternion();
- wargs->image.indexed_cell = cell_rotate(qargs->cell, orientation);
+ cell = cell_rotate(qargs->cell, orientation);
+ crystal_set_cell(cr, cell);
+ cell_free(cell);
snprintf(wargs->image.filename, 255, "dummy.h5");
- wargs->image.reflections = find_intersections(&wargs->image,
- wargs->image.indexed_cell);
+ reflections = find_intersections(&wargs->image, cr);
+ crystal_set_reflections(cr, reflections);
for ( i=0; i<NBINS; i++ ) {
wargs->n_ref[i] = 0;
@@ -265,14 +282,14 @@ static void run_job(void *vwargs, int cookie)
wargs->p_max[i] = 0.0;
}
- calculate_partials(wargs->image.reflections, osf, qargs->full,
+ calculate_partials(cr, qargs->full,
qargs->sym, qargs->random_intensities,
&qargs->full_lock,
wargs->n_ref, wargs->p_hist, wargs->p_max,
- qargs->max_q, wargs->image.indexed_cell);
+ qargs->max_q);
/* Give a slightly incorrect cell in the stream */
- mess_up_cell(wargs->image.indexed_cell, qargs->cnoise);
+ mess_up_cell(cr, qargs->cnoise);
}
@@ -282,7 +299,7 @@ static void finalise_job(void *vqargs, void *vwargs)
struct queue_args *qargs = vqargs;
int i;
- write_chunk(qargs->stream, &wargs->image, NULL, STREAM_INTEGRATED);
+ write_chunk(qargs->stream, &wargs->image, NULL, 0, 1);
for ( i=0; i<NBINS; i++ ) {
qargs->n_ref[i] += wargs->n_ref[i];
@@ -295,8 +312,7 @@ static void finalise_job(void *vqargs, void *vwargs)
qargs->n_done++;
progress_bar(qargs->n_done, qargs->n_to_do, "Simulating");
- reflist_free(wargs->image.reflections);
- cell_free(wargs->image.indexed_cell);
+ crystal_free(wargs->crystal);
free(wargs);
}
@@ -315,7 +331,7 @@ int main(int argc, char *argv[])
char *sym_str = NULL;
SymOpList *sym;
UnitCell *cell = NULL;
- FILE *ofh;
+ Stream *stream;
int n = 2;
int random_intensities = 0;
char *save_file = NULL;
@@ -496,13 +512,12 @@ int main(int argc, char *argv[])
ERROR("You must give a filename for the output.\n");
return 1;
}
- ofh = fopen(output_file, "w");
- if ( ofh == NULL ) {
+ stream = open_stream_for_write(output_file);
+ if ( stream == NULL ) {
ERROR("Couldn't open output file '%s'\n", output_file);
return 1;
}
free(output_file);
- write_stream_header(ofh, argc, argv);
image.det = det;
image.width = det->max_fs;
@@ -511,7 +526,6 @@ int main(int argc, char *argv[])
image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
image.div = beam->divergence;
image.bw = beam->bandwidth;
- image.profile_radius = beam->profile_radius;
image.filename = malloc(256);
image.copyme = NULL;
@@ -528,7 +542,7 @@ int main(int argc, char *argv[])
qargs.random_intensities = random_intensities;
qargs.cell = cell;
qargs.template_image = &image;
- qargs.stream = ofh;
+ qargs.stream = stream;
qargs.cnoise = cnoise;
qargs.max_q = largest_q(&image);
@@ -574,7 +588,7 @@ int main(int argc, char *argv[])
}
- fclose(ofh);
+ close_stream(stream);
cell_free(cell);
free_detector_geometry(det);
free(beam);