Commit 8da20816 authored by Thomas White's avatar Thomas White
Browse files

render_hkl: Also do 2D

parent 8b54a98a
......@@ -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 );
......
......@@ -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;
}
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