aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2010-06-10 20:07:36 -0700
committerThomas White <taw@bitwiz.org.uk>2010-06-10 20:07:36 -0700
commit8da20816630b540a1a65be3f02159d44c7574a7c (patch)
tree84a8b3821b6fc306a7f505e76d8d08e13c08a434
parent8b54a98a5e644e1b81f3d68ef49e37112d3e3e2d (diff)
render_hkl: Also do 2D
-rw-r--r--src/reflections.c8
-rw-r--r--src/render_hkl.c242
2 files changed, 201 insertions, 49 deletions
diff --git a/src/reflections.c b/src/reflections.c
index 4b48c38f..90406bee 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -115,14 +115,18 @@ double *read_reflections(const char *filename, unsigned int *counts,
char line[1024];
signed int h, k, l;
float intensity, ph;
+ char phs[1024];
int r;
rval = fgets(line, 1023, fh);
- r = sscanf(line, "%i %i %i %f %f", &h, &k, &l, &intensity, &ph);
+ r = sscanf(line, "%i %i %i %f %s", &h, &k, &l, &intensity, phs);
if ( r != 5 ) continue;
set_intensity(ref, h, k, l, intensity);
- if ( phases != NULL ) set_phase(phases, h, k, l, ph);
+ if ( phases != NULL ) {
+ ph = atof(phs);
+ set_phase(phases, h, k, l, ph);
+ }
if ( counts != NULL ) set_count(counts, h, k, l, 1);
} while ( rval != NULL );
diff --git a/src/render_hkl.c b/src/render_hkl.c
index c4555242..e18ae90d 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -22,6 +22,8 @@
#include <getopt.h>
#include <sys/types.h>
#include <sys/wait.h>
+#include <cairo.h>
+#include <cairo-pdf.h>
#include "utils.h"
#include "reflections.h"
@@ -37,76 +39,152 @@ static void show_help(const char *s)
"Render intensity lists using POV-ray.\n"
"\n"
" -h, --help Display this help message.\n"
+" --povray Render a 3D animation using POV-ray.\n"
" -j <n> Run <n> instances of POV-ray in parallel.\n"
+" --zone-axis Render a 2D zone axis pattern.\n"
"\n");
}
-int main(int argc, char *argv[])
+static void render_za(UnitCell *cell, double *ref, unsigned int *c)
{
- int c;
- UnitCell *cell;
- char *infile;
- double *ref1;
- signed int h, k, l;
- unsigned int *c1;
- FILE *fh;
+ cairo_surface_t *surface;
+ cairo_t *dctx;
+ double max_u, max_v, max_res, max_intensity, scale;
+ double sep_u, sep_v, max_r;
+ double as, bs, theta;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
- pid_t pids[MAX_PROC];
- float max;
- int nproc = 1;
- int i;
+ signed int h, k;
+ float wh, ht;
- /* Long options */
- const struct option longopts[] = {
- {"help", 0, NULL, 'h'},
- {0, 0, NULL, 0}
- };
+ wh = 1024;
+ ht = 1024;
- /* Short options */
- while ((c = getopt_long(argc, argv, "hj:", longopts, NULL)) != -1) {
+ surface = cairo_pdf_surface_create("za.pdf", wh, ht);
- switch (c) {
- case 'h' : {
- show_help(argv[0]);
- return 0;
- }
+ if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) {
+ fprintf(stderr, "Couldn't create Cairo surface\n");
+ cairo_surface_destroy(surface);
+ return;
+ }
- case 'j' : {
- nproc = atoi(optarg);
- break;
- }
+ dctx = cairo_create(surface);
- case 0 : {
- break;
- }
+ /* Black background */
+ cairo_rectangle(dctx, 0.0, 0.0, wh, ht);
+ cairo_set_source_rgb(dctx, 0.0, 0.0, 0.0);
+ cairo_fill(dctx);
- default : {
- return 1;
- }
+ max_u = 0.0; max_v = 0.0; max_intensity = 0.0;
+ max_res = 0.0;
+
+ /* Work out reciprocal lattice spacings and angles for this cut */
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ theta = angle_between(asx, asy, asz, bsx, bsy, bsz);
+ as = modulus(asx, asy, asz) / 1e9;
+ bs = modulus(bsx, bsy, bsz) / 1e9;
+ STATUS("theta=%f\n", rad2deg(theta));
+
+ for ( h=-INDMAX; h<INDMAX; h++ ) {
+ for ( k=-INDMAX; k<INDMAX; k++ ) {
+
+ double u, v, intensity, res;
+ int ct;
+
+ ct = lookup_count(c, h, k, 0);
+ if ( ct < 1 ) continue;
+
+ intensity = lookup_intensity(ref, h, k, 0) / (float)ct;
+
+ res = resolution(cell, h, k, 0);
+ if ( res > max_res ) max_res = res;
+
+ if ( intensity != 0 ) {
+ u = (double)h*as*sin(theta);
+ v = (double)h*as*cos(theta) + k*bs;
+ if ( fabs(u) > fabs(max_u) ) max_u = fabs(u);
+ if ( fabs(v) > fabs(max_v) ) max_v = fabs(v);
+ if ( fabs(intensity) > fabs(max_intensity) )
+ max_intensity = fabs(intensity);
}
}
+ }
- if ( (nproc > MAX_PROC) || (nproc < 1) ) {
- ERROR("Number of processes is invalid.\n");
- return 1;
+ max_res /= 1e9;
+ max_u /= 0.5;
+ max_v /= 0.5;
+ printf("Maximum resolution is %f nm^-1\n", max_res);
+
+ if ( max_intensity <= 0.0 ) {
+ max_r = 4.0;
+ goto out;
}
- infile = argv[optind];
+ /* Choose whichever scaling factor gives the smallest value */
+ scale = ((double)wh-50.0) / (2*max_u);
+ if ( ((double)ht-50.0) / (2*max_v) < scale ) {
+ scale = ((double)ht-50.0) / (2*max_v);
+ }
- cell = load_cell_from_pdb("molecule.pdb");
- c1 = new_list_count();
- ref1 = read_reflections(infile, c1, NULL);
- if ( ref1 == NULL ) {
- ERROR("Couldn't open file '%s'\n", infile);
- return 1;
+ sep_u = as * scale * cos(theta);
+ sep_v = bs * scale;
+ max_r = ((sep_u < sep_v)?sep_u:sep_v);
+
+ for ( h=-INDMAX; h<INDMAX; h++ ) {
+ for ( k=-INDMAX; k<INDMAX; k++ ) {
+
+ double u, v, intensity, val;
+ int ct;
+
+ ct = lookup_count(c, h, k, 0);
+ if ( ct < 1 ) continue;
+
+ intensity = lookup_intensity(ref, h, k, 0) / (float)ct;
+ val = 10.0*intensity/max_intensity;
+
+ u = (double)h*as*sin(theta);
+ v = (double)h*as*cos(theta) + k*bs;
+
+ cairo_arc(dctx, ((double)wh/2)+u*scale*2,
+ ((double)ht/2)+v*scale*2, max_r, 0, 2*M_PI);
+
+ cairo_set_source_rgb(dctx, val, val, val);
+ cairo_fill(dctx);
+
+ }
}
+out:
+ /* Centre marker */
+ cairo_arc(dctx, (double)wh/2,
+ (double)ht/2, max_r, 0, 2*M_PI);
+ cairo_set_source_rgb(dctx, 1.0, 0.0, 0.0);
+ cairo_fill(dctx);
+
+ cairo_surface_finish(surface);
+ cairo_destroy(dctx);
+}
+
+
+static void povray_render_animation(UnitCell *cell, double *ref,
+ unsigned int *c, unsigned int nproc)
+{
+ FILE *fh;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ pid_t pids[MAX_PROC];
+ float max;
+ int i;
+ signed int h, k, l;
+
fh = fopen("render.pov", "w");
- fprintf(fh, "/* POV-Ray scene written by Synth3D */\n\n");
+ fprintf(fh, "/* POV-Ray scene written by CrystFEL */\n\n");
fprintf(fh, "#include \"colors.inc\"\n");
fprintf(fh, "#include \"textures.inc\"\n\n");
fprintf(fh, "global_settings {\n");
@@ -232,9 +310,9 @@ int main(int argc, char *argv[])
int s;
float val, p, r, g, b, trans;
- if ( !lookup_count(c1, h, k, l) ) continue;
+ if ( !lookup_count(c, h, k, l) ) continue;
- val = lookup_intensity(ref1, h, k, l);
+ val = lookup_intensity(ref, h, k, l);
val = max-val;
@@ -344,6 +422,76 @@ int main(int argc, char *argv[])
int r;
waitpid(pids[i], &r, 0);
}
+}
+
+
+int main(int argc, char *argv[])
+{
+ int c;
+ UnitCell *cell;
+ char *infile;
+ double *ref;
+ unsigned int *cts;
+ int config_povray = 0;
+ int config_zoneaxis = 0;
+ unsigned int nproc = 1;
+
+ /* Long options */
+ const struct option longopts[] = {
+ {"help", 0, NULL, 'h'},
+ {"povray", 0, &config_povray, 1},
+ {"zone-axis", 0, &config_zoneaxis, 1},
+ {0, 0, NULL, 0}
+ };
+
+ /* Short options */
+ while ((c = getopt_long(argc, argv, "hj:", longopts, NULL)) != -1) {
+
+ switch (c) {
+ case 'h' : {
+ show_help(argv[0]);
+ return 0;
+ }
+
+ case 'j' : {
+ nproc = atoi(optarg);
+ break;
+ }
+
+ case 0 : {
+ break;
+ }
+
+ default : {
+ return 1;
+ }
+ }
+
+ }
+
+ if ( (nproc > MAX_PROC) || (nproc < 1) ) {
+ ERROR("Number of processes is invalid.\n");
+ return 1;
+ }
+
+ infile = argv[optind];
+
+ cell = load_cell_from_pdb("molecule.pdb");
+ cts = new_list_count();
+ ref = read_reflections(infile, cts, NULL);
+ if ( ref == NULL ) {
+ ERROR("Couldn't open file '%s'\n", infile);
+ return 1;
+ }
+
+ if ( config_povray ) {
+ povray_render_animation(cell, ref, cts, nproc);
+ } else if ( config_zoneaxis ) {
+ render_za(cell, ref, cts);
+ } else {
+ ERROR("Try again with either --povray or --zone-axis.\n");
+ }
+
return 0;
}