From fbef32528d3b4b2b2f694aa16d60c53511b1d402 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 9 Jul 2020 12:01:43 +0200 Subject: Add wavelength and electron_voltage, plus units, to geometry file --- libcrystfel/src/datatemplate.c | 111 +++++++++++++++++++++++++++++++++++- libcrystfel/src/datatemplate_priv.h | 13 +++++ libcrystfel/src/image.c | 42 ++++++++------ libcrystfel/src/utils.h | 16 +++++- 4 files changed, 162 insertions(+), 20 deletions(-) (limited to 'libcrystfel') diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c index e9d98639..e3b44207 100644 --- a/libcrystfel/src/datatemplate.c +++ b/libcrystfel/src/datatemplate.c @@ -628,6 +628,102 @@ static int parse_field_bad(struct dt_badregion *badr, const char *key, } +static int parse_electron_voltage(const char *val, + char **p_from, + enum wavelength_unit *punit) +{ + char *valcpy; + char *sp; + + valcpy = strdup(val); + if ( valcpy == NULL ) return 1; + + /* "electron_voltage" directive must have explicit units */ + sp = strchr(valcpy, ' '); + if ( sp == NULL ) { + free(valcpy); + return 1; + } + + if ( strcmp(sp+1, "V") == 0 ) { + *punit = WAVELENGTH_ELECTRON_V; + } else if ( strcmp(sp+1, "kV") == 0 ) { + *punit = WAVELENGTH_ELECTRON_KV; + } else { + free(valcpy); + return 1; + } + + sp[0] = '\0'; + *p_from = valcpy; + return 0; +} + + +static int parse_wavelength(const char *val, + char **p_from, + enum wavelength_unit *punit) +{ + char *valcpy; + char *sp; + + valcpy = strdup(val); + if ( valcpy == NULL ) return 1; + + /* "wavelength" directive must have explicit units */ + sp = strchr(valcpy, ' '); + if ( sp == NULL ) { + free(valcpy); + return 1; + } + + if ( strcmp(sp+1, "m") == 0 ) { + *punit = WAVELENGTH_M; + } else if ( strcmp(sp+1, "A") == 0 ) { + *punit = WAVELENGTH_A; + } else { + free(valcpy); + return 1; + } + + sp[0] = '\0'; + *p_from = valcpy; + return 0; +} + + +static int parse_photon_energy(const char *val, + char **p_from, + enum wavelength_unit *punit) +{ + char *valcpy; + char *sp; + + valcpy = strdup(val); + if ( valcpy == NULL ) return 1; + + /* "photon_energy" is the only one of the wavelength + * directives which is allowed to not have units */ + sp = strchr(valcpy, ' '); + if ( sp == NULL ) { + *punit = WAVELENGTH_PHOTON_EV; + } else if ( strcmp(sp+1, "eV") == 0 ) { + *punit = WAVELENGTH_PHOTON_EV; + sp[0] = '\0'; + } else if ( strcmp(sp+1, "keV") == 0 ) { + *punit = WAVELENGTH_PHOTON_KEV; + sp[0] = '\0'; + } else { + /* Unit specified, but unrecognised */ + free(valcpy); + return 1; + } + + *p_from = valcpy; + return 0; +} + + static int parse_toplevel(DataTemplate *dt, const char *key, const char *val, @@ -664,8 +760,19 @@ static int parse_toplevel(DataTemplate *dt, defaults->cnz_offset = atof(val); } else if ( strcmp(key, "photon_energy") == 0 ) { - /* Will be expanded when image is loaded */ - dt->wavelength_from = strdup(val); + return parse_photon_energy(val, + &dt->wavelength_from, + &dt->wavelength_unit); + + } else if ( strcmp(key, "electron_voltage") == 0 ) { + return parse_electron_voltage(val, + &dt->wavelength_from, + &dt->wavelength_unit); + + } else if ( strcmp(key, "wavelength") == 0 ) { + return parse_wavelength(val, + &dt->wavelength_from, + &dt->wavelength_unit); } else if ( strcmp(key, "peak_list") == 0 ) { dt->peak_list = strdup(val); diff --git a/libcrystfel/src/datatemplate_priv.h b/libcrystfel/src/datatemplate_priv.h index e94daf4f..75492db5 100644 --- a/libcrystfel/src/datatemplate_priv.h +++ b/libcrystfel/src/datatemplate_priv.h @@ -45,6 +45,17 @@ enum adu_per_unit }; +enum wavelength_unit +{ + WAVELENGTH_M, + WAVELENGTH_A, + WAVELENGTH_ELECTRON_KV, + WAVELENGTH_ELECTRON_V, + WAVELENGTH_PHOTON_KEV, + WAVELENGTH_PHOTON_EV +}; + + struct dt_rigid_group { char *name; @@ -179,6 +190,8 @@ struct _datatemplate char *wavelength_from; double photon_energy_bandwidth; /* Eww */ + enum wavelength_unit wavelength_unit; + unsigned int mask_bad; unsigned int mask_good; diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index cce29f73..47cbfea7 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -366,25 +366,31 @@ static double get_length(struct image *image, const char *from) } -static double get_wavelength(struct image *image, const char *from) +static double convert_to_m(double val, int units) { - char *units; - double value; + switch ( units ) { + + case WAVELENGTH_M : + return val; + + case WAVELENGTH_A : + return val * 1e-10; + + case WAVELENGTH_PHOTON_EV : + return ph_eV_to_lambda(val); + + case WAVELENGTH_PHOTON_KEV : + return ph_eV_to_lambda(val*1e3); + + case WAVELENGTH_ELECTRON_V : + return el_V_to_lambda(val); + + case WAVELENGTH_ELECTRON_KV : + return el_V_to_lambda(val*1e3); - units = get_value_and_units(image, from, &value); - if ( units == NULL ) { - /* Default unit is eV */ - return ph_eV_to_lambda(value); - } else { - if ( strcmp(units, "A") == 0 ) { - free(units); - return value * 1e-10; - } else { - ERROR("Invalid wavelength unit '%s'\n", units); - free(units); - return NAN; - } } + + return NAN; } @@ -491,7 +497,9 @@ struct image *image_read(DataTemplate *dtempl, const char *filename, if ( image == NULL ) return NULL; /* Wavelength might be needed to create detgeom (adu_per_eV) */ - image->lambda = get_wavelength(image, dtempl->wavelength_from); + image->lambda = convert_to_m(get_value(image, + dtempl->wavelength_from), + dtempl->wavelength_unit); create_detgeom(image, dtempl); diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index a43263e0..04208530 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -54,9 +54,12 @@ /* -------------------------- Fundamental constants ------------------------ */ -/* Electron charge in C */ +/* Electron charge (Coulombs) */ #define ELECTRON_CHARGE (1.6021773e-19) +/* Electron rest mass (kg) */ +#define ELECTRON_MASS (9.1093837015e-31) + /* Planck's constant (Js) */ #define PLANCK (6.62606896e-34) @@ -194,6 +197,17 @@ static inline int within_tolerance(double a, double b, double percent) /* Photon energy (eV) to k (1/m) */ #define ph_eV_to_k(a) ((a)*ELECTRON_CHARGE/PLANCK/C_VACUO) +/* Electron accelerating voltage (V) to wavelength (m) */ +static inline double el_V_to_lambda(double E) +{ + double Estar; + + /* Relativistically corrected accelerating voltage */ + Estar = E * (1.0 + E * ELECTRON_CHARGE/(2.0*ELECTRON_MASS*C_VACUO*C_VACUO)); + + return PLANCK / sqrt(2.0*ELECTRON_MASS*ELECTRON_CHARGE*Estar); +} + /* ------------------------------ Message logging ---------------------------- */ -- cgit v1.2.3