aboutsummaryrefslogtreecommitdiff
path: root/src/main.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/main.c')
-rw-r--r--src/main.c471
1 files changed, 471 insertions, 0 deletions
diff --git a/src/main.c b/src/main.c
new file mode 100644
index 0000000..a2fff85
--- /dev/null
+++ b/src/main.c
@@ -0,0 +1,471 @@
+/*
+ * main.c
+ *
+ * The Top Level Source File
+ *
+ * (c) 2006-2008 Thomas White <taw27@cam.ac.uk>
+ *
+ * synth2d - Two-Dimensional Crystallographic Fourier Synthesis
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <gtk/gtk.h>
+#include <fftw3.h>
+#include <string.h>
+#include <math.h>
+#include <gsl/gsl_errno.h>
+#include <assert.h>
+
+#include "displaywindow.h"
+#include "data.h"
+#include "main.h"
+#include "reflist.h"
+#include "normalise.h"
+#include "dpsynth.h"
+#include "argand.h"
+#include "symmetry.h"
+#include "gsf.h"
+#include "cflip.h"
+#include "cdm.h"
+#include "clean.h"
+#include "model.h"
+#include "geometry.h"
+#include "statistics.h"
+#include "elements.h"
+#include "refine.h"
+#include "superlattice.h"
+#include "options.h"
+
+
+ReflectionList *main_reflections = NULL;
+
+/* A rather shameless glue layer to avoid sharing information */
+void main_show_patterson() { displaywindow_show_patterson(main_reflections); }
+void main_show_pattersone() { displaywindow_show_pattersone(main_reflections); }
+void main_show_knownphases() { displaywindow_show_knownphases(main_reflections); }
+void main_show_calcphases() { displaywindow_show_calcphases(main_reflections); }
+void main_show_difference() {
+ ReflectionList *reflections = reflist_copy(main_reflections);
+ model_calculate_difference_coefficients(reflections);
+ displaywindow_show_difference(reflections);
+ reflist_free(reflections);
+}
+void main_show_refsyn() {
+ ReflectionList *reflections = reflist_copy(main_reflections);
+ model_calculate_refinement_coefficients(reflections);
+ displaywindow_show_refsyn(reflections);
+ reflist_free(reflections);
+}
+void main_show_diffpatt() {
+ ReflectionList *reflections = reflist_copy(main_reflections);
+ model_calculate_difference_coefficients(reflections);
+ displaywindow_show_diffpatt(reflections);
+ reflist_free(reflections);
+}
+void main_wilsonplot() { normalise_wilsonplot(main_reflections); }
+void main_falloffplot() { normalise_falloffplot(main_reflections); }
+void main_dpsynth() { dpsynth_main_open(main_reflections); }
+void main_dpsynth_update() { dpsynth_main_update(main_reflections); }
+void main_argand() { argand_open(main_reflections); }
+void main_argand_update() { argand_update(main_reflections); }
+void main_dethermalise(double level) { normalise_dethermalise(main_reflections, level); }
+void main_normalise_exponential(double a, double b, double c) { normalise_exponential(main_reflections, a, b, c); }
+void main_normalise(double level) { normalise_execute(main_reflections, level); }
+void main_gsf_initialise() { gsf_initialise(main_reflections); }
+void main_gsf_reset() { gsf_reset(main_reflections); }
+void main_aperture_open() { clean_aperture_open(main_reflections); }
+unsigned int main_cdm_tangentexpansion(CDMContext *cdm) {
+ cdm->reflections = main_reflections;
+ return cdm_tangentexpansion(cdm);
+}
+void main_display_phasing_solution(CDMContext *cdm, PhasingSolution *sol) {
+ cdm_display_phasing_solution(cdm, sol, main_reflections);
+}
+void main_symmetrise(Symmetry symmetry) {
+
+ char string[256];
+ double rsym;
+ SymFlags flags = 0;
+
+ if ( displaywindow_mode() == DWV_KNOWNPHASE ) {
+ flags = flags | SYMFLAG_PHASES_KNOWN;
+ } else if ( displaywindow_mode() == DWV_CALCPHASE ) {
+ flags = flags | SYMFLAG_PHASES_CALC;
+ } /* else amplitudes only */
+
+ rsym = symmetry_symmetrise(main_reflections, symmetry, flags);
+
+ main_displayr();
+ main_dpsynth_update();
+
+ displaywindow_switchview();
+ snprintf(string, 255, "Rsym=%.2f%%", rsym*100);
+ displaywindow_statusbar(string);
+
+}
+void main_geometry_correct(GeometryCorrectionType correction_type, double wavenumber, double circleradius) {
+ geometry_correct(main_reflections, correction_type, wavenumber, circleradius);
+}
+
+void main_superlattice_split(unsigned int xc, unsigned int yc) {
+
+ char string[256];
+
+ main_substitutereflections(superlattice_split(main_reflections, xc, yc));
+ displaywindow_switchview();
+ snprintf(string, 255, "Split %ix%i superlattice", xc, yc);
+ displaywindow_statusbar(string);
+
+}
+
+/* Cut out reflections with |g|<2.5 nm^-1 */
+void main_hpfilter() {
+
+ unsigned int i;
+ double a, b, c;
+ char string[256];
+ unsigned int n_del, n_orig;
+
+ a = data_a(); b = data_b(); c = data_c();
+ n_del = 0; n_orig = main_reflections->n_reflections;
+
+ for ( i=1; i<main_reflections->n_reflections; i++ ) {
+
+ double g;
+ signed int h = main_reflections->refs[i].h;
+ signed int k = main_reflections->refs[i].k;
+ signed int l = main_reflections->refs[i].l;
+
+ g = sqrt(((h*h)/(a*a)) + ((k*k)/(b*b)) + ((l*l)/(c*c)));
+
+ if ( g < 2.5 ) {
+ reflist_delref(main_reflections, h, k, l);
+ n_del++;
+ }
+ }
+
+ displaywindow_switchview();
+ snprintf(string, 255, "%0.2f%% of reflections filtered out (|g|<2.5nm^-1)", ((double)n_del/n_orig)*100);
+ displaywindow_statusbar(string);
+
+}
+
+void main_savereflections(const char *filename) {
+
+ FILE *fh;
+ unsigned int i;
+ int sim;
+
+ fh = fopen(filename, "w");
+
+ fprintf(fh, "a %f\n", data_a());
+ fprintf(fh, "b %f\n", data_b());
+ fprintf(fh, "c %f\n", data_c());
+ fprintf(fh, "angle %f\n", 180*data_gamma()/M_PI);
+ fprintf(fh, "scale %i\n", data_get_image_scale());
+
+ switch ( displaywindow_mode() ) {
+ case DWV_PATTERSON : sim = 0; break;
+ case DWV_PATTERSONE : sim = 0; break;
+ case DWV_KNOWNPHASE : sim = 0; break;
+ /* DWV_REALSPACE should never get here */
+ case DWV_DIFFERENCE : sim = 1; break;
+ case DWV_DIFFPATT : sim = 1; break;
+ case DWV_REFSYN : sim = 1; break;
+ case DWV_MODEL : sim = 1; break;
+ case DWV_SIMPATT : sim = 1; break;
+ case DWV_SIMFOLZPATT : sim = 1; break;
+ case DWV_EXITWAVE : sim = 1; break;
+ default : sim = 0; break;
+ }
+
+ if ( sim ) {
+
+ ReflectionList *calc;
+
+ switch ( displaywindow_mode() ) {
+ case DWV_DIFFERENCE : calc = reflist_copy(main_reflections);
+ model_calculate_difference_coefficients(calc); break;
+ case DWV_DIFFPATT : calc = reflist_copy(main_reflections);
+ model_calculate_difference_coefficients(calc); break;
+ case DWV_REFSYN : calc = reflist_copy(main_reflections);
+ model_calculate_refinement_coefficients(calc); break;
+ case DWV_MODEL : calc = model_calculate_f(NULL, NULL, 0); break;
+ case DWV_SIMPATT : calc = model_calculate_f(NULL, NULL, 0); break;
+ case DWV_SIMFOLZPATT : calc = model_calculate_f(NULL, NULL, 1); break;
+ case DWV_EXITWAVE : calc = model_calculate_f(NULL, NULL, 0); break;
+ default : calc = model_calculate_f(NULL, NULL, 0); break;
+ }
+
+ for ( i=1; i<calc->n_reflections; i++ ) {
+ signed int h = calc->refs[i].h;
+ signed int k = calc->refs[i].k;
+ signed int l = calc->refs[i].l;
+ double am = calc->refs[i].amplitude;
+ if ( (displaywindow_mode() != DWV_DIFFPATT) && (displaywindow_mode() != DWV_SIMPATT) && (displaywindow_mode() != DWV_SIMFOLZPATT) ) {
+ fprintf(fh, "%3i %3i %3i %8f %8f\n", h, k, l, am, calc->refs[i].phase_known);
+ } else {
+ fprintf(fh, "%3i %3i %3i %8f\n", h, k, l, am);
+ }
+ }
+ reflist_free(calc);
+
+ } else {
+
+ for ( i=1; i<main_reflections->n_reflections; i++ ) {
+
+ signed int h = main_reflections->refs[i].h;
+ signed int k = main_reflections->refs[i].k;
+ signed int l = main_reflections->refs[i].l;
+ double am = main_reflections->refs[i].amplitude;
+
+ switch ( displaywindow_mode() ) {
+
+ case DWV_PATTERSON : fprintf(fh, "%3i %3i %3i %8f\n", h, k, l, am); break;
+ case DWV_PATTERSONE : fprintf(fh, "%3i %3i %3i %8f\n", h, k, l, am); break;
+ case DWV_KNOWNPHASE : if ( main_reflections->refs[i].phase_known_set )
+ fprintf(fh, "%3i %3i %3i %8f %8f\n", h, k, l, am, main_reflections->refs[i].phase_known);
+ break;
+ case DWV_CALCPHASE : if ( main_reflections->refs[i].phase_calc_set )
+ fprintf(fh, "%3i %3i %3i %8f %8f\n", h, k, l, am, main_reflections->refs[i].phase_calc);
+ break;
+ default : break;
+
+ }
+
+ }
+
+ }
+
+ fclose(fh);
+
+}
+
+static void main_reset() {
+
+ unsigned int i;
+ unsigned int phases_known;
+
+ if ( main_reflections ) {
+ free(main_reflections);
+ }
+ main_reflections = reflist_new();
+ memcpy(main_reflections, data_getreflections(), sizeof(ReflectionList));
+ data_free();
+
+ /* Decide whether to start with known phases or Patterson */
+ phases_known = 0;
+ for ( i=0; i<main_reflections->n_reflections; i++ ) {
+ if ( main_reflections->refs[i].phase_known_set ) {
+ phases_known++;;
+ }
+ }
+ printf("MA: Have phase values for %i of %i reflections\n", phases_known, main_reflections->n_reflections);
+ if ( phases_known == 0 ) {
+ displaywindow_show_patterson(main_reflections);
+ printf("MA: No phase values known - displaying Patterson map instead.\n");
+ displaywindow_statusbar("No phase values known - displaying Patterson map instead");
+ } else {
+ displaywindow_show_knownphases(main_reflections);
+ displaywindow_statusbar("Displaying known phases from file");
+ }
+
+ displaywindow_brightness_auto(NULL, NULL);
+
+}
+
+int main(int argc, char *argv[]) {
+
+ /* Sort out configuration */
+ options_load();
+
+ g_thread_init(NULL);
+ gtk_init(&argc, &argv);
+
+ if ( argc != 2 ) {
+ fprintf(stderr, "Syntax: %s <data file>\n", argv[0]);
+ return 1;
+ }
+
+ gsl_set_error_handler_off();
+
+ /* Read data */
+ if ( data_read(argv[1]) ) {
+ return 1;
+ }
+
+ elements_initialise();
+ model_default();
+
+ /* Open main window */
+ displaywindow_open(argv[1]);
+
+ /* Execute first transform */
+ main_reset();
+
+ gtk_main();
+
+ options_save();
+
+ return 0;
+
+}
+
+void main_substitutereflections(ReflectionList *new) {
+ memcpy(main_reflections, new, sizeof(ReflectionList));
+ displaywindow_switchview();
+}
+
+void main_displayr() {
+
+ ReflectionList *model_reflections;
+ char r_string[256];
+
+ model_reflections = model_calculate_f(main_reflections, NULL, 69);
+ dpsynth_simdp_update(model_reflections);
+ if ( !model_current_is_blank() ) {
+ snprintf(r_string, 255, "R1=%.2f%%, R2=%.2f%%", 100*stat_r(main_reflections, model_reflections),
+ 100*stat_r2(main_reflections, model_reflections));
+ displaywindow_statusbar(r_string);
+ }
+ reflist_free(model_reflections);
+
+}
+
+/* Caution! */
+ReflectionList *main_reflist() { return main_reflections; }
+
+unsigned int main_max_h() {
+ unsigned int max, i;
+ max = 0;
+ for ( i=1; i<main_reflections->n_reflections; i++ ) {
+ if ( abs(main_reflections->refs[i].h) > max ) max = abs(main_reflections->refs[i].h);
+ }
+ return max;
+}
+
+unsigned int main_max_k() {
+ unsigned int max, i;
+ max = 0;
+ for ( i=1; i<main_reflections->n_reflections; i++ ) {
+ if ( abs(main_reflections->refs[i].k) > max ) max = abs(main_reflections->refs[i].k);
+ }
+ return max;
+}
+
+void main_stripzero() {
+
+ int i;
+
+ i = reflist_inlist(main_reflections, 0, 0, 0);
+ assert(i == 0);
+ main_reflections->refs[i].amplitude = 0.0;
+
+}
+
+/* Anti-alias the pattern by restricting resolution to a circle */
+void main_antialias() {
+
+ unsigned int i;
+ double max;
+ double ph, pk, mh, mk;
+
+ /* Find the resolution in each direction */
+ max = 0.0;
+ for ( i=1; i<main_reflections->n_reflections; i++ ) {
+ double res;
+ /* Find reflections on centre line */
+ if ( main_reflections->refs[i].k != 0 ) continue;
+ if ( main_reflections->refs[i].h < 0 ) continue;
+ res = resolution(main_reflections->refs[i].h,
+ main_reflections->refs[i].k,
+ main_reflections->refs[i].l,
+ data_a(), data_b(), data_c(), data_gamma());
+ /* Find highest resolution in this direction */
+ if ( res > max ) max = res;
+ }
+ ph = max;
+
+ max = 0.0;
+ for ( i=1; i<main_reflections->n_reflections; i++ ) {
+ double res;
+ /* Find reflections on centre line */
+ if ( main_reflections->refs[i].k != 0 ) continue;
+ if ( main_reflections->refs[i].h > 0 ) continue;
+ res = resolution(main_reflections->refs[i].h,
+ main_reflections->refs[i].k,
+ main_reflections->refs[i].l,
+ data_a(), data_b(), data_c(), data_gamma());
+ /* Find highest resolution in this direction */
+ if ( res > max ) max = res;
+ }
+ mh = max;
+
+ max = 0.0;
+ for ( i=1; i<main_reflections->n_reflections; i++ ) {
+ double res;
+ /* Find reflections on centre line */
+ if ( main_reflections->refs[i].h != 0 ) continue;
+ if ( main_reflections->refs[i].k < 0 ) continue;
+ res = resolution(main_reflections->refs[i].h,
+ main_reflections->refs[i].k,
+ main_reflections->refs[i].l,
+ data_a(), data_b(), data_c(), data_gamma());
+ /* Find highest resolution in this direction */
+ if ( res > max ) max = res;
+ }
+ pk = max;
+
+ max = 0.0;
+ for ( i=1; i<main_reflections->n_reflections; i++ ) {
+ double res;
+ /* Find reflections on centre line */
+ if ( main_reflections->refs[i].h != 0 ) continue;
+ if ( main_reflections->refs[i].k > 0 ) continue;
+ res = resolution(main_reflections->refs[i].h,
+ main_reflections->refs[i].k,
+ main_reflections->refs[i].l,
+ data_a(), data_b(), data_c(), data_gamma());
+ /* Find highest resolution in this direction */
+ if ( res > max ) max = res;
+ }
+ mk = max;
+
+ /* Find the smallest */
+ double hm, km, m;
+ if ( mh < ph ) hm = mh; else hm = ph;
+ if ( mk < pk ) km = mk; else km = pk;
+ if ( hm < km ) m = hm; else m=km;
+
+ printf("Data resolution is %f nm^-1\n", ph);
+
+ i = 1;
+ while ( i < main_reflections->n_reflections ) {
+
+ double res;
+ signed int h, k, l;
+
+ h = main_reflections->refs[i].h;
+ k = main_reflections->refs[i].k;
+ l = main_reflections->refs[i].l;
+ res = resolution(h, k, l, data_a(), data_b(), data_c(),
+ data_gamma());
+ //printf("Resolution of %3i %3i %3i is %f - ", h, k, l, res);
+ if ( res >= m ) {
+ reflist_delref(main_reflections, h, k, l);
+ // printf("eliminated\n");
+ } else {
+ // main_reflections->refs[i].amplitude = 1000000000;
+ // printf("maximised\n");
+ i++;
+ }
+
+ }
+
+}