code changes for release 2.2.0

.gitignore

@@ -2,5 +2,6 @@
 ~*
 *.bak
 *.ipfT*
+doc/html/*
 doc/latex/*
 

README.md

@@ -14,8 +14,13 @@ PEARL Procedures should be installed according to the regular Igor Pro guideline
 - Find the `HDF5.XOP` (`HDF5-64.xop` for Igor 7 64-bit) extension in the `Igor Pro Folder` under `More Extensions/File Loaders` (`More Extensions (64-bit)/File Loaders`), create a shortcut, and move the shortcut to the `Igor Extensions` folder next to your `User Procedures` folder.
 - Find the `HDF5 Help.ihf` next to `HDF5.XOP`, create a shortcut, and move the shortcut to the `Igor Help Files` folder next to your `User Procedures` folder.
 
-PEARL Procedures are compatible with Igor Pro versions 6.37 and 8.04, 32-bit and 64-bit.
+PEARL Procedures are tested on Igor 8.04, 64-bit.
 Please make sure to use the latest release version.
+
+While most of the code remains compatible with Igor 6.37, it is not tested and not supported.
+Importing recent PShell data files may requires Igor 8 due to changes in the HDF5 library.
+Igor 7 contains some bugs which affect PEARL Procedures and should not be used.
+
 As long as no Igor 8 specific features are used (long object names), the produced experiment files remain compatible with Igor 6.
 
 
@@ -34,12 +39,16 @@ Matthias Muntwiler, <mailto:matthias.muntwiler@psi.ch>
 Copyright
 ---------
 
-Copyright 2009-2020 by [Paul Scherrer Institut](http://www.psi.ch)
+Copyright 2009-2021 by [Paul Scherrer Institut](http://www.psi.ch)
 
 
 Release Notes
 =============
 
+## rev-distro-2.2.0
+
+- Updates, bugfixes and performance improvements in angle scan processing.
+
 ## rev-distro-2.1.0
 
 - Check compatibility of major features with Igor 8.
@@ -53,8 +62,3 @@ Release Notes
 
 - The interface of data reduction functions has changed to make data reduction more efficient in multi-peak fits. The supplied reduction functions and dialogs have been refactored. If you want to use your own reduction functions written for pre-2.0, you have to adapt them to the new interface.
 
-## rev-distro-1.1.1 
-
-- If you have upgraded PEARL Procedures from pre-1.1.1 and Igor breaks in pearl-elog.ipf while opening an experiment, please delete the ELOG preferences file `pearl-elog/preferences.pxp`. (Check the Igor Help to find the package preferences folder on your system.)
-
-

pearl/fermi-edge-analysis.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #include "pearl-area-profiles"
 
@@ -297,24 +298,3 @@ function show_shift(data)
 	shift_x += ecenter
 end
 
-/// calculate the energy resolution of the analyser
-///
-/// @param	epass	pass energy in eV
-/// @param	slit		analyser entrance slit in mm
-///
-/// @return energy resolution (FWHM)
-function analyser_energy_resolution(epass, slit)
-	variable epass
-	variable slit
-
-	variable respow
-	if (epass < 4)
-		respow = 1110
-	elseif (epass < 8)
-		respow = 1400
-	else
-		respow = 1750
-	endif
-	
-	return epass * max(0.2, slit) / 0.2 / respow	
-end

pearl/pearl-anglescan-panel.ipf

@@ -1,9 +1,12 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma version = 1.8
 #pragma IgorVersion = 6.2
 #pragma ModuleName = PearlAnglescanPanel
 #include "pearl-anglescan-process"
 #include "pearl-pmsco-import"
+#include "pearl-scienta-preprocess"
+#include "pearl-area-display"
 
 // copyright (c) 2018-20 Paul Scherrer Institut
 //
@@ -430,9 +433,9 @@ static function delete_rows(rows, data, theta, tilt, phi)
 	extract /free idx, idx, idx >= 0
 	variable nx = dimsize(data, 0)
 	variable ny = numpnts(idx)
-	theta = theta[idx]
+	theta[0,ny-1] = theta[idx]
-	tilt = tilt[idx]
+	tilt[0,ny-1] = tilt[idx]
-	phi = phi[idx]
+	phi[0,ny-1] = phi[idx]
 	redimension /n=(ny) theta, tilt, phi
 	duplicate /free data, data_copy
 	redimension /n=(nx,ny) data

pearl/pearl-anglescan-process.ipf

@@ -1,12 +1,13 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
-#pragma version = 1.8
+#pragma version = 1.9
 #pragma IgorVersion = 6.2
 #pragma ModuleName = PearlAnglescanProcess
 #include "pearl-vector-operations"
 #include "pearl-polar-coordinates"
 #include <New Polar Graphs>
 
-// copyright (c) 2013-17 Paul Scherrer Institut
+// copyright (c) 2013-21 Paul Scherrer Institut
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 //   you may not use this file except in compliance with the License.
@@ -64,7 +65,7 @@
 ///
 /// @author matthias muntwiler, matthias.muntwiler@psi.ch
 ///
-/// @copyright 2013-17 Paul Scherrer Institut @n
+/// @copyright 2013-21 Paul Scherrer Institut @n
 /// Licensed under the Apache License, Version 2.0 (the "License"); @n
 ///   you may not use this file except in compliance with the License. @n
 ///   You may obtain a copy of the License at
@@ -80,6 +81,70 @@
 /// PearlAnglescanProcess is declared in @ref pearl-anglescan-process.ipf.
 
 
+/// append an angle scan strip to another one
+///
+/// concatenate two angle scan strips including matching attribute waves
+/// and replace the first strip with the resulting waves.
+/// this is useful if a scan was interrupted and continues in a second data file.
+///
+/// all accompanying 1D waves which have a matching length and exist in both source and destination folders
+/// are concatenated and stored in the destination.
+/// 'accompanying waves' are those in the same folder as the 2D strip wave and those in the :attr sub-folder.
+///
+/// @attention this function modifies all matching waves in the data and attr folders of strip1!
+///	            consider a backup before calling the function, or work on a copy of the original data
+///            (cf. Igor's DuplicateDataFolder operation)!
+///
+/// @param[in,out]	strip1		2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan.
+///									this is the first source wave and destination wave.
+///
+/// @param[in]		strip2		2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan.
+///									this is the second source wave.
+///
+/// @return (string) semicolon-separated list of modified wave names (without folder path).
+///
+function /s strip_append(strip1, strip2)
+	wave strip1
+	wave strip2
+
+	dfref df1 = GetWavesDataFolderDFR(strip1)
+	dfref df2 = GetWavesDataFolderDFR(strip2)
+	variable ny1 = dimsize(strip1, 1)
+	variable ny2 = dimsize(strip2, 1)
+
+	concatenate /np=1 {strip2}, strip1
+	string modified = AddListItem(NameOfWave(strip1), "")
+
+	variable idf = 0
+	do
+		variable iw = 0
+		do
+			wave /z w1 = WaveRefIndexedDFR(df1, iw)
+			if (!WaveExists(w1))
+				break
+			endif
+
+			wave /z w2 = df2:$(NameOfWave(w1))
+			if (WaveExists(w1) && WaveExists(w2))
+				if ((DimSize(w1, 0) == ny1) && (DimSize(w1, 1) == 0) && (DimSize(w2, 0) == ny2) && (DimSize(w2, 1) == 0))
+					concatenate /np=0 {w2}, w1
+					modified = AddListItem(NameOfWave(w1), modified, "", Inf)
+				endif
+			endif
+			iw += 1
+		while(1)
+		
+		df1 = df1:attr
+		df2 = df2:attr
+		if ((DataFolderRefStatus(df1) != 1) || (DataFolderRefStatus(df2) != 1))
+			break
+		endif
+		idf += 1
+	while(idf < 2)
+	
+	return modified
+end
+
 /// delete a contiguous range of frames from a strip.
 ///
 /// this can be used to remove a region of bad frames due to, e.g., measurement problems.
@@ -299,7 +364,9 @@ function normalize_strip_phi(strip, theta, phi, [theta_offset, theta_range, chec
 	endif
 
 	// average over analyser angles
-	wave dist = ad_profile_y(strip, -inf, inf, "")
+	duplicate /free strip, strip_copy
+	MatrixFilter NanZapMedian strip_copy
+	wave dist = ad_profile_y(strip_copy, -inf, inf, "")
 	
 	// smooth distribution function
 	duplicate /free dist, dist_smoo
@@ -385,7 +452,9 @@ function normalize_strip_theta(strip, theta, [theta_offset, smooth_method, smoot
 	endif
 
 	// average over analyser angles
-	wave dist = ad_profile_y(strip, -inf, inf, "")
+	duplicate /free strip, strip_copy
+	MatrixFilter NanZapMedian strip_copy
+	wave dist = ad_profile_y(strip_copy, -inf, inf, "")
 	
 	// smooth distribution function
 	duplicate /free dist, dist_smoo
@@ -484,7 +553,9 @@ function normalize_strip_thetaphi(strip, theta, phi, [theta_offset, smooth_metho
 	endif
 
 	// average over analyser angles
-	wave dist = ad_profile_y(strip, -inf, inf, "")
+	duplicate /free strip, strip_copy
+	MatrixFilter NanZapMedian strip_copy
+	wave dist = ad_profile_y(strip_copy, -inf, inf, "")
 	
 	// smooth distribution function
 	duplicate /free dist, dist_smoo
@@ -542,7 +613,9 @@ function normalize_strip_theta_scans(strip, theta, [theta_offset, smooth_method,
 	endif
 
 	// average over analyser angles
-	wave dist = ad_profile_y(strip, -inf, inf, "")
+	duplicate /free strip, strip_copy
+	MatrixFilter NanZapMedian strip_copy
+	wave dist = ad_profile_y(strip_copy, -inf, inf, "")
 	
 	// smooth distribution function
 	duplicate /free dist, dist_smoo
@@ -635,6 +708,7 @@ function normalize_strip_2d(strip, theta, [theta_offset, smooth_method, smooth_f
 	variable ny = dimsize(strip, 1)
 	
 	duplicate /free strip, dist, alpha_int, theta_int
+	MatrixFilter NanZapMedian dist
 	theta_int = theta[q] - theta_offset
 	alpha_int = dimoffset(strip, 0) + p * dimdelta(strip, 0)
 	redimension /n=(nx * ny) dist, alpha_int, theta_int
@@ -1955,7 +2029,7 @@ static function /s display_polar_graph(graphname, [angle_offset, do_ticks])
 		WMPolarGraphSetVar(graphname, "doAngleTickLabelSubRange", 1)
 		WMPolarGraphSetVar(graphname, "angleTickLabelRangeStart", 0)
 		WMPolarGraphSetVar(graphname, "angleTickLabelRangeExtent", 90)
-		WMPolarGraphSetStr(graphname, "angleTickLabelNotation", "%g°")
+		WMPolarGraphSetStr(graphname, "angleTickLabelNotation", "%g°")
   
 		WMPolarGraphSetVar(graphname, "doPolarGrids", 0)
 		WMPolarGraphSetVar(graphname, "doRadiusTickLabels", 0)
@@ -2083,9 +2157,9 @@ static function /s draw_hemi_axes(graphname, [do_grids])
 		SetDrawEnv /W=$graphname textxjust= 1,textyjust= 2
 		SetDrawEnv /W=$graphname save
 		radi = calc_graph_radius(30, projection=projection)
-		DrawText /W=$graphname radi, -0.1, "30°"
+		DrawText /W=$graphname radi, -0.1, "30°"
 		radi = calc_graph_radius(60, projection=projection)
-		DrawText /W=$graphname radi, -0.1, "60°"
+		DrawText /W=$graphname radi, -0.1, "60°"
 	endif
 	
 	setdatafolder savedf
@@ -2205,7 +2279,7 @@ function /s display_scanlines(nickname, alpha_lo, alpha_hi, m_theta, m_tilt, m_p
 	make /n=1 /d /free d_polar, d_azi
 	variable n_alpha = round(alpha_hi - alpha_lo) + 1
 	make /n=(n_alpha) /d /free analyser
-	setscale /i x alpha_lo, alpha_hi, "°", analyser
+	setscale /i x alpha_lo, alpha_hi, "°", analyser
 	analyser = x
 	
 	convert_angles_ttpa2polar(m_theta, loc_m_tilt, m_phi, analyser, d_polar, d_azi)
@@ -2526,50 +2600,78 @@ function set_polar_graph_cursor(nickname, cursorname, polar_angle, azim_angle, [
 	endif
 end
 
-/// add an arbitrary angle scan to a hemispherical scan grid.
+
+/// add arbitrary angle scan data to a hemispherical scan grid.
 ///
-/// the hemi grid must have been created in the current data folder by the make_hemi_grid function.
+/// the function fills the input data into bins defined by the hemi scan grid.
-/// the function determines the bin size at the given polar angle,
+/// it sums up the values and weights of the data points which fall into each bin,
-/// and adds all data points which fall into a bin.
+/// and adds the results to the totals and weights waves of the existing hemi grid.
-/// a point which lies exactly on the upper boundary falls into the next bin.
+/// finally, it updates the values wave (values divided by weights).
-/// this function does not clear previous values before adding new data.
+///
+/// the hemi grid must have been created in the current data folder by the make_hemi_grid() function.
+/// the function does not clear previous values before adding new data.
 /// values are added to the _tot wave, weights to the _wt wave.
-/// the intensity (_i) wave is calculated as _tot / _wt.
+/// the intensity (_i/values) wave is calculated as _tot divided by _wt.
+///
+///	 @param nickname			name prefix of holo waves.
+///								empty if waves are in current data folder.
+/// @param values			counts/intensity values at the positions given in the polara nd azi waves.
+///								one- or two-dimensional.
+///								NaN values are ignored.
+/// @param polar				polar angles (in degrees) of each data point.
+///								allowed range 0 <= theta <= 90.
+///								no specific ordering required.
+/// @param azi				azimuthal angles (in degrees) of each data point.
+///								allowed range -360 <= phi < +360.
+///								no specific order required.
+/// @param weights			weight or accumulation time of each point of values.
+///								weights must be positive. values with weight 0 are ignored.
+///								defaults to 1 if not specified.
+///
+/// the values, weights, polar and azi waves must have the same dimensions (one- or two-dimensional).
+/// no specific order is required, the function sorts (copies of) the arrays internally.
+///
+/// the actual binning is delegated to the thread-safe add_anglescan_worker() function
+/// under Igor's automatic multi-threading facility.
+///
 function hemi_add_anglescan(nickname, values, polar, azi, [weights])
-	string nickname // name prefix of holo waves.
+	string nickname
-		// may be empty.
+	wave values
-	wave values	// intensity values
+	wave polar
-		// the wave can be one- or two-dimensional.
+	wave azi
-		// no specific order required, the function sorts the arrays internally
+	wave weights
-	wave polar	// polar coordinates. allowed range 0 <= theta <= 90
+
-		// dimensions corresponding to value.
+	dfref savedf = GetDataFolderDFR()
-	wave azi	// azimuthal coordinates. allowed range -360 <= phi < +360
-		// dimensions corresponding to value.
-	wave weights	// total accumulation time of each point of values. default = 1
 
 	if (ParamIsDefault(weights))
 		duplicate /free values, weights
 		weights = 1
 	endif
-	
+
-	// quick check whether hemi grid is existing
 	string s_prefix = ""
 	string s_int = "values"
 	dfref df = find_hemi_data(nickname, s_prefix, s_int)
 
+	string s_totals = s_prefix + "tot"
+	string s_weights = s_prefix + "wt"
 	string s_polar = s_prefix + "pol"
 	string s_azim = s_prefix + "az"
+	string s_index = s_prefix + "index"
 	string s_theta = s_prefix + "th"
-
+	string s_dphi = s_prefix + "dphi"
-	wave /sdfr=df /z w_values = $s_int
+	string s_nphis = s_prefix + "nphis"
-	wave /sdfr=df /z w_azim = $s_azim
-	wave /sdfr=df /z w_polar = $s_polar
-	wave /sdfr=df /z w_theta = $s_theta
-	if (!waveexists(w_values) || !waveexists(w_azim) || !waveexists(w_polar))
-		abort "Missing hemispherical scan grid. Please call make_hemi_grid() first."
-	endif
 	
-	// make internal copies, one-dimensional, ordered in theta
+	wave /sdfr=df w_polar = $s_polar
+	wave /sdfr=df w_azim = $s_azim
+	wave /sdfr=df w_values = $s_int
+	wave /sdfr=df w_totals = $s_totals
+	wave /sdfr=df w_weights = $s_weights
+	wave /sdfr=df w_index = $s_index
+	wave /sdfr=df w_theta = $s_theta
+	wave /sdfr=df w_dphi = $s_dphi
+	wave /sdfr=df w_nphis = $s_nphis
+	
+	// make internal copies of input, one-dimensional, ordered in theta
 	duplicate /free values, values_copy
 	duplicate /free polar, polar_copy
 	duplicate /free azi, azi_copy
@@ -2577,49 +2679,186 @@ function hemi_add_anglescan(nickname, values, polar, azi, [weights])
 	variable nn = dimsize(values, 0) * max(dimsize(values, 1), 1)
 	redimension /n=(nn) values_copy, polar_copy, azi_copy, weights_copy
 	sort /r polar_copy, polar_copy, azi_copy, values_copy, weights_copy
+
+	make /n=(numpnts(w_theta)) /free /df dfw
+	// for debugging: remove the MultiThread keyword and the ThreadSafe keywords of sub-functions
+	MultiThread dfw = add_anglescan_worker(p, values_copy, weights_copy, polar_copy, azi_copy, w_polar, w_azim, w_theta, w_index, w_dphi, w_nphis)
+
+	variable pp
+	for (pp = 0; pp < numpnts(dfw); pp += 1)
+		dfref tdf= dfw[pp]
+		wave df_totals = tdf:w_totals
+		wave df_weights = tdf:w_weights
+		w_totals += df_totals
+		w_weights += df_weights
+	endfor
+	w_values = w_weights > 0 ? w_totals / w_weights : nan
 	
-	variable pol
+	SetDataFolder savedf
+end
+
+/// thread worker for hemi_add_anglescan
+///
+/// this function extracts one azimuthal scan from the input data and adds it to an existing holo scan.
+/// it should be considered as a part of hemi_add_anglescan and not used elsewhere,
+/// as its interface may change in the future.
+///
+/// the function takes as input the entire input data, an existing hemi grid and the index of a polar angle to work on.
+/// the results are a w_totals and w_weights wave that can be added to the hemi scan.
+/// the two waves are returned in a free data folder referenced by the return value of the function.
+/// the function does not change global data.
+///
+threadsafe static function /df add_anglescan_worker(ith, values, weights, polar, azi, w_polar, w_azim, w_theta, w_index, w_dphi, w_nphis)
+	variable ith		// index into w_theta
+	wave values		// input data: intensity/counts
+	wave weights		// input data: weights/dwell time
+	wave polar		// input data: polar angles
+	wave azi			// input data: azimuthal angles
+	wave w_polar		// hemi grid
+	wave w_azim		// hemi grid
+	wave w_theta		// hemi grid
+	wave w_index		// hemi grid
+	wave w_dphi		// hemi grid
+	wave w_nphis		// hemi grid
+	
+	dfref savedf= GetDataFolderDFR()
+	dfref freedf= NewFreeDataFolder()
+	SetDataFolder freedf
+	make /n=(numpnts(w_polar)) /d w_totals, w_weights
+	
+	variable pol = w_theta[ith]
 	variable pol_st = abs(w_theta[1] - w_theta[0])
-	variable pol1, pol2
+	variable pol1 = pol - pol_st / 2
+	variable pol2 = pol + pol_st / 2
+
+	extract /free /indx polar, sel, (pol1 < polar) && (polar <= pol2) && (numtype(values) == 0) && (weights > 0)
+	if (numpnts(sel) > 0)
+		duplicate /free /r=[0, numpnts(sel)-1] azi, azi_slice
+		duplicate /free /r=[0, numpnts(sel)-1] values, values_slice
+		duplicate /free /r=[0, numpnts(sel)-1] weights, weights_slice
+		azi_slice = azi[sel]
+		values_slice = values[sel]
+		weights_slice = weights[sel]
+		add_aziscan_core(values_slice, weights_slice, pol, azi_slice, w_theta, w_azim, w_index, w_dphi, w_totals, w_weights)
+	endif
+
+	SetDataFolder savedf
+	return freedf
+end
+
+/// thread worker for hemi_add_anglescan and hemi_add_aziscan
+///
+/// this function adds one azimuthal scan to an existing holo scan.
+/// it should be considered as a part of hemi_add_anglescan and hemi_add_aziscan and not used elsewhere,
+/// as its interface may change in the future.
+///
+/// the function takes as input an azimuthal scan and an existing hemi grid.
+/// the results are added to w_totals and w_weights waves.
+///
+/// @attention the function sorts the input arrays by azimuthal angle!
+///	            these waves must not refer to global objects if multi-threading is used!
+///
+threadsafe static function add_aziscan_core(values, weights, polar, azi, w_theta, w_azim, w_index, w_dphi, w_totals, w_weights)
+	wave values		// input data: intensity/counts
+	wave weights		// input data: weights/dwell time
+	variable polar	// input data: polar angle
+	wave azi			// input data: angle positions of the azimuthal scan
+						// acceptable range: >= -360 and < +360
+						// no specific order required, the function sorts the array in place (!)
+	wave w_theta		// hemi grid
+	wave w_azim		// hemi grid
+	wave w_index		// hemi grid
+	wave w_dphi		// hemi grid
+	wave w_totals	// output data: total counts in hemi grid order
+	wave w_weights	// output data: total weights in hemi grid order
 	
-	duplicate /free azi_copy, azi_slice
+	// destination slice coordinates
-	duplicate /free values_copy, values_slice
+	variable ipol = BinarySearch(w_theta, polar)
-	duplicate /free weights_copy, weights_slice
+	if (ipol < 0)
-	for (pol = 90; pol >= 0; pol -= pol_st)
+		return -1
-		pol1 = pol - pol_st / 2
+	endif
-		pol2 = pol + pol_st / 2
+	
-		extract /free /indx polar_copy, sel, (pol1 < polar_copy) && (polar_copy <= pol2)
+	variable d1, d2
-		if (numpnts(sel) > 0)
+	if (ipol >= 1)
-			redimension /n=(numpnts(sel)) azi_slice, values_slice, weights_slice
+		d1 = w_index[ipol - 1]
-			azi_slice = azi_copy[sel]
+	else
-			values_slice = values_copy[sel]
+		d1 = 0
-			weights_slice = weights_copy[sel]
+	endif
-			hemi_add_aziscan(nickname, values_slice, pol, azi_slice, weights=weights_slice)
+	d2 = w_index[ipol] - 1
+	variable nd = d2 - d1 + 1
+	variable dphi = w_dphi[ipol]
+	make /n=(nd+1) /free bin_index
+	setscale /i x w_azim[d1] - dphi/2, w_azim[d2] + dphi/2, "deg", bin_index
+	
+	// source slice coordinates	
+	// order the slice from -dphi/2 to 360-dphi/2
+	azi = azi < 0 ? azi + 360 : azi
+	azi = azi >= 360 - dphi/2 ? azi - 360 : azi
+	sort azi, values, weights, azi
+	setscale /p x 0, 1, "", values, weights, azi
+	
+	bin_index = BinarySearch(azi, x) + 1
+	bin_index = bin_index == -2 ? 0 : bin_index[p]
+	bin_index = bin_index == -1 ? numpnts(azi) : bin_index[p]
+	bin_index[nd] = numpnts(azi)
+	
+	// loop over destination
+	variable id
+	variable v1, v2, w1, w2
+	for (id = 0; id < nd; id += 1)
+		if (bin_index[id+1] > bin_index[id])
+			v1 = w_totals[d1 + id]
+			w1 = w_weights[d1 + id]
+			if ((numtype(v1) == 2) || (w1 <= 0))
+				v1 = 0
+				w1 = 0
+			endif
+			v2 = sum(values, bin_index[id], bin_index[id+1] - 1)
+			w2 = sum(weights, bin_index[id], bin_index[id+1] - 1)
+			w_totals[d1 + id] = v1 + v2
+			w_weights[d1 + id] = w1 + w2
 		endif
 	endfor
 end
 
-/// add an azimuthal scan to a hemispherical scan grid.
+/// add azimuthal data to a hemispherical scan grid.
 ///
-/// the hemi grid must have been created in the current data folder by the make_hemi_grid function.
+/// the hemi grid must have been created in the current data folder by the make_hemi_grid() function.
 /// the function determines the bin size at the given polar angle,
-/// and calculates the mean values of the data points which fall into a bin.
+/// sums up the values and weights of the data points which fall into each bin,
-/// a point which lies exactly on the upper boundary falls into the next bin.
+/// and adds the results to the totals and weights waves of the existing hemi grid.
+/// finally, it updates the values wave (values divided by weights).
+///
+///	 @param nickname			name prefix of holo waves.
+///								empty if waves are in current data folder.
+/// @param values			counts/intensity values of the azimuthal scan at the positions given in the azi parameter.
+/// @param polar				polar angle (in degrees) where to add the azi scan.
+/// @param azi				angle positions of the azimuthal scan.
+///								acceptable range: >= -360 and < +360.
+///								no specific order required, the function sorts the array internally.
+/// @param weights			weight or accumulation time of each point of values.
+///								defaults to 1 if not specified.
+///
+/// the actual binning is delegated to the thread-safe add_aziscan_core() function shared with hemi_add_anglescan().
+///
 function hemi_add_aziscan(nickname, values, polar, azi, [weights])
-	string nickname // name prefix of holo waves.
+	string nickname
-		// may be empty.
+	wave values
-	wave values // intensity values of the azimuthal scan at the positions given in the azi parameter
+	variable polar
-	variable polar // polar angle where to add the azi scan
+	wave azi
-	wave azi // angle positions of the azimuthal scan
+	wave weights
-		// acceptable range: >= -360 and < +360
-		// no specific order required, the function sorts the array internally
-	wave weights	// total accumulation time of each point of values. default = 1
 
+	dfref savedf = GetDataFolderDFR()
+
+	duplicate /free values, values_copy
+	duplicate /free azi, azi_copy
 	if (ParamIsDefault(weights))
-		duplicate /free values, weights
+		duplicate /free values, weights_copy
-		weights = 1
+		weights_copy = 1
+	else
+		duplicate /free weights, weights_copy
 	endif
-		
+
 	// hemi grid waves
 	string s_prefix = ""
 	string s_int = "values"
@@ -2644,56 +2883,11 @@ function hemi_add_aziscan(nickname, values, polar, azi, [weights])
 	wave /sdfr=df w_dphi = $s_dphi
 	wave /sdfr=df w_nphis = $s_nphis
 	
-	// destination slice coordinates
+	add_aziscan_core(values_copy, weights_copy, polar, azi_copy, w_theta, w_azim, w_index, w_dphi, w_totals, w_weights)
-	//polar = round(polar)
+
-	//variable ipol = 90 - polar
+	w_values = w_weights > 0 ? w_totals / w_weights : nan
-	variable ipol = BinarySearch(w_theta, polar)
+
-	if (ipol < 0)
+	SetDataFolder savedf
-		abort "assertion failed in hemi_add_aziscan(): polar angle not found in grid."
-	endif
-	
-	variable d1, d2
-	if (ipol >= 1)
-		d1 = w_index[ipol - 1]
-	else
-		d1 = 0
-	endif
-	d2 = w_index[ipol] - 1
-	variable nd = d2 - d1 + 1
-	variable dphi = w_dphi[ipol]
-	variable az1, az2
-	
-	// source slice coordinates	
-	// order the slice from -dphi/2 to 360-dphi/2
-	azi = azi < 0 ? azi + 360 : azi
-	azi = azi >= 360 - dphi/2 ? azi - 360 : azi
-	duplicate /free values, sel_values
-	duplicate /free weights, sel_weights
-	
-	// loop over destination
-	variable id
-	variable v1, v2, w1, w2
-	for (id = 0; id < nd; id += 1)
-		az1 = (id - 0.5) * dphi
-		az2 = (id + 0.5) * dphi
-		extract /free /indx azi, sel, (az1 <= azi) && (azi < az2)
-		if (numpnts(sel) > 0)
-			redimension /n=(numpnts(sel)) sel_values, sel_weights
-			sel_values = values[sel]
-			sel_weights = weights[sel]
-			v1 = w_totals[d1 + id]
-			w1 = w_weights[d1 + id]
-			if ((numtype(v1) == 2) || (w1 <= 0))
-				v1 = 0
-				w1 = 0
-			endif
-			v2 = sum(sel_values)
-			w2 = sum(sel_weights)
-			w_totals[d1 + id] = v1 + v2
-			w_weights[d1 + id] = w1 + w2
-		endif
-	endfor
-	w_values[d1, d1 + nd - 1] = w_totals[p] / w_weights[p]
 end
 
 /// interpolate a hemispherical scan onto a rectangular grid
@@ -3011,8 +3205,9 @@ function import_tpi_scan(nickname, theta, phi, intensity, [folding, npolar, nogr
 	variable ifold
 	duplicate /free phi, fold_phi
 	for (ifold = 0; ifold < folding; ifold += 1)
+		fold_phi = fold_phi >= 360 ? fold_phi - 360 : fold_phi
 		hemi_add_anglescan(nickname, intensity, theta, fold_phi)
-		fold_phi = fold_phi >= 180 ? fold_phi + 360 / folding - fold_phi : fold_phi + 360 / folding
+		fold_phi += 360 / folding
 	endfor
 
 	if (nograph==0)
@@ -3199,7 +3394,7 @@ function /wave hemi_azi_cut(nickname, pol)
 		make /n=(nsel) /o $s_cut
 		wave w_cut = $s_cut
 		w_cut = w_values[sel]
-		setscale /i x w_azim[sel[0]], w_azim[sel[nsel-1]], "°", w_cut
+		setscale /i x w_azim[sel[0]], w_azim[sel[nsel-1]], "°", w_cut
 		setdatafolder savedf
 		return w_cut
 	else
@@ -3209,23 +3404,37 @@ function /wave hemi_azi_cut(nickname, pol)
 	setdatafolder savedf
 end
 
-static function check_contrast(values, pcmin, pcmax, vmin, vmax)
+static function check_contrast(values, pcmin, pcmax, vmin, vmax, sym)
 	wave values
 	variable pcmin
 	variable pcmax
 	variable &vmin
 	variable &vmax
+	variable sym
 	
 	dfref save_df = GetDataFolderDFR()
 	dfref dfr = NewFreeDataFolder()
 	setdatafolder dfr
 	StatsQuantiles /inan /iw /q /z values
-	wave index = w_quantilesindex
+	wave /z index = w_quantilesindex
-	variable imin = round(numpnts(index) * pcmin / 100)
-	variable imax = round(numpnts(index) * (100 - pcmax) / 100)
-	vmin = values[index[imin]]
-	vmax = values[index[imax]]
 	setdatafolder save_df
+
+	if (waveexists(index))
+		variable imin = round(numpnts(index) * pcmin / 100)
+		variable imax = round(numpnts(index) * (100 - pcmax) / 100)
+		vmin = values[index[imin]]
+		vmax = values[index[imax]]
+		if (sym)
+			variable d = vmax - vmin
+			if ((vmax >= d/4) && (-vmin >= d/4))
+				vmax = min(abs(vmin), abs(vmax))
+				vmin = -vmax
+			endif
+		endif
+	else
+		vmin = wavemin(values)
+		vmax = wavemax(values)
+	endif
 end
 
 /// set the pseudocolor contrast by percentile.
@@ -3243,12 +3452,21 @@ end
 /// @param	pcmax		percentile above the maximum color (0-100).
 /// @param	graphname	name of graph. default: top graph.
 /// @param	colortable	name of new colortable. default: keep current table.
+/// @param	reversecolors	reverse colors of new colorable.
+///							takes effect only if colortable argument is defined.
+///							@arg 0 (default) normal colors,
+///							@arg 1 reverse color table
+/// @param	symmetric	make scale symmetric about zero (for modulation functions, e.g.).
+///							@arg 0 (default) do not enforce symmetry.
+///							@arg 1 try symmetric scale if "reasonable".
 ///
-function set_contrast(pcmin, pcmax, [graphname, colortable])
+function set_contrast(pcmin, pcmax, [graphname, colortable, reversecolors, symmetric])
 	variable pcmin
 	variable pcmax
 	string graphname
 	string colortable
+	variable reversecolors
+	variable symmetric
 	
 	if (ParamIsDefault(graphname))
 		graphname = ""
@@ -3256,6 +3474,12 @@ function set_contrast(pcmin, pcmax, [graphname, colortable])
 	if (ParamIsDefault(colortable))
 		colortable = ""
 	endif
+	if (ParamIsDefault(reversecolors))
+		reversecolors = 0
+	endif
+	if (ParamIsDefault(symmetric))
+		symmetric = 0
+	endif
 	
 	dfref save_df = GetDataFolderDFR()
 	
@@ -3285,9 +3509,12 @@ function set_contrast(pcmin, pcmax, [graphname, colortable])
 				rev = str2num("0" + StringFromList(4, info, ","))
 				if (strlen(colortable) > 0)
 					ctab = colortable
+					rev = reversecolors
+				endif
+				check_contrast(w, pcmin, pcmax, vmin, vmax, symmetric)
+				if (vmax > vmin)
+					ModifyGraph /w=$graphname zColor($objname)={w, vmin, vmax, $ctab, rev}
 				endif
-				check_contrast(w, pcmin, pcmax, vmin, vmax)
-				ModifyGraph /w=$graphname zColor($objname)={w, vmin, vmax, $ctab, rev}
 			endif
 		endif
 	endfor
@@ -3307,9 +3534,12 @@ function set_contrast(pcmin, pcmax, [graphname, colortable])
 				rev = str2num("0" + StringFromList(3, info, ","))
 				if (strlen(colortable) > 0)
 					ctab = colortable
+					rev = reversecolors
+				endif
+				check_contrast(w, pcmin, pcmax, vmin, vmax, symmetric)
+				if (vmax > vmin)
+					ModifyImage /w=$graphname $objname ctab={vmin, vmax, $ctab, rev}
 				endif
-				check_contrast(w, pcmin, pcmax, vmin, vmax)
-				ModifyImage /w=$graphname $objname ctab={vmin, vmax, $ctab, rev}
 			endif
 		endif
 	endfor

pearl/pearl-anglescan-tracker.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "Windows-1252"
 #pragma rtGlobals=3
 #pragma version = 1.4
 #pragma IgorVersion = 6.2

pearl/pearl-area-display.ipf

@@ -1,8 +1,10 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.2
 #pragma ModuleName = PearlAreaDisplay
 #pragma version = 1.04
 #include "pearl-compat"
+#include "pearl-area-profiles"
 
 /// @file
 /// @brief visualization tools for 2D and 3D data.

pearl/pearl-area-import.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "Windows-1252"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.2
 #pragma ModuleName = PearlAreaImport
@@ -5,7 +6,7 @@
 #include "pearl-compat"
 #include "pearl-gui-tools"
 
-// copyright (c) 2013-18 Paul Scherrer Institut
+// copyright (c) 2013-20 Paul Scherrer Institut
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 //   you may not use this file except in compliance with the License.
@@ -26,7 +27,7 @@
 ///
 /// @author matthias muntwiler, matthias.muntwiler@psi.ch
 ///
-/// @copyright 2013-18 Paul Scherrer Institut @n
+/// @copyright 2013-20 Paul Scherrer Institut @n
 /// Licensed under the Apache License, Version 2.0 (the "License"); @n
 ///   you may not use this file except in compliance with the License. @n
 ///   You may obtain a copy of the License at
@@ -1331,7 +1332,7 @@ function adh5_reduce_brick(source, reduction_func, reduction_param, result_prefi
 			for (iw = 0; iw < func_result; iw += 1)
 				sw = "redw_" + num2str(iw)
 				wave profile = dfr:$sw
-				sw = "ReducedData" + num2str(iw+1)
+				sw = result_prefix + num2str(iw+1)
 				make /n=(dimsize(profile, 0), nz, nt) /d /o $sw
 				wave data = $sw
 				setdimlabel 0, -1, $getdimlabel(profile, 0, -1), data

pearl/pearl-area-profiles.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.1
 #pragma ModuleName = PearlAreaProfiles

pearl/pearl-arpes.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.1
 #pragma ModuleName = PearlArpes
@@ -18,7 +19,7 @@
 #include "pearl-epics" // EPICS access under Igor
 #include "pearl-arpes-scans" // run ARPES scans under Igor
 #include "pearl-sample-tracker" // live tracking and adjustment of sample position
-#include "pearl-scienta-countrate"
+#include "pearl-scienta-live"
 #endif
 
 /// @file
@@ -84,4 +85,4 @@ end
 function UnloadPearlArpesPackage()
 	execute /p/q/z "DELETEINCLUDE \"pearl-arpes\""
 	execute /p/q/z "COMPILEPROCEDURES "
-end
+end

pearl/pearl-compat.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.1
 #pragma ModuleName = PearlCompat

pearl/pearl-data-explorer.ipf

@@ -1,7 +1,8 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
-#pragma IgorVersion = 6.1
+#pragma IgorVersion = 6.36
 #pragma ModuleName = PearlDataExplorer
-#pragma version = 1.50
+#pragma version = 1.60
 #include "pearl-area-import"
 #include "pearl-area-profiles"
 #include "pearl-area-display"
@@ -11,7 +12,7 @@
 #include "pearl-matrix-import"
 #endif
 
-// copyright (c) 2013-16 Paul Scherrer Institut
+// copyright (c) 2013-20 Paul Scherrer Institut
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 //   you may not use this file except in compliance with the License.
@@ -258,6 +259,9 @@ static function preview_file(filename)
 	string filename
 	
 	dfref saveDF = GetDataFolderDFR()
+	dfref previewDF = $package_path
+	
+	killStrings /z authors, pgroup, proposal, proposer, sample
 	
 	variable ft = pearl_file_type(filename)
 	switch(ft)
@@ -288,6 +292,35 @@ static function preview_file(filename)
 				sprintf cmd, "PearlElog#set_panel_graphs(\"\", \"%s\")", graphname
 				execute /Q/Z cmd
 			endif
+			svar /sdfr=previewDF /z authors
+			if (svar_Exists(authors))
+				if (strlen(authors)>=1)
+					sprintf cmd, "PearlElog#set_panel_attributes(\"\", \"author=%s\")", authors
+					execute /Q/Z cmd
+				endif
+			endif
+			svar /sdfr=previewDF /z pgroup
+			if (svar_Exists(pgroup))
+				if (strlen(pgroup)>=1)
+					sprintf cmd, "PearlElog#set_panel_attributes(\"\", \"p-group=%s\")", pgroup
+					execute /Q/Z cmd
+				endif
+			endif
+			svar /sdfr=previewDF /z proposal
+			if (svar_Exists(proposal))
+				if (strlen(proposal)>=1)
+					sprintf cmd, "PearlElog#set_panel_attributes(\"\", \"project=%s\")", proposal
+					execute /Q/Z cmd
+				endif
+			endif
+			svar /sdfr=previewDF /z proposer
+			svar /sdfr=previewDF /z sample
+			if (svar_Exists(sample))
+				if (strlen(sample)>=1)
+					sprintf cmd, "PearlElog#set_panel_attributes(\"\", \"sample=%s\")", sample
+					execute /Q/Z cmd
+				endif
+			endif
 		endif
 	endif
 	
@@ -337,6 +370,8 @@ static function /wave preview_pshell_file(filename)
 	
 	if (strlen(s_preview_file) > 0)
 		s_file_info = psh5_load_info("pearl_explorer_filepath", filename)
+		setdatafolder previewDF
+		psh5_load_general_group("pearl_explorer_filepath", filename)
 	else
 		s_file_info = ""
 	endif
@@ -960,15 +995,22 @@ static function /s display_preview_trace(xtrace, ytrace)
 		lab = "X"
 	endif
 	Label /w=$graphname bottom lab + " (\\U)"
-	lab = StringByKey("AxisLabelD", labels, "=", "\r")
+	lab = StringByKey("Dataset", labels, "=", "\r")
 	if (!strlen(lab))
 		lab = "value"
 	endif
 	Label /w=$graphname left lab + " (\\U)"
-	
+
 	return s_name
 end
 
+/// load the selected files
+///
+/// load the files that are selected in the data explorer panel.
+/// the files are loaded using the load_file() function.
+///
+/// @note this function may change the current data folder!
+///
 static function load_selected_files([options])
 	string options
 	
@@ -981,6 +1023,7 @@ static function load_selected_files([options])
 	variable ii
 	for (ii = 0; ii < nn; ii += 1)
 		if (wSelectedFiles[ii])
+			setdatafolder saveDF
 			if (ParamIsDefault(options))
 				load_file(wtFiles[ii])
 			else
@@ -990,15 +1033,19 @@ static function load_selected_files([options])
 	endfor
 	
 	update_datasets()
-	setdatafolder saveDF
 end
 
+/// load one file
+///
+/// this can be a PShell, HDF5, ITX or MTRX file.
+/// (HDF5 and MTRX files require the corresponding XOP to be loaded - cf. file documentation)
+///
+/// @note this function may change the current data folder!
+///
 static function load_file(filename, [options])
 	string filename
 	string options
 	
-	dfref saveDF = GetDataFolderDFR()
-
 	variable ft = pearl_file_type(filename)
 	switch(ft)
 		case 1:
@@ -1024,8 +1071,6 @@ static function load_file(filename, [options])
 		default:
 			break
 	endswitch
-
-	setdatafolder saveDF
 end
 
 static function prompt_hdf_options(options)
@@ -1087,6 +1132,27 @@ function prompt_func_params(func_name, func_param)
 	endif
 end
 
+/// load a pshell file
+///
+/// load a pshell hdf5 file (complete or reduced).
+///
+/// if options is not specified, the complete file is loaded.
+/// if options is an empty string, the package default options are used.
+///
+/// the only supported options is `mode:load_reduced`.
+/// in this case, the name of the reduction function must also be given under the `reduction_func` key.
+/// the reduction parameters are prompted for if a prompt function for the reduction function is found.
+/// the default reduction parameters are the most recent parameters `s_reduction_params` stored in the package data folder.
+///
+/// @param	options		`key:value;` list of load options.
+///								by default, load complete, using psh5_load_complete().
+///								empty string, use options from `s_hdf_options.
+///								@arg `mode:load_reduced` load reduced, using psh5_load_reduced().
+///								@arg `reduction_func:...` name of the reduction function.
+///
+/// @note after the function returns,
+///       the current data folder points to the loaded data (scan1).
+///
 static function /df load_pshell_file(filename, [options])
 	string filename
 	string options
@@ -1122,7 +1188,7 @@ static function /df load_pshell_file(filename, [options])
 					print reduction_func, reduction_params
 					psh5_load_reduced(nickname, "pearl_explorer_filepath", filename, $reduction_func, reduction_params)
 					svar s_filepath
-					loaded_filename = s_filepath
+					loaded_filename = filename
 				endif
 				break
 		endswitch
@@ -1132,10 +1198,12 @@ static function /df load_pshell_file(filename, [options])
 	if (strlen(loaded_filename) > 0)
 		setdatafolder $("root:" + nickname)
 		dataDF = GetDataFolderDFR()
+		setdatafolder $(":scan1")
 		string /g pearl_explorer_import = "load_pshell_file"
+	else
+		setdatafolder saveDF
 	endif
 	
-	setdatafolder saveDF
 	return dataDF
 end
 
@@ -1183,9 +1251,10 @@ static function /df load_hdf_file(filename, [options])
 		setdatafolder $("root:" + nickname)
 		dataDF = GetDataFolderDFR()
 		string /g pearl_explorer_import = "load_hdf_file"
+	else
+		setdatafolder saveDF
 	endif
 	
-	setdatafolder saveDF
 	return dataDF
 end
 

pearl/pearl-elog.ipf

@@ -1,10 +1,11 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
-#pragma version = 1.41
+#pragma version = 1.50
-#pragma IgorVersion = 6.2
+#pragma IgorVersion = 6.36
 #pragma ModuleName = PearlElog
 
 // author: matthias.muntwiler@psi.ch
-// Copyright (c) 2013-17 Paul Scherrer Institut
+// Copyright (c) 2013-20 Paul Scherrer Institut
 
 // Licensed under the Apache License, Version 2.0 (the "License");
 //   you may not use this file except in compliance with the License.
@@ -71,7 +72,7 @@
 ///
 /// @author matthias muntwiler, matthias.muntwiler@psi.ch
 ///
-/// @copyright 2013-17 Paul Scherrer Institut @n
+/// @copyright 2013-20 Paul Scherrer Institut @n
 /// Licensed under the Apache License, Version 2.0 (the "License"); @n
 ///   you may not use this file except in compliance with the License. @n
 ///   You may obtain a copy of the License at
@@ -272,14 +273,14 @@ function elog_init_pearl_templates()
 	
 	// attributes (persistent)
 	// available attributes
-	string /g attributes = "author;project;p-group;sample;source;task;technique;file;valid;"
+	string /g attributes = "author;project;pgroup;sample;source;task;technique;file;valid;"
 	// controls corresponding to attributes
 	// prefix determines the control type: sv_ = setvariable (string), pm_ = popup menu, cb = check box
 	string /g controls = "sv_author;sv_project;sv_pgroup;sv_sample;pm_source;pm_task;pm_technique;sv_file;cb_valid;"
 	// attributes with fixed options, value item declares the options string
 	string /g options = "source=sources;task=tasks;technique=techniques"
 	// attributes which must be defined
-	string /g required_attributes = "author;project;p-group;sample;source;task;technique;valid"
+	string /g required_attributes = "author;project;pgroup;sample;source;task;technique;valid"
 
 	// option lists
 	string /g sources = "Manual Entry;PShell;Scienta Data;SScan Data;Prosilica Data;OTF Data;Beamline Status;LEED Data;QMS Data;Matrix Data;Igor Pro;Other"
@@ -292,7 +293,7 @@ function elog_init_pearl_templates()
 	
 	// attributes (persistent)
 	// available attributes
-	string /g attributes = "author;project;p-group;sample;program;revision;machine;job;experiment;source path;result path;valid"
+	string /g attributes = "author;project;pgroup;sample;program;revision;machine;job;experiment;source path;result path;valid"
 	// controls corresponding to attributes
 	// prefix determines the control type: sv_ = setvariable (string), pm_ = popup menu, cb = check box
 	string /g controls = "sv_author;sv_project;sv_pgroup;sv_sample;pm_program;sv_revision;pm_machine;sv_job;sv_experiment;sv_sourcepath;sv_resultpath;cb_valid"

pearl/pearl-fitfuncs.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.2
 #pragma ModuleName = PearlFitFuncs
@@ -124,6 +125,51 @@ threadsafe function DoubletGaussLinBG_AO(pw, yw, xw) : FitFunc
 	yw += pw[2] * pw[3] * exp( -( (xw[p] - pw[4] + pw[5] /2) / pw[6] / pw[7] )^2 )
 end
 
+/// two doublet gaussian peaks on a linear background fit function (all at once).
+///
+/// this fits four gaussian peaks.
+/// peak positions are specified by center, splitting and shift rather than individually.
+/// amplitudes are specified as absolute values for peaks 1 and 3,
+/// and relative values for peaks 2 and 4.
+///
+/// can be used if a spin-orbit doublet is split by a chemical shift
+/// which affects both spin-orbit peaks equally.
+///
+/// @note FWHM = width * 2 * sqrt(ln(2)) = width * 1.665
+///
+/// @param	pw		shape parameters.
+///		@arg	pw[0] = constant coefficient of background
+///		@arg	pw[1] = linear coefficient of background
+///		@arg	pw[2] = amplitude of peak 1
+///		@arg	pw[3] = amplitude of peak 2
+///		@arg	pw[4] = amplitude of peak 3, relative to peak 1
+///		@arg	pw[5] = amplitude of peak 4, relative to peak 2
+///		@arg	pw[6] = position of peak 1
+///		@arg	pw[7] = splitting (distance between peaks 1 and 3)
+///		@arg	pw[8] = shift (distance between peaks 1 and 2)
+///		@arg	pw[9] = width of peaks 1 and 3
+///		@arg	pw[10] = width of peaks 2 and 4
+///
+/// @param	yw		y (dependent) values.
+///
+/// @param	xw		x (independent) independent values.
+///
+threadsafe function DblDoubletGaussLinBG_AO(pw, yw, xw) : FitFunc
+	wave pw
+	wave yw
+	wave xw
+
+	yw = pw[0] + xw[p] * pw[1]
+	// peak 1
+	yw += pw[2] * exp( -( (xw[p] - pw[6]) / pw[9] )^2 )
+	// peak 2
+	yw += pw[3] * exp( -( (xw[p] - pw[6] - pw[8]) / pw[10] )^2 )
+	// peak 3
+	yw += pw[2] * pw[4] * exp( -( (xw[p] - pw[6] - pw[7]) / pw[9] )^2 )
+	// peak 4
+	yw += pw[3] * pw[5] * exp( -( (xw[p] - pw[6] - pw[7] - pw[8]) / pw[10] )^2 )
+end
+
 //------------------------------------------------------------------------------
 // Voigt shapes
 //------------------------------------------------------------------------------
@@ -141,7 +187,7 @@ end
 ///	@arg	w[5 + (i-1) * 4] = shape of peak i
 /// @param	x		independent variable
 ///
-function MultiVoigtLinBG(w,x) : FitFunc
+threadsafe function MultiVoigtLinBG(w,x) : FitFunc
 	wave w
 	variable x
 
@@ -156,6 +202,39 @@ function MultiVoigtLinBG(w,x) : FitFunc
 	return v
 end
 
+/// multiple voigt peaks on a linear background fit function.
+///
+///
+/// this is the all-at-once version of @ref MultiVoigtLinBG.
+/// it runs slightly faster compared to the point-by-point function.
+///
+/// @param	pw		shape parameters.
+///					the length of the wave defines the number of peaks.
+///	@arg	pw[0] = constant coefficient of background
+///	@arg	pw[1] = linear coefficient of background
+///	@arg	pw[2 + (i-1) * 4] = amplitude of peak i
+///	@arg	pw[3 + (i-1) * 4] = position of peak i
+///	@arg	pw[4 + (i-1) * 4] = width of peak i
+///	@arg	pw[5 + (i-1) * 4] = shape of peak i
+///
+/// @param	yw		y (dependent) values.
+///
+/// @param	xw		x (independent) independent values.
+///
+threadsafe function MultiVoigtLinBG_AO(pw, yw, xw) : FitFunc
+	wave pw
+	wave yw
+	wave xw
+
+	variable np = numpnts(pw)
+	variable ip
+	
+	yw = pw[0] + xw[p] * pw[1]
+	for (ip = 2; ip < np; ip += 4)
+		yw += pw[ip] * VoigtFunc((xw[p] - pw[ip+1]) / pw[ip+2], pw[ip+3])
+	endfor
+end
+
 
 //------------------------------------------------------------------------------
 // Doniach-Sunjic shapes
@@ -198,7 +277,7 @@ end
 ///	@arg	w[5 + (i-1) * 4] = singularity index (0...1) of peak i
 /// @param	x		independent variable
 ///
-function MultiDoniachSunjicLinBG(w,x) : FitFunc
+threadsafe function MultiDoniachSunjicLinBG(w,x) : FitFunc
 	wave w
 	variable x
 

pearl/pearl-gui-tools.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.1
 #pragma ModuleName = PearlGuiTools

pearl/pearl-matrix-import.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3
 #pragma version = 1.00
 #pragma IgorVersion = 6.36

pearl/pearl-menu.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=1		// Use modern global access method.
 #pragma ModuleName = PearlMenu
 #pragma version = 1.02

pearl/pearl-pmsco-import.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.2
 #pragma ModuleName = PearlPmscoImport
@@ -174,6 +175,137 @@ static function /s save_scan_helper(destname, value, template, destdfr, wavename
 	return wavenames
 end
 
+/// load a PMSCO scan file into the current data folder.
+///
+/// the function loads all columns from the file.
+/// the waves are named with two-letter names according to the file extension.
+/// existing waves are overwritten.
+///
+/// the file extension must be `etpais` or a subset of it, e.g., `etpi`.
+/// the wave names will be `en` (energy), `th` (theta), `ph` (phi), `al` (alpha), 
+/// `in` (intensity) or `mo` (modulation), and `si` (sigma).
+///
+/// @param pathname		name of igor symbolic path to destination folder.
+///						prompt user if empty.
+///
+/// @param filename		requested file name.
+///						prompt user if empty.
+///						the extension must be a string of characters indicating the data of each column.
+///						it must be "etpais" or any substring of it, and the columns must be ordered accordingly.
+///						if the name contains `.modf`, the intensity wave is named as `mo` rather than `in`.
+///						this behaviour can be overridden by the `is_modulation` flag.
+///
+/// @param is_modulation	select whether the intensity column is named `mo` rather than `in`.
+///						@arg 0		(default) decide based on existens of `.modf` in the file name.
+///						@arg > 0	use `mo` regardless of file name.
+///						@arg < 0	use `in` regardless of file name.
+///
+/// @param	quiet (optional)
+///						@arg 0		(default)	print the file name and wave names to the history.
+///						@arg 1		do not print messages to the history.
+///
+function /s load_pmsco_scan(pathname, filename, [is_modulation, quiet])
+	string pathname	// name of a symbolic path
+	string filename
+	variable is_modulation
+	variable quiet
+	
+	if (ParamIsDefault(quiet))
+		quiet = 0
+	endif
+
+	loadwave /p=$pathname /a /g /o /q filename
+
+	if (ParamIsDefault(is_modulation))
+		is_modulation = StringMatch(s_filename, "*.modf.*")
+	else
+		is_modulation = is_modulation > 0
+	endif
+
+	string fileext
+	string waves = ""
+	if (v_flag > 0)
+		fileext = StringFromList(ItemsInList(s_filename, ".") - 1, s_filename, ".")
+		variable nw = ItemsInList(s_wavenames)
+		variable iw
+		string sw1, sw2
+		for (iw = 0; iw < nw; iw += 1)
+			sw1 = StringFromlist(iw, s_wavenames)
+			strswitch(fileext[iw])
+				case "e":
+					sw2 = "en"
+					break
+				case "t":
+					sw2 = "th"
+					break
+				case "p":
+					sw2 = "ph"
+					break
+				case "a":
+					sw2 = "al"
+					break
+				case "i":
+					if (is_modulation)
+						sw2 = "mo"
+					else
+						sw2 = "in"
+					endif
+					break
+				case "s":
+					sw2 = "si"
+					break
+			endswitch
+			duplicate /o $sw1, $sw2
+			killwaves /z $sw1
+			waves = AddListItem(sw2, waves, ",", inf)
+		endfor
+		
+		// Sort {en,th,ph, al} en,th,ph,al,int,sig
+		
+		if (!quiet)
+			print "load_pmsco_scan ",  s_filename, ": ", waves
+		endif
+		
+		return s_filename
+	else
+		return ""
+	endif
+end
+
+/// load a PMSCO result file into the current data folder.
+///
+/// result files have the extension dat or tasks.dat.
+/// this will overwrite existing waves.
+/// the function loads all columns.
+///
+/// @param	pathname		name of a symbolic path
+/// @param	filename			file name
+/// @param	quiet (optional)	@arg 0 (default)	print the file name and wave names to the history.
+///							@arg 1			do not print messages to the history.
+///
+function /s load_pmsco_result(pathname, filename, [quiet])
+	string pathname	// name of a symbolic path
+	string filename
+	variable quiet
+	
+	if (ParamIsDefault(quiet))
+		quiet = 0
+	endif
+	
+	if (quiet)
+		loadwave /p=$pathname /a /w /g /o /q filename
+	else
+		loadwave /p=$pathname /a /w /g /o filename
+	endif
+	
+	if (v_flag > 0)
+		return s_filename
+	else
+		return ""
+	endif
+end
+
+
 /// load an xyz cluster file
 ///
 /// load an xyz cluster file into the current data folder

pearl/pearl-polar-coordinates.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3
 #pragma version = 1.1
 #pragma IgorVersion = 6.1

pearl/pearl-pshell-import.ipf

@@ -1,12 +1,14 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.36
 #pragma ModuleName = PearlPShellImport
+#pragma version = 1.11
 #include <HDF5 Browser>
 #include "pearl-compat"
 #include "pearl-gui-tools"
 #include "pearl-area-import"
 
-// copyright (c) 2013-18 Paul Scherrer Institut
+// copyright (c) 2013-21 Paul Scherrer Institut
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 //   you may not use this file except in compliance with the License.
@@ -50,7 +52,7 @@
 ///
 /// @author matthias muntwiler, matthias.muntwiler@psi.ch
 ///
-/// @copyright 2013-18 Paul Scherrer Institut @n
+/// @copyright 2013-21 Paul Scherrer Institut @n
 /// Licensed under the Apache License, Version 2.0 (the "License"); @n
 ///   you may not use this file except in compliance with the License. @n
 ///   You may obtain a copy of the License at
@@ -75,7 +77,7 @@ strconstant kScanDimLabel = "scan"
 strconstant kDataDimLabel = "data"
 
 /// List of preferred datasets to load for preview
-strconstant kPreviewDatasets = "ScientaImage;ScientaSpectrum;ImageAngleDistribution;ImageEnergyDistribution;Counts;SampleCurrent;"
+strconstant kPreviewDatasets = "ImageEnergyDistribution;ScientaSpectrum;ScientaImage;Counts;SampleCurrent;"
 
 /// List of datasets that must be loaded to determine the axis scaling of a Scienta image
 strconstant kScientaScalingDatasets = "LensMode;ScientaChannelBegin;ScientaChannelEnd;ScientaSliceBegin;ScientaSliceEnd;"
@@ -336,6 +338,81 @@ function /s psh5_load_preview(APathName, AFileName, [load_data, load_attr, pref_
 	return dataname
 end
 
+/// load organizational metadata from the general group.
+///
+/// the general group contains the following datasets:
+/// authors, pgroup, proposal, proposer, sample.
+///
+/// data is loaded into the current data folder.
+/// all items are loaded into strings, authors is a comma-separated list.
+/// missing items default to empty strings.
+///
+/// @param APathName		igor symbolic path name. can be empty if the path is specified in FileName or a dialog box should be displayed
+///
+/// @param AFileName		if empty a dialog box shows up
+///
+/// @return					semicolon-separated list of the objects.
+///
+function /s psh5_load_general_group(APathName, AFileName)
+	string APathName
+	string AFileName
+
+	variable fileID
+
+	HDF5OpenFile /P=$APathName /R /Z fileID as AFileName
+	if (v_flag == 0)
+		string obj_names = "authors;pgroup;proposal;proposer;sample;"
+		variable nn = ItemsInList(obj_names, ";")
+		variable ii
+		string name
+	
+		for (ii = 0; ii < nn; ii += 1)
+			name = StringFromList(ii, obj_names, ";")
+			psh_load_general_string(fileID, name)
+		endfor
+	
+		return obj_names
+	else
+		return ""
+	endif
+end
+
+/// load a string from the general group.
+///
+/// the general group contains the following datasets:
+/// authors, pgroup, proposal, proposer, sample.
+///
+/// data is loaded into a global string in the current data folder.
+/// arrays with multiple items are loaded into a comma-separated list.
+/// a missing item defaults to the empty string.
+///
+/// @param	fileID			ID of open HDF5 file from psh5_open_file().
+///
+/// @return					comma-separated list of values.
+///
+function /s psh_load_general_string(fileID, name)
+	variable fileID
+	string name
+	
+	string path = "/general/" + name
+	HDF5LoadData /O /Q /Z /N=wt_load_general /TYPE=1 fileID, path
+	string values = ""
+	if (!v_flag)
+		wave /t wt_load_general
+		variable nn = numpnts(wt_load_general)
+		variable ii
+		for (ii = 0; ii < nn; ii += 1)
+			values = AddListItem(wt_load_general[ii], values, ",", inf)
+		endfor
+		killwaves /z wt_load_general
+		if (strlen(values) >= 1)
+			values = values[0,strlen(values)-2]
+		endif
+	endif
+	string /g $name = values
+	return values
+end
+
 /// load all data of a selected scan from a PShell data file.
 ///
 /// data is loaded into the current data folder.
@@ -1667,13 +1744,13 @@ function ps_detect_scale(ax, lo, hi, un)
 			case "Angular45":
 				lo[%$kAngleDimLabel] = -45/2
 				hi[%$kAngleDimLabel] = +45/2
-				un[%$kAngleDimLabel] = "°"
+				un[%$kAngleDimLabel] = "°"
 				ax[%$kAngleDimLabel] = "angle"
 				break
 			case "Angular60":
 				lo[%$kAngleDimLabel] = -60/2
 				hi[%$kAngleDimLabel] = +60/2
-				un[%$kAngleDimLabel] = "°"
+				un[%$kAngleDimLabel] = "°"
 				ax[%$kAngleDimLabel] = "angle"
 				break
 			case "Transmission":
@@ -1719,12 +1796,12 @@ function ps_detect_scale(ax, lo, hi, un)
 					un[%$kScanDimLabel] = "mm"
 					break
 				case "ExitSlit":
-					un[%$kScanDimLabel] = "µm"
+					un[%$kScanDimLabel] = "µm"
 					break
 				case "ManipulatorTheta":
 				case "ManipulatorTilt":
 				case "ManipulatorPhi":
-					un[%$kScanDimLabel] = "°"
+					un[%$kScanDimLabel] = "°"
 					break
 				case "FocusXRot":
 				case "FocusYRot":
@@ -1890,6 +1967,11 @@ end
 ///
 /// @param reduction_param	parameter string for the reduction function.
 ///
+/// @param dataset		name of dataset to load, optionally including group path relative to scan (scan 1).
+///							by default, the function looks for a ScientaImage dataset.
+///							in a multi-region scan, this will be region 1.
+///							to select region 2, e.g., use `dataset="region2/ScientaImage"`.
+///
 /// @param progress			progress window.
 ///	@arg	1				(default) show progress window
 ///	@arg	0				do not show progress window
@@ -1908,12 +1990,13 @@ end
 ///
 /// @return					global string s_scanpaths in new data folder contains a list of scan groups inside the file.
 ///
-function /s psh5_load_reduced(ANickName, APathName, AFileName, reduction_func, reduction_param, [progress, nthreads])
+function /s psh5_load_reduced(ANickName, APathName, AFileName, reduction_func, reduction_param, [dataset, progress, nthreads])
 	string ANickName
 	string APathName
 	string AFileName
 	funcref adh5_default_reduction reduction_func
 	string reduction_param
+	string dataset
 	variable progress
 	variable nthreads
 		
@@ -1972,7 +2055,9 @@ function /s psh5_load_reduced(ANickName, APathName, AFileName, reduction_func, r
 		
 		setdatafolder dataDF
 		string datasets = psh5_list_scan_datasets(fileID, scanpath, include_regions=1)
-		string dataset = select_dataset(datasets, "ScientaImage")
+		if (ParamIsDefault(dataset))
+			dataset = select_dataset(datasets, "ScientaImage")
+		endif
 		wavenames = psh5_load_dataset_reduced(fileID, scanpath, dataset, reduction_func, reduction_param, progress=progress, nthreads=nthreads)
 
 		psh5_close_file(fileID)
@@ -2503,3 +2588,48 @@ static function /s wave2list(w, format, sep)
 	
 	return list
 end
+
+/// kill any waves matching a pattern in the experiment
+///
+/// this may be used to kill big waves of original data before saving
+///
+function /s kill_matching_waves(dfr, pattern, recurse, [killed])
+	DFREF dfr
+	string pattern
+	variable recurse
+	string killed
+	
+	if (ParamIsDefault(killed))
+		killed = ""
+	endif
+	
+	string s
+	string r
+	variable index = 0
+	do
+		Wave/Z w = WaveRefIndexedDFR(dfr, index)
+		if (!WaveExists(w))
+			break
+		endif
+
+		s = NameOfWave(w)
+		if (stringmatch(s, pattern))
+			killwaves /z w
+			killed = AddListItem(s, killed, ";", Inf)
+		endif
+
+		index += 1
+	while(1)
+	
+	if (recurse)
+		Variable numChildDataFolders = CountObjectsDFR(dfr, 4)
+		Variable i
+		for(i=0; i<numChildDataFolders; i+=1)
+			String childDFName = GetIndexedObjNameDFR(dfr, 4, i)
+			DFREF childDFR = dfr:$childDFName
+			killed = kill_matching_waves(childDFR, pattern, 1, killed=killed)
+		endfor
+	endif
+	
+	return killed
+End

pearl/pearl-scienta-live.ipf

@@ -1,6 +1,7 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.35
-#pragma ModuleName = PearlScientaCountrate
+#pragma ModuleName = PearlScientaLive
 #include "pearl-area-display"
 
 // Copyright (c) 2019 Paul Scherrer Institut
@@ -11,10 +12,13 @@
 //   http://www.apache.org/licenses/LICENSE-2.0
 
 /// @file
-/// @brief count rate functions for Scienta detector images.
+/// @brief utility functions for operating the Scienta analyser.
 /// @ingroup ArpesPackage
 ///
-/// this procedure contains functions for working with true count rates.
+/// this procedure contains various utility functions for working with the Scienta analyser, e.g.
+///
+/// - statistical analysis of true count rates
+/// - estimate energy and angle resolution
 ///
 /// @author matthias muntwiler, matthias.muntwiler@psi.ch
 ///
@@ -24,10 +28,10 @@
 ///   You may obtain a copy of the License at
 ///   http://www.apache.org/licenses/LICENSE-2.0
 
-/// @namespace PearlScientaCountrate
+/// @namespace PearlScientaLive
-/// @brief count rate functions for Scienta detector images.
+/// @brief utility functions for operating the Scienta analyser.
 ///
-/// PearlScientaCountrate is declared in @ref pearl-scienta-countrate.ipf.
+/// PearlScientaLive is declared in @ref pearl-scienta-live.ipf.
 
 /// open live display of most recent scienta measurement
 ///
@@ -211,3 +215,25 @@ function check_exposure_opt(image, outmask, dwelltime, [calc_df])
 	
 	SetDataFolder save_df
 end
+
+/// calculate the energy resolution of the analyser
+///
+/// @param	epass	pass energy in eV
+/// @param	slit		analyser entrance slit in mm
+///
+/// @return energy resolution (FWHM)
+function analyser_energy_resolution(epass, slit)
+	variable epass
+	variable slit
+
+	variable respow
+	if (epass < 4)
+		respow = 1110
+	elseif (epass < 8)
+		respow = 1400
+	else
+		respow = 1750
+	endif
+	
+	return epass * max(0.2, slit) / 0.2 / respow
+end

pearl/pearl-scienta-preprocess.ipf

@@ -1,7 +1,9 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma IgorVersion = 6.1
 #pragma ModuleName = PearlScientaPreprocess
 #include "pearl-fitfuncs"
+#include "pearl-area-import"
 
 // Copyright (c) 2013-18 Paul Scherrer Institut
 //
@@ -804,6 +806,7 @@ threadsafe function /wave gauss4_reduction(source, param)
 		p1 -= ceil((ybox - 1) / 2)
 	endif
 	variable V_FitNumIters
+	variable V_FitError
 	
 	for (pp = p0; pp <= p1; pp += 1)
 		// box average
@@ -820,10 +823,11 @@ threadsafe function /wave gauss4_reduction(source, param)
 		wmax = wavemax(xprof)
 		w_coef[0] = wmin
 		w_coef[1] = 0
-
 		for (ipk=0; ipk < npeaks; ipk += 1)
 			w_coef[2 + ipk*3] = wmax - wmin
 		endfor
+		
+		V_FitError = 0
 		FuncFit /H=hold /Q /NTHR=1 /N /W=2 MultiGaussLinBG_AO w_coef xprof[pl,ph] /C={cmat, cvec} /I=1 /W=xprof_sig[pl,ph]
 		wave w_sigma
 	
@@ -948,6 +952,237 @@ function /s find_gauss4_reduction_params(spectrum, peakpos)
 	return param
 end
 
+/// fit horizontal cuts of an image with up to four voigtian peaks on a linear background
+///
+/// the function fits each horizontal profile (EDC) with four voigtian peaks on a linear background.
+/// the position, width and shape of the peaks is kept fixed according to input parameters.
+/// the peak amplitude is constrained to positive value.
+///
+/// the width and shape parameters are defined as in Igor's VoigtFunc:
+/// width is the width of the gaussian component,
+/// shape * width is the width of the lorentzian component.
+///
+/// @param source			source wave.
+///							two-dimensional distribution of counts.
+///							for correct weighting and error estimation it is important 
+///							that the source wave contains actual counts (Poisson statistics).
+///
+/// @param param			(in, out) semicolon-separated key=value list of processing parameters.
+///							this is a pass-by-reference argument.
+///							the following parameters are required.
+///							position, width and limit parameters are on the x (energy) scale.
+///			@arg rngl		low limit of fit interval
+///			@arg rngh		high limit of fit interval
+///			@arg pos1		position of peak 1
+///			@arg wid1		width of peak 1
+///			@arg shp1		shape of peak 1
+///			@arg pos2		position of peak 2
+///			@arg wid2		width of peak 2
+///			@arg shp2		shape of peak 2
+///			@arg pos3		position of peak 3
+///			@arg wid3		width of peak 3
+///			@arg shp3		shape of peak 3
+///			@arg pos4		position of peak 4
+///			@arg wid4		width of peak 4
+///			@arg shp4		shape of peak 4
+///			@arg npeaks		number of peaks to fit: 1...4
+///							the others are held at amplitude 0.
+///
+/// @return	free wave containing references of the result waves.
+///			the number of waves is two times the number of peaks that are fit.
+///			the first npeaks waves contain the peak amplitudes,
+///			the second npeaks waves the corresponding error estimates.
+///
+threadsafe function /wave voigt4_reduction(source, param)
+	wave source
+	string &param
+	
+	dfref orig_dfr = GetDataFolderDFR()
+	
+	variable nx = dimsize(source, 0)
+	variable ny = dimsize(source, 1)
+	
+	// read parameters
+	variable rngl = NumberByKey("rngl", param, "=", ";")
+	variable rngh = NumberByKey("rngh", param, "=", ";")
+	variable pos1 = NumberByKey("pos1", param, "=", ";")
+	variable wid1 = NumberByKey("wid1", param, "=", ";")
+	variable shp1 = NumberByKey("shp1", param, "=", ";")
+	variable pos2 = NumberByKey("pos2", param, "=", ";")
+	variable wid2 = NumberByKey("wid2", param, "=", ";")
+	variable shp2 = NumberByKey("shp2", param, "=", ";")
+	variable pos3 = NumberByKey("pos3", param, "=", ";")
+	variable wid3 = NumberByKey("wid3", param, "=", ";")
+	variable shp3 = NumberByKey("shp3", param, "=", ";")
+	variable pos4 = NumberByKey("pos4", param, "=", ";")
+	variable wid4 = NumberByKey("wid4", param, "=", ";")
+	variable shp4 = NumberByKey("shp4", param, "=", ";")
+	variable npeaks = NumberByKey("npeaks", param, "=", ";")
+	
+	// prepare curve fit
+	variable ipk
+	make /free xprof
+	adh5_setup_profile(source, xprof, 0)
+	duplicate /free xprof, xprof_sig
+	variable pl = max(x2pnt(xprof, rngl), 0)
+	variable ph = min(x2pnt(xprof, rngh), numpnts(xprof) - 1)
+	
+	make /free /n=(npeaks) peak_coef
+	variable coef_per_peak = 4
+	peak_coef = p * coef_per_peak + 2
+	variable n_coef = npeaks * coef_per_peak + 2
+	make /free /d /n=18 w_coef, W_sigma
+	w_coef[0] = {0, 0, 1, pos1, wid1, shp1, 1, pos2, wid2, shp2, 1, pos3, wid3, shp3, 1, pos4, wid4, shp4}
+	redimension /n=(n_coef)  w_coef, w_sigma
+
+	// text constraints cannot be used in threadsafe functions.
+	// the following matrix-vector formulation enforces all peak amplitudes to be positive.
+	make /free /n=(npeaks, numpnts(w_coef)) cmat
+	make /free /n=(npeaks) cvec
+	cmat = 0
+	cvec = 0
+
+	string hold = "00"
+	for (ipk=0; ipk < npeaks; ipk += 1)
+		hold += "0111"
+		cmat[ipk][2 + ipk * coef_per_peak] = -1
+	endfor
+
+	// prepare output
+	make /free /n=(npeaks * 2) /wave result_waves
+	string s_note
+	for (ipk = 0; ipk < npeaks; ipk += 1)
+		make /free /n=0 pk_int
+		adh5_setup_profile(source, pk_int, 1)
+		pk_int = nan
+		sprintf s_note, "AxisLabelD=peak %u integral", ipk+1
+		Note pk_int, s_note
+		sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * coef_per_peak]
+		Note pk_int, s_note
+		result_waves[ipk] = pk_int
+		
+		make /free /n=0 pk_sig
+		adh5_setup_profile(source, pk_sig, 1)
+		pk_sig = nan
+		sprintf s_note, "AxisLabelD=peak %u sigma", ipk+1
+		Note pk_sig, s_note
+		sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * coef_per_peak]
+		Note pk_sig, s_note
+		result_waves[ipk + npeaks] = pk_sig
+		
+		waveclear pk_int, pk_sig
+	endfor
+
+	// loop over angle scale
+	variable p0 = 0
+	variable p1 = dimsize(source, 1) - 1
+	variable pp
+	variable wmin
+	variable wmax
+	variable V_FitNumIters
+	variable V_FitError
+	
+	for (pp = p0; pp <= p1; pp += 1)
+		xprof = source[p][pp]
+		xprof_sig = max(sqrt(xprof), 1)
+		
+		// generate guess
+		wmin = wavemin(xprof)
+		wmax = wavemax(xprof)
+		w_coef[0] = wmin
+		w_coef[1] = 0
+
+		for (ipk=0; ipk < npeaks; ipk += 1)
+			w_coef[2 + ipk * coef_per_peak] = wmax - wmin
+		endfor
+		
+		V_FitError = 0
+		FuncFit /H=hold /Q /N /W=2 MultiVoigtLinBG_AO w_coef xprof[pl,ph] /C={cmat, cvec} /I=1 /W=xprof_sig[pl,ph]
+		wave w_sigma
+		
+		// retrieve results, leave them at nan if the fit did not converge
+		if (V_FitNumIters < 40)
+			for (ipk = 0; ipk < npeaks; ipk += 1)
+				wave val = result_waves[ipk]
+				wave sig = result_waves[ipk + npeaks]
+				val[pp] = max(w_coef[peak_coef[ipk]], 0)
+				sig[pp] = max(w_sigma[peak_coef[ipk]], 0)
+			endfor
+		endif
+	endfor
+
+	SetDataFolder orig_dfr
+	return result_waves
+end
+
+
+/// prompt for the voigt4_reduction parameters
+///
+///
+function prompt_voigt4_reduction(param)
+	string &param
+
+	variable rngl = NumberByKey("rngl", param, "=", ";")
+	variable rngh = NumberByKey("rngh", param, "=", ";")
+	variable pos1 = NumberByKey("pos1", param, "=", ";")
+	variable wid1 = NumberByKey("wid1", param, "=", ";")
+	variable shp1 = NumberByKey("shp1", param, "=", ";")
+	variable pos2 = NumberByKey("pos2", param, "=", ";")
+	variable wid2 = NumberByKey("wid2", param, "=", ";")
+	variable shp2 = NumberByKey("shp2", param, "=", ";")
+	variable pos3 = NumberByKey("pos3", param, "=", ";")
+	variable wid3 = NumberByKey("wid3", param, "=", ";")
+	variable shp3 = NumberByKey("shp3", param, "=", ";")
+	variable pos4 = NumberByKey("pos4", param, "=", ";")
+	variable wid4 = NumberByKey("wid4", param, "=", ";")
+	variable shp4 = NumberByKey("shp4", param, "=", ";")
+	variable npeaks = NumberByKey("npeaks", param, "=", ";")
+	variable dummy = nan
+
+	prompt rngl, "range low"
+	prompt rngh, "range high"
+	prompt pos1, "position 1"
+	prompt wid1, "width 1"
+	prompt shp1, "shape 1"
+	prompt pos2, "position 2"
+	prompt wid2, "width 2"
+	prompt shp2, "shape 2"
+	prompt pos3, "position 3"
+	prompt wid3, "width 3"
+	prompt shp3, "shape 3"
+	prompt pos4, "position 4"
+	prompt wid4, "width 4"
+	prompt shp4, "shape 4"
+	prompt npeaks, "number of peaks"
+	prompt dummy, "(not used)"
+	
+	doprompt "voigt4_reduction reduction parameters (1/2)", rngl, rngh, npeaks, dummy, pos1, pos2, pos3, pos4
+	if (v_flag == 0)
+		param = ReplaceNumberByKey("rngl", param, rngl, "=", ";")
+		param = ReplaceNumberByKey("rngh", param, rngh, "=", ";")
+		param = ReplaceNumberByKey("npeaks", param, npeaks, "=", ";")
+		param = ReplaceNumberByKey("pos1", param, pos1, "=", ";")
+		param = ReplaceNumberByKey("pos2", param, pos2, "=", ";")
+		param = ReplaceNumberByKey("pos3", param, pos3, "=", ";")
+		param = ReplaceNumberByKey("pos4", param, pos4, "=", ";")
+	
+		doprompt "voigt4_reduction reduction parameters (2/2)", wid1, shp1, wid2, shp2, wid3, shp3, wid4, shp4
+		if (v_flag == 0)
+			param = ReplaceNumberByKey("wid1", param, wid1, "=", ";")
+			param = ReplaceNumberByKey("shp1", param, shp1, "=", ";")
+			param = ReplaceNumberByKey("wid2", param, wid2, "=", ";")
+			param = ReplaceNumberByKey("shp2", param, shp2, "=", ";")
+			param = ReplaceNumberByKey("wid3", param, wid3, "=", ";")
+			param = ReplaceNumberByKey("shp3", param, shp3, "=", ";")
+			param = ReplaceNumberByKey("wid4", param, wid4, "=", ";")
+			param = ReplaceNumberByKey("shp4", param, shp4, "=", ";")
+		endif
+	endif
+	
+	return v_flag	
+end
+
+
 /// apply the Shockley_anglefit function to a single image
 ///
 /// useful for testing or manual processing

pearl/pearl-tools.ipf

@@ -1,3 +1,4 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3		// Use modern global access method and strict wave access.
 #pragma version = 1.00
 #pragma IgorVersion = 6.2

pearl/pearl-vector-operations.ipf

@@ -1,9 +1,10 @@
+#pragma TextEncoding = "UTF-8"
 #pragma rtGlobals=3
-#pragma version = 2.1
+#pragma version = 2.2
 #pragma IgorVersion = 6.1
 #pragma ModuleName = PearlVectorOperations
 
-// copyright (c) 2011-17 Paul Scherrer Institut
+// copyright (c) 2011-21 Paul Scherrer Institut
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 //   you may not use this file except in compliance with the License.
@@ -21,7 +22,7 @@
 ///
 /// @author matthias muntwiler, matthias.muntwiler@psi.ch
 ///
-/// @copyright 2011-17 Paul Scherrer Institut @n
+/// @copyright 2011-21 Paul Scherrer Institut @n
 /// Licensed under the Apache License, Version 2.0 (the "License"); @n
 ///   you may not use this file except in compliance with the License. @n
 ///   You may obtain a copy of the License at
@@ -176,16 +177,17 @@ function rotate_x_wave(inout, angle)
 	wave inout
 	variable angle
 
-	wave m_rotation_x = create_rotation_matrix_free()
+	wave m_rotation = create_rotation_matrix_free()
-	make /n=3/d/free w_temp_rotate_x
+	set_rotation_x(m_rotation, angle)
-	variable ivec, nvec
+
-	nvec = max(DimSize(inout, 1), 1)
+	duplicate /free inout, out
-	for (ivec = 0; ivec < nvec; ivec += 1)
+	out = 0
-		set_rotation_x(m_rotation_x, angle)
+	variable j
-		w_temp_rotate_x = inout[p][ivec]
+	for (j = 0; j < 3; j += 1)
-		matrixop /free w_temp_rotate_x_result = m_rotation_x x w_temp_rotate_x
+		out += m_rotation[p][j] * inout[j][q]
-		inout[0,2][ivec] = w_temp_rotate_x_result[p]
 	endfor
+	
+	inout = out
 end
 
 /// rotates a wave of 3-vectors about the y axis
@@ -204,16 +206,17 @@ function rotate_y_wave(inout, angle)
 	wave inout
 	variable angle
 
-	wave m_rotation_y = create_rotation_matrix_free()
+	wave m_rotation = create_rotation_matrix_free()
-	make /n=3/d/free w_temp_rotate_y
+	set_rotation_y(m_rotation, angle)
-	variable ivec, nvec
+
-	nvec = max(DimSize(inout, 1), 1)
+	duplicate /free inout, out
-	for (ivec = 0; ivec < nvec; ivec += 1)
+	out = 0
-		set_rotation_y(m_rotation_y, angle)
+	variable j
-		w_temp_rotate_y = inout[p][ivec]
+	for (j = 0; j < 3; j += 1)
-		matrixop /free w_temp_rotate_y_result = m_rotation_y x w_temp_rotate_y
+		out += m_rotation[p][j] * inout[j][q]
-		inout[0,2][ivec] = w_temp_rotate_y_result[p]
+	endfor
-	endfor	
+	
+	inout = out
 end
 
 /// rotates a wave of 3-vectors about the z axis
@@ -232,14 +235,15 @@ function rotate_z_wave(inout, angle)
 	wave inout
 	variable angle
 
-	wave m_rotation_z = create_rotation_matrix_free()
+	wave m_rotation = create_rotation_matrix_free()
-	make /n=3/d/free w_temp_rotate_z
+	set_rotation_z(m_rotation, angle)
-	variable ivec, nvec
+
-	nvec = max(DimSize(inout, 1), 1)
+	duplicate /free inout, out
-	for (ivec = 0; ivec < nvec; ivec += 1)
+	out = 0
-		set_rotation_z(m_rotation_z, angle)
+	variable j
-		w_temp_rotate_z = inout[p][ivec]
+	for (j = 0; j < 3; j += 1)
-		matrixop /free w_temp_rotate_z_result = m_rotation_z x w_temp_rotate_z
+		out += m_rotation[p][j] * inout[j][q]
-		inout[0,2][ivec] = w_temp_rotate_z_result[p]
+	endfor
-	endfor	
+	
+	inout = out
 end