Add wavelength and electron_voltage, plus units, to geometry file

4 files changed, 162 insertions, 20 deletions

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 ---------------------------- */