#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.
#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-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.
//   You may obtain a copy of the License at
//   http:///www.apache.org/licenses/LICENSE-2.0
//
// Please acknowledge the use of this code.

/// @file
/// @brief processing and holographic mapping of angle scanned XPD data.
/// @ingroup ArpesPackage
///
/// the functions in this file map angle scanned data measured at PEARL
/// onto a hemispherical angle grid which is compatible with XPDplot.
/// the resulting data are in a canonical polar coordinate system
/// (normal emission <=> polar angle = 0, azimuthal axis right-handed)
/// which is anchored in the sample surface.
/// the orientation of polar graphs (phi = 0 at 3 o'clock) created by this procedure
/// corresponds to a top view of the sample surface at normal emission,
/// and the handle of the sample plate pointing to the left (phi = 180).
/// this is the canonical orientation of a spherical coordinate system 
/// where phi = 0 corresponds to the positive part of the x axis.
///
/// @note the orientation of the sample coordinate system has changed in version 1.6.
/// the change was necessary for compatibility with other data analysis software
/// and calculation programs.
/// @note data imported with version 1.5 and earlier, must be offset by 180 deg in phi
/// to be compatible with the new version.
/// data imported and displayed by the same code version will give the same picture
/// but with different azimuthal axis.
/// the new graph functions shows a warning if they are applied to code imported with earlier versions.
///
/// the measurement geometry is hard-coded but may be parametrized in the future.
/// the theta rotation axis is perpendicular to the scattering plane.
/// the angle dispersive axis of the analyser is parallel to the theta rotation axis.
/// the tilt rotation axis is in the scattering plane.
/// it rotates with theta.
/// at normal emission it is perpendicular to the axis of the lens stack of the analyser.
/// the phi rotation axis corresponds to the surface normal of the sample.
/// it rotates with theta and with tilt.
/// at normal emission it is parallel to the axis of the lens stack of the analyser.
///
/// coordinate transformations: (to be revised - v1.6)
///@verbatim
///     theta_sample = theta_manipulator - theta_offset
///     phi_sample = phi_manipulator - phi_offset
///@endverbatim
///
/// valid for theta_manipulator = normal emission only: (to be revised - v1.6)
///@verbatim
///     theta_sample = | -(tilt_manipulator - tilt_offset) |
///     phi_sample = 270 if tilt_manipulator - tilt_offset > 0
///     phi_sample = 90 if tilt_manipulator - tilt_offset < 0
///@endverbatim
///
/// @author matthias muntwiler, matthias.muntwiler@psi.ch
///
/// @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
///   http://www.apache.org/licenses/LICENSE-2.0
///
/// @version 1.8
/// canonical orientation of spherical coordinate system.
///

/// @namespace PearlAnglescanProcess
/// @brief processing and holographic mapping of angle scanned XPD data.
///
/// 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.
/// the function operates on 2D intensity data and manipulator coordinates at the same time.
///
/// @param[in,out]	strip			2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan.
///								the result is written to the original wave.
///
/// @param[in,out]	theta		1D data, manipulator scan.
///								the result is written to the original wave.
///
/// @param[in,out]	tilt			1D data, manipulator scan.
///								the result is written to the original wave.
///
/// @param[in,out]	phi			1D data, manipulator scan.
///								the result is written to the original wave.
///
/// @param[in]		qlo			point index of first frame to delete.
///
/// @param[in]		qhi			point index of last frame to delete.
///								qhi must be greater or equal than qlo.
///
function strip_delete_frames(strip, qlo, qhi, theta, tilt, phi)
	wave strip // 2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan
	variable qlo
	variable qhi
	wave theta
	wave tilt
	wave phi
	
	if (qlo > qhi)
		return -1
	endif
	
	// source indices
	variable snx = dimsize(strip, 0)
	variable sny = dimsize(strip, 1)
	variable sq1lo = 0
	variable sq1hi = max(qlo-1, 0)
	variable sq2lo = min(qhi+1, sny - 1)
	variable sq2hi = dimsize(strip, 1) - 1
	
	// dest indices
	variable dnx = snx
	variable dny = sny - (sq2lo - sq1hi + 1)
	variable dq1lo = 0
	variable dq1hi = sq1hi
	variable dq2lo = dq1hi + 1
	variable dq2hi = dny - 1
	variable q1ofs = sq1lo - dq1lo
	variable q2ofs = sq2lo - dq2lo
	
	duplicate /free strip, strip_copy
	redimension /n=(dnx,dny) strip
	strip[][dq1lo,dq1hi] = strip_copy[p][q + q1ofs]
	strip[][dq2lo,dq2hi] = strip_copy[p][q + q2ofs]
	
	duplicate /free theta, theta_copy
	redimension /n=(dny) theta
	theta[dq1lo,dq1hi] = theta_copy[p + q1ofs]
	theta[dq2lo,dq2hi] = theta_copy[p + q2ofs]
	
	duplicate /free tilt, tilt_copy
	redimension /n=(dny) tilt
	tilt[dq1lo,dq1hi] = tilt_copy[p + q1ofs]
	tilt[dq2lo,dq2hi] = tilt_copy[p + q2ofs]
	
	duplicate /free phi, phi_copy
	redimension /n=(dny) phi
	phi[dq1lo,dq1hi] = phi_copy[p + q1ofs]
	phi[dq2lo,dq2hi] = phi_copy[p + q2ofs]
	
	return 0
end

/// divide the strip by the average X distribution.
///
/// this is a simple way to remove the effect of the angle-dependence of the analyser transmission function.
/// the strip is normalized in place, previous data is overwritten.
///
/// the function can handle sparse NaNs.
///
/// @param[in,out] strip 2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan
///
/// @param[in] smooth_method smoothing method
/// @arg 0 none
/// @arg 1 binomial, see Igor's Smooth operation
/// @arg 2 boxcar, see Igor's Smooth operation
/// @arg 3 scienta_ang_transm() function fit
/// @arg 4 (default) LOESS smoothing, see Igor's Loess operation
///
/// @param[in] smooth_factor num parameter of Igor's Smooth operation.
/// the default value depends on smooth_method.
/// it is 0.5 for LOESS smoothing, 2 otherwise.
///
/// @param[in] check enable output of intermediate results
/// @arg 0 (default) don't create additional waves
/// @arg 1 create check waves in the current folder
/// @arg 2 calculate check waves only, do not modify strip
///
/// @return if check waves are enabled, the following waves are created (overwritten if existing):
/// @arg check_dist average X distribution
/// @arg check_smoo smoothed distribution used to normalize the strip
///
function normalize_strip_x(strip, [smooth_method, smooth_factor, check])
	wave strip
	variable smooth_method
	variable smooth_factor
	variable check
	
	if (ParamIsDefault(smooth_method))
		smooth_method = 4
	endif
	if (ParamIsDefault(smooth_factor))
		switch(smooth_method)
			case 4:
				smooth_factor = 0.5
				break
			default:
				smooth_factor = 2
		endswitch
	endif
	if (ParamIsDefault(check))
		check = 0
	endif
	
	// average over all scan positions
	wave raw_dist = ad_profile_x(strip, -inf, inf, "")
	
	// remove nans
	extract /free /indx raw_dist, clean_index, numtype(raw_dist) == 0
	duplicate /free raw_dist, dist, dist_x
	redimension /n=(numpnts(clean_index)) dist, dist_x
	dist = raw_dist[clean_index[p]]
	dist_x = pnt2x(raw_dist, clean_index[p])
	variable div = mean(dist)
	dist /= div

	if (check)
		duplicate /o raw_dist, check_dist
		check_dist = numtype(raw_dist) == 0 ? interp(x, dist_x, dist) : nan
	endif
	
	// smooth distribution function
	switch(smooth_method)
		case 1:
			Smooth /B /E=3 smooth_factor, dist
			break
		case 2:
			Smooth /E=3 smooth_factor, dist
			break
		case 3:
			make /n=1 /d /free fit_params
			fit_scienta_ang_transm(raw_dist, fit_params)
			duplicate /free raw_dist, dist, dist_x
			dist_x = x
			dist = scienta_ang_transm(fit_params, x)
			break
		case 4:
			loess /smth=(smooth_factor) srcWave=dist, factors={dist_x}
			break
	endswitch
	
	if (check)
		duplicate /o raw_dist, check_smoo
		check_smoo = interp(x, dist_x, dist)
	endif

	// divide
	if (check != 2)
		strip /= interp(x, dist_x, dist)
	endif
end

/// divide the strip by a sine function in phi (wobble correction).
///
/// the sine function is a curve fit to the intensity integrated over detector angle
/// with a period of 360&deg;.
///
/// this normalization may be useful if the intensity varies with a 360&deg; periodicity in the azimuthal angle, 
/// e.g. due to misalignment of the surface normal and the azimuthal rotation axis of the manipulator (wobble).
/// note, however, that this function does not correct other effects of wobble such as angle shifts.
///
/// the strip is normalized in place, previous data is overwritten.
///
/// @param[in,out]	strip				2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan
/// @param[in]		theta			polar manipulator angle.
/// @param[in]		phi				azimuthal manipulator angle, arbitrary offset.
/// @param[in]		theta_offset		theta value corresponding to normal emission (default 0).
/// @param[in]		theta_range		maximum (offset corrected) theta to consider in the sine fit (default 10).
///
/// @param check enable output of intermediate results
///	@arg 0 (default)  don't create additional waves
/// 	@arg 1 create check waves in the current folder
///	@arg 2 calculate check waves only, do not modify strip
///
/// @return if check waves are enabled, the following waves are created (overwritten if existing):
///	@arg check_dist average theta distribution
///	@arg check_smoo smoothed distribution used to normalize the strip
///
function normalize_strip_phi(strip, theta, phi, [theta_offset, theta_range, check])
	wave strip
	wave theta
	wave phi
	variable theta_offset
	variable theta_range
	variable check
	
	if (ParamIsDefault(check))
		check = 0
	endif
	if (ParamIsDefault(theta_offset))
		theta_offset = 0
	endif
	if (ParamIsDefault(theta_range))
		theta_offset = 10
	endif

	// average over analyser angles
	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
	duplicate /free theta, theta_int
	theta_int = theta - theta_offset
	duplicate /free phi, phi_int
	setscale /p x phi_int[0], phi_int[1] - phi_int[0], waveunits(phi, -1), dist, dist_smoo

	extract /free /indx dist, red_idx, theta_int < theta_range
	duplicate /free red_idx, red_dist, red_phi
	red_dist = dist[red_idx]
	red_phi = phi_int[red_idx]
	
	variable wavg = mean(red_dist)
	make /n=4 /d /free coef
	coef[0] = {wavg, wavg/100, pi/180, 0}
	CurveFit /q /h="0010" /g /w=2 sin, kwcWave=coef, red_dist /x=red_phi
	dist_smoo = coef[0] + coef[1] * sin(coef[2] * phi_int[p] + coef[3])
	
	// divide
	if (check != 2)
		strip = strip / dist_smoo[q] * coef[0]
	endif
	
	// check
	if (check)
		duplicate /o dist, check_dist
		duplicate /o dist_smoo, check_smoo
		setscale /p x dimoffset(dist,0), dimdelta(dist,0), waveunits(dist,0), check_dist, check_smoo
	endif
end

/// divide the strip by the average polar distribution.
///
/// this is a simple way to remove the polar angle dependence.
/// the strip is normalized in place, previous data is overwritten.
/// @param[in,out] strip 2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan
/// @param[in] theta polar manipulator angle, 0 = normal emission, 90 = grazing emission
/// @param[in] theta_offset
/// @param[in] smooth_method smoothing method
///	@arg 0 none
///	@arg 1 binomial (requires monotonic theta), see Igor's Smooth operation
///	@arg 2 boxcar (requires monotonic theta), see Igor's Smooth operation
///	@arg 3 polynomial fit per slice
///	@arg 4 (default) Loess, see Igor's Loess operation
///
/// caution: binomial and boxcar smoothing are not aware of theta.
/// this may give unpredictable results if theta is non-monotonic.
///
/// @param[in] smooth_factor smoothing parameter, depends on smooth_method
///	@arg binomial/boxcar: see Igor's Smooth operation
///	@arg loess: see Igor's Loess operation, 0 <= smooth_factor <= 1, default 0.5
///	@arg polynomial fit: polynomial degree, 1 = linear (default), 2 = quadratic
///
/// @param check enable output of intermediate results
///	@arg 0 (default)  don't create additional waves
/// 	@arg 1 create check waves in the current folder
///	@arg 2 calculate check waves only, do not modify strip
///
/// @return if check waves are enabled, the following waves are created (overwritten if existing):
///	@arg check_dist average theta distribution
///	@arg check_smoo smoothed distribution used to normalize the strip
///
function normalize_strip_theta(strip, theta, [theta_offset, smooth_method, smooth_factor, check])
	wave strip
	wave theta
	variable theta_offset
	variable smooth_method
	variable smooth_factor
	variable check
	
	if (ParamIsDefault(check))
		check = 0
	endif
	if (ParamIsDefault(theta_offset))
		theta_offset = 0
	endif
	if (ParamIsDefault(smooth_method))
		smooth_method = 4
	endif
	if (ParamIsDefault(smooth_factor))
		smooth_factor = 0.5
	endif

	// average over analyser angles
	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
	duplicate /free theta, theta_int
	theta_int = theta - theta_offset
	setscale /p x theta_int[0], theta_int[1] - theta_int[0], waveunits(theta,-1), dist, dist_smoo
	variable nx = dimsize(strip, 0)
	variable ix
	
	switch(smooth_method)
		case 1:
			Smooth /B /E=3 smooth_factor, dist_smoo
			break
		case 2:
			Smooth /E=3 smooth_factor, dist_smoo
			break
		case 4:
			loess /dest=dist_smoo /smth=(smooth_factor) srcWave=dist, factors={theta_int}
			break
		case 3:
			for (ix = 0; ix < nx; ix += 1)
				dist = strip[ix][p]
				if (smooth_factor > 1)
					CurveFit /nthr=0 /q /w=2 poly smooth_factor+1, dist /x=theta_int /d=dist_smoo
				else
					CurveFit /nthr=0 /q /w=2 line, dist /x=theta_int /d=dist_smoo
				endif
				strip[ix,ix][] /= dist_smoo[q]
			endfor
			dist_smoo = 1
			break
	endswitch
	
	// divide
	if (check != 2)
		strip /= dist_smoo[q]
	endif
	
	// check
	if (check)
		duplicate /o dist, check_dist
		duplicate /o dist_smoo, check_smoo
		setscale /p x dimoffset(dist,0), dimdelta(dist,0), waveunits(dist,0), check_dist, check_smoo
	endif
end

/// divide the strip by a smooth polar-azimuthal distribution.
///
/// this is a simple way to remove the polar angle dependence.
/// in contrast to @ref normalize_strip_theta this function also removes a smooth variation over azimuthal angles.
///
/// the strip is normalized in place, previous data is overwritten.
///
/// @warning experimental. this function is under development.
///
/// @param[in,out] strip 2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan
/// @param[in] theta polar manipulator angle, 0 = normal emission, 90 = grazing emission
/// @param[in] phi azimuthal manipulator angle.
/// @param[in] theta_offset
/// @param[in] smooth_method smoothing method
///	@arg 0 none
///	@arg 4 (default) Loess, see Igor's Loess operation
///
/// @param[in] smooth_factor smoothing parameter, depends on smooth_method
///	@arg loess: see Igor's Loess operation, 0 <= smooth_factor <= 1, default 0.5
///
/// @param check enable output of intermediate results
///	@arg 0 (default)  don't create additional waves
/// 	@arg 1 create check waves in the current folder
///	@arg 2 calculate check waves only, do not modify strip
///
/// @return if check waves are enabled, the following waves are created (overwritten if existing):
///	@arg check_dist average theta distribution
///	@arg check_smoo smoothed distribution used to normalize the strip
///
function normalize_strip_thetaphi(strip, theta, phi, [theta_offset, smooth_method, smooth_factor, check])
	wave strip
	wave theta
	wave phi
	variable theta_offset
	variable smooth_method
	variable smooth_factor
	variable check
	
	if (ParamIsDefault(check))
		check = 0
	endif
	if (ParamIsDefault(theta_offset))
		theta_offset = 0
	endif
	if (ParamIsDefault(smooth_method))
		smooth_method = 4
	endif
	if (ParamIsDefault(smooth_factor))
		smooth_factor = 0.5
	endif

	// average over analyser angles
	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
	duplicate /free theta, theta_int
	theta_int = theta - theta_offset
	setscale /p x theta_int[0], theta_int[1] - theta_int[0], waveunits(theta,-1), dist, dist_smoo
	variable nx = dimsize(strip, 0)
	variable ix
	
	switch(smooth_method)
		case 4:
			loess /dest=dist_smoo /smth=(smooth_factor) srcWave=dist, factors={theta_int, phi}
			break
		default:
			abort "smooth method not supported"
	endswitch
	
	// divide
	if (check != 2)
		strip /= dist_smoo[q]
	endif
	
	// check
	if (check)
		duplicate /o dist, check_dist
		duplicate /o dist_smoo, check_smoo
		setscale /p x dimoffset(dist,0), dimdelta(dist,0), waveunits(dist,0), check_dist, check_smoo
	endif
end

/// divide the strip piecewise by a smooth polar distribution.
///
/// @warning experimental. this function is under development.
///
///
function normalize_strip_theta_scans(strip, theta, [theta_offset, smooth_method, smooth_factor, check])
	wave strip
	wave theta
	variable theta_offset
	variable smooth_method
	variable smooth_factor
	variable check
	
	if (ParamIsDefault(check))
		check = 0
	endif
	if (ParamIsDefault(theta_offset))
		theta_offset = 0
	endif
	if (ParamIsDefault(smooth_method))
		smooth_method = 4
	endif
	if (ParamIsDefault(smooth_factor))
		smooth_factor = 0.5
	endif

	// average over analyser angles
	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
	duplicate /free theta, theta_int
	theta_int = theta - theta_offset
	setscale /p x theta_int[0], theta_int[1] - theta_int[0], waveunits(theta,-1), dist, dist_smoo

	// analyse scanning scheme
	duplicate /free theta_int, d1_theta, d2_theta
	Differentiate /METH=2 theta_int /D=d1_theta
	Differentiate /METH=2 d1_theta /D=d2_theta
	d2_theta =  abs(d2_theta)
	make /free w_levels
	FindLevels /edge=1 /p /q /d=w_levels d2_theta, 0.1
	if (v_flag != 1)
		abort "unrecognized scanning scheme"
	endif
	w_levels = ceil(w_levels)
	InsertPoints 0, 1, w_levels
	w_levels[0] = 0
	InsertPoints numpnts(w_levels), 1, w_levels
	w_levels[numpnts(w_levels)-1] = numpnts(theta_int)
	
	variable n_scans = numpnts(w_levels) - 1
	variable i_scan
	variable p1, p2
	for (i_scan = 0; i_scan < n_scans; i_scan += 1)
		p1 = w_levels[i_scan]
		p2 = w_levels[i_scan+1] - 1
		duplicate /free /r=[p1, p2] dist, dist_piece, smooth_piece
		duplicate /free /r=[p1, p2] theta_int, theta_piece
		switch(smooth_method)
			case 4:
				loess /dest=smooth_piece /smth=(smooth_factor) srcWave=dist_piece, factors={theta_piece}
				break
			default:
				abort "smooth method not supported"
		endswitch
		dist_smoo[p1, p2] = smooth_piece[p - p1]
	endfor
	
	// divide
	if (check != 2)
		strip /= dist_smoo[q]
	endif

	// check
	if (check)
		duplicate /o dist, check_dist
		duplicate /o dist_smoo, check_smoo
		setscale /p x dimoffset(dist,0), dimdelta(dist,0), waveunits(dist,0), check_dist, check_smoo
	endif
end

/// divide the strip by a two-dimensional normalization function.
///
/// @warning experimental. this function is under development.
///
/// @param check enable output of intermediate results
///	@arg 0 (default)  don't create additional waves
/// 	@arg 1 create check waves in the current folder
///	@arg 2 calculate check waves only, do not modify strip
///
/// @return if check waves are enabled, the following waves are created (overwritten if existing):
///	@arg check_dist average theta distribution
///	@arg check_smoo smoothed distribution used to normalize the strip
///
function normalize_strip_2d(strip, theta, [theta_offset, smooth_method, smooth_factor, check])
	wave strip
	wave theta
	variable theta_offset
	variable smooth_method
	variable smooth_factor
	variable check
	
	if (ParamIsDefault(check))
		check = 0
	endif
	if (ParamIsDefault(theta_offset))
		theta_offset = 0
	endif
	if (ParamIsDefault(smooth_method))
		smooth_method = 4
	endif
	if (ParamIsDefault(smooth_factor))
		smooth_factor = 0.5
	endif

	variable nx = dimsize(strip, 0)
	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
	
	switch(smooth_method)
		case 4:
			loess /dest=dist_smoo /smth=(smooth_factor) srcWave=dist, factors={alpha_int, theta_int}
			redimension /n=(nx, ny) dist_smoo
			break
		default:
			Abort "undefined smooth method"
			break
	endswitch
	
	// divide
	if (check != 2)
		strip /= dist_smoo
	endif
	
	// check
	if (check)
		//duplicate /o dist, check_dist
		duplicate /o dist_smoo, check_smoo
	endif
end

/// crop a strip at the sides.
///
/// the strip is cropped in place, data outside the region of interest is lost.
///
/// @param[in,out] 	strip		2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan
/// @param[in]		xlo		lowest analyser angle to keep (will be rounded to nearest existing point)
/// @param[in]		xhi		highest analyser angle to keep (will be rounded to nearest existing point)
///
/// @remark cropping should be done after smoothing and normalization operations to reduce artefacts.
function crop_strip(strip, xlo, xhi)
	wave strip
	variable xlo
	variable xhi
	
	variable plo = round((xlo - dimoffset(strip, 0)) / dimdelta(strip, 0))
	variable phi = round((xhi - dimoffset(strip, 0)) / dimdelta(strip, 0))
	xlo = plo * dimdelta(strip, 0) + dimoffset(strip, 0)
	xhi = phi * dimdelta(strip, 0) + dimoffset(strip, 0)
	variable nx = phi - plo + 1
	variable ny = dimsize(strip, 1)
	
	duplicate /free strip, strip_copy
	redimension /n=(nx,ny) strip
	strip = strip_copy[p + plo][q]
	setscale /i x xlo, xhi, waveunits(strip, 0), strip
end

/// crop a strip in theta.
///
/// the strip is cropped in place, data outside the region of interest is lost.
///
/// @param[in,out] 	strip		2D data, X-axis = analyser angle, Y-axis = arbitrary manipulator scan
/// @param[in]		ylo			lowest polar angle to keep (will be rounded to nearest existing point)
/// @param[in]		yhi			highest polar angle to keep (will be rounded to nearest existing point)
/// @param[in,out]	theta		polar angle along the Y dimension of strip.
///									this wave is modified: cropped rows are deleted.
/// @param[in,out]	tilt		tilt angle along the Y dimension of strip.
///									this wave is modified: cropped rows are deleted.
/// @param[in,out]	phi			azimuthal angle along the Y dimension of strip.
///									this wave is modified: cropped rows are deleted.
///
function crop_strip_theta(strip, theta_lo, theta_hi, theta, tilt, phi)
	wave strip
	variable theta_lo
	variable theta_hi
	wave theta
	wave tilt
	wave phi
	
	extract /indx /free theta, idx, (theta >= theta_lo) && (theta <= theta_hi)
	variable nx = dimsize(strip, 0)
	variable ny = numpnts(idx)
	
	theta[0, ny-1] = theta[idx]
	tilt[0, ny-1] = tilt[idx]
	phi[0, ny-1] = phi[idx]
	redimension /n=(ny) theta, tilt, phi
	
	duplicate /free strip, strip_copy
	redimension /n=(nx,ny) strip
	strip = strip_copy[p][idx[q]]
end

/// create a pizza plot from a measured (energy-integrated) data strip
///
/// accepts angle-scan data as returned by adh5_load_reduced(),
/// maps them onto a hemispherical scan grid,
/// and displays a polar graph.
///
/// @param	data			2D intensity wave, see requirements above
///	@arg X-axis analyser angle
///	@arg Y-axis manipulator scan.
///	no specific ordering required.
///	manipulator angle waves (ManipulatorTheta, ManipulatorTilt, ManipulatorPhi)
///	must be in the subfolder _attr_ below the data wave.
///
/// @param	nickname	nick name for output data
///	@arg in default mode, this will become the name of a child folder containing the output.
///	@arg in XPDplot mode, this will become a prefix of the generated data in the root folder.
///
/// @param	theta_offset	manipulator theta angle corresponding to normal emission.
///                                    the offset is subtracted from the ManipulatorTheta wave before processing.
///
/// @param	tilt_offset	manipulator tilt angle corresponding to normal emission
///                                    the offset is subtracted from the ManipulatorTilt wave before processing.
///
/// @param	phi_offset	manipulator phi angle corresponding to phi_result = 0
///                                    the offset is subtracted from the ManipulatorPhi wave before processing.
///
/// @param	npolar		number of polar angles, determines polar and azimuthal step size.
///						default = 91 (1 degree steps)
///
/// @param	folding		rotational averaging, default = 1
///
/// @param	nograph		display a new graph window?
///	@arg	0			(default) display a new polar graph
///	@arg	1			don't display a new graph
///
/// @param	xpdplot		XPDplot compatibility
///	@arg	0			(default) create waves in child folder $nickname
///	@arg	1			create waves in root folder (compatible with XPDplot)
///
/// @attention if you modify the structure of the data wave, e.g. delete some angles,
/// this function cannot be used because the manipulator settings do not correspond to the original manipulator waves!
/// instead, create your own manipulator waves and use pizza_service_2().
function pizza_service(data, nickname, theta_offset, tilt_offset, phi_offset, [npolar, nograph, folding, xpdplot])
	wave data
	string nickname
	variable theta_offset
	variable tilt_offset
	variable phi_offset
	variable npolar
	variable nograph
	variable folding
	variable xpdplot
		
	if (ParamIsDefault(npolar))
		npolar = 91
	endif
	if (ParamIsDefault(nograph))
		nograph = 0
	endif
	if (ParamIsDefault(folding))
		folding = 1
	endif
	if (ParamIsDefault(xpdplot))
		xpdplot = 0
	endif
	
	// sort out data folder structure
	dfref saveDF = GetDataFolderDFR()
	dfref dataDF = GetWavesDataFolderDFR(data)
	setdatafolder dataDF
	if (DataFolderExists(":attr"))
		setdatafolder :attr
	endif
	dfref attrDF = GetDataFolderDFR()
	
	wave /sdfr=attrDF ManipulatorTheta
	wave /sdfr=attrDF ManipulatorTilt
	wave /sdfr=attrDF ManipulatorPhi

	if ((dimsize(ManipulatorTheta, 0) != dimsize(data, 1)) || (dimsize(ManipulatorTilt, 0) != dimsize(data, 1)) || (dimsize(ManipulatorPhi, 0) != dimsize(data, 1)))
		Abort "Warning: The dimension size of the manipulator waves does not match the Y dimension of the data wave!\rIf you restructured the data wave, please use pizza_service_2 with properly scaled manipulator waves."
	endif

	duplicate /free ManipulatorTheta, m_theta
	duplicate /free ManipulatorTilt, m_tilt
	duplicate /free ManipulatorPhi, m_phi

	m_theta -= theta_offset
	m_tilt -= tilt_offset
	m_phi -= phi_offset
	
	pizza_service_2(data, nickname, m_theta, m_tilt, m_phi, npolar=npolar, nograph=nograph, folding=folding, xpdplot=xpdplot)
	
	setdatafolder saveDF
end

/// create a pizza plot from a measured (energy-integrated) data strip
///
/// accepts angle-scan data as returned by adh5_load_reduced(),
/// maps them onto a hemispherical scan grid,
/// and displays a polar graph.
///
/// the behaviour of this function is the same as pizza_service()
/// except that the manipulator waves are specified explicitly.
///
/// @param	data			2D intensity wave, see requirements above
///	@arg X-axis analyser angle
///	@arg Y-axis manipulator scan.
///	no specific ordering required.
///	manipulator angle waves (ManipulatorTheta, ManipulatorTilt, ManipulatorPhi)
///	must be in the subfolder _attr_ below the data wave.
///
/// @param	nickname	nick name for output data
///	@arg in default mode, this will become the name of a child folder containing the output.
///	@arg in XPDplot mode, this will become a prefix of the generated data in the root folder.
///
/// @param	m_theta		manipulator theta angles, 0 = normal emission. size = dimsize(data, 1)
///
/// @param	m_tilt		manipulator tilt angles, 0 = normal emission. size = dimsize(data, 1)
///
/// @param	m_phi		manipulator phi angles, 0 = azimuthal origin. size = dimsize(data, 1)
///
/// @param	npolar		number of polar angles, determines polar and azimuthal step size.
///						default = 91 (1 degree steps)
///
/// @param	folding		rotational averaging, default = 1
///
/// @param	nograph		display a new graph window?
///	@arg	0			(default) display a new polar graph
///	@arg	1			don't display a new graph
///
/// @param	xpdplot		XPDplot compatibility
///	@arg	0			(default) create waves in child folder $nickname
///	@arg	1			create waves in root folder (compatible with XPDplot)
///
function pizza_service_2(data, nickname, m_theta, m_tilt, m_phi, [npolar, nograph, folding, xpdplot])
	wave data
	string nickname
	wave m_theta
	wave m_tilt
	wave m_phi
	variable npolar
	variable nograph
	variable folding
	variable xpdplot
		
	if (ParamIsDefault(npolar))
		npolar = 91
	endif
	if (ParamIsDefault(nograph))
		nograph = 0
	endif
	if (ParamIsDefault(folding))
		folding = 1
	endif
	if (ParamIsDefault(xpdplot))
		xpdplot = 0
	endif
	
	if ((dimsize(m_theta, 0) != dimsize(data, 1)) || (dimsize(m_tilt, 0) != dimsize(data, 1)) || (dimsize(m_phi, 0) != dimsize(data, 1)))
		Abort "Warning: The dimension size of the manipulator waves does not match the Y dimension of the data wave!"
	endif
	
	string graphname = "graph_" + nickname
	string outprefix = nickname
	
	// sort out data folder structure
	dfref saveDF = GetDataFolderDFR()
	dfref dataDF = GetWavesDataFolderDFR(data)
	setdatafolder dataDF

	if (xpdplot)
		setdatafolder root:
		outprefix = nickname
	else
		setdatafolder dataDF
		newdatafolder /s/o $nickname
		outprefix = ""
	endif
	dfref destDF = GetDataFolderDFR()
	
	// performance monitoring
	variable timerRefNum
	variable /g pol_perf_secs
	timerRefNum = startMSTimer
	
	duplicate /free m_tilt, corr_tilt
	duplicate /free m_phi, corr_phi
	corr_tilt = -m_tilt	// checked 140702
	corr_phi = m_phi	// checked 140702
	
	make /n=1/d/free d_polar, d_azi
	
	convert_angles_ttpd2polar(m_theta, corr_tilt, corr_phi, data, d_polar, d_azi)
	d_azi += 180	// changed 151030 (v1.6)
	make_hemi_grid(npolar, outprefix, xpdplot=xpdplot)
	variable ifold
	for (ifold = 0; ifold < folding; ifold += 1)
		d_azi = d_azi >= 360 ? d_azi - 360 : d_azi
		hemi_add_anglescan(outprefix, data, d_polar, d_azi)
		d_azi += 360 / folding
	endfor
	
	// normalize folding
	if (strlen(outprefix))
		string s_prefix = outprefix + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	if (folding > 1)
		wave values = $s_int
		values /= folding
	endif
	
	if (!nograph)
		display_hemi_scan(outprefix, graphname = graphname)
	endif
	
	if (timerRefNum >= 0)
		pol_perf_secs = stopMSTimer(timerRefNum) / 1e6
	endif

	setdatafolder saveDF
end

/// calculate and display the line seen by the analyser for a specific emission angle
///
/// this can be used to compare to an hemispherical plot and check the manipulator angle.
///
/// @param	theta		manipulator theta angle
/// @param	tilt			manipulator tilt angle
/// @param	phi			manipulator phi angle
///
/// @param	theta_offset	manipulator theta angle corresponding to normal emission
/// @param	tilt_offset	manipulator tilt angle corresponding to normal emission
/// @param	phi_offset	manipulator phi angle corresponding to phi_result = 0
///
/// @param	npolar		number of polar angles, determines polar and azimuthal step size.
///						default = 91 (1 degree steps)
///
/// @param	nograph		display a new graph window?
///	@arg	0			(default) display a new polar graph
///	@arg	1			don't display a new graph
///
/// @param	xpdplot		XPDplot compatibility
///	@arg	0			(default) create waves in child folder $nickname
///	@arg	1			create waves in root folder (compatible with XPDplot)
///
/// @remark the function creates angle scan data under the nickname _analyser_.
///
function show_analyser_line(theta, tilt, phi, theta_offset, tilt_offset, phi_offset, [npolar, nograph, xpdplot])
	variable theta
	variable tilt
	variable phi
	variable theta_offset
	variable tilt_offset
	variable phi_offset
	variable npolar
	variable nograph
	variable xpdplot
		
	string nickname = "analyser"

	if (ParamIsDefault(npolar))
		npolar = 91
	endif
	if (ParamIsDefault(nograph))
		nograph = 0
	endif
	if (ParamIsDefault(xpdplot))
		xpdplot = 0
	endif
	string graphname = "graph_" + nickname
	string outprefix = nickname
	
	// sort out data folder structure
	dfref saveDF = GetDataFolderDFR()
	dfref dataDF = saveDF
	if (xpdplot)
		setdatafolder root:
		outprefix = nickname
	else
		setdatafolder dataDF
		newdatafolder /s/o $nickname
		outprefix = ""
	endif
	dfref destDF = GetDataFolderDFR()
	
	make /n=1 /free m_theta
	make /n=1 /free m_tilt
	make /n=1 /free m_phi
	m_theta = theta - theta_offset
	m_tilt = tilt - tilt_offset
	m_tilt *= -1 // checked 140702
	m_phi = phi - phi_offset
	//m_phi *= -1 // checked 140702
	
	make /n=60 /free data
	setscale /i x -30, 30, data
	data = x
	make /n=1/d/free d_polar, d_azi
	
	convert_angles_ttpa2polar(m_theta, m_tilt, m_phi, data, d_polar, d_azi)
	d_azi += 180	// changed 151030 (v1.6)
	d_azi = d_azi >= 360 ? d_azi - 360 : d_azi
	make_hemi_grid(npolar, outprefix, xpdplot=xpdplot)
	hemi_add_anglescan(outprefix, data, d_polar, d_azi)
	
	if (!nograph)
		display_hemi_scan(outprefix, graphname = graphname)
	endif
	
	setdatafolder saveDF
end

/// convert angles from TTPA (theta-tilt-phi-analyser) scheme to polar coordinates.
///
/// similar to convert_angles_ttpa2polar()
/// but reads the analyser angles from the X scale of data

function convert_angles_ttpd2polar(theta, tilt, phi, data, polar, azi)
	wave theta, tilt, phi // see convert_angles_ttpa2polar
	wave data // in, 1D or 2D
		// X-scale must be set to analyser angle scale
	wave polar, azi // see convert_angles_ttpa2polar
	
	make /n=(dimsize(data, 0)) /d /free ana
	setscale /p x dimoffset(data, 0), dimdelta(data, 0), waveunits(data, 0), ana
	ana = x
	convert_angles_ttpa2polar(theta, tilt, phi, ana, polar, azi)
end

/// convert angles from TTPA (theta-tilt-phi-analyser) scheme to polar coordinates.
///
/// the angles are in the manipulator coordinate system.
///
/// @param[in]	theta		offset-corrected theta angle, normal emission = 0, grazing emission = 90.
///							one dimensional wave.
///
/// @param[in]	tilt			offset-corrected tilt angle, normal emission = 0
///							same dimension size and scale as theta
///
/// @param[in]	phi			phi angle, range -360 < phi < +360
///							offset correction is optional as long as the angles lie in the accepted range.
///							same dimension size and scale as theta
///
/// @param[in]	analyser		analyser angle scale corresponding to the slices scale of Scienta.
///							one dimensional wave.
///							this values are constant regardless of manipulator angle.
///
/// @param[out]	polar		wave to receive the polar coordinates.
///
/// @param[out]	azi			wave to receive the azimuthal coordinates.
///
/// for the output parameters polar and azi, you need to pass in existing numeric waves.
/// dimension size does not matter, the waves are redimensioned by the function
/// so that they have the same dimensions as the intensity data set.
/// X dimension = analyser scale, Y dimension = manipulator scan.
///
function convert_angles_ttpa2polar(theta, tilt, phi, analyser, polar, azi)
	wave theta
	wave tilt
	wave phi
	wave analyser
	wave polar, azi
		
	variable nn = numpnts(theta)
	variable na = numpnts(analyser)
	redimension /n=(na, nn) polar, azi

	variable radius = 1 // don't need to specify - everything is scalable
	
	// step 1: calculate cartesian detection vectors at normal emission
	// this is simply a polar-cartesian mapping, independent of the manipulator
	// phi=0 is in the polar rotation plane
	make /n=(3,na) /d /free w_orig_polar, w_orig_cart, w_rot_cart, w_rot_polar
	w_orig_polar[0][] = radius
	w_orig_polar[1][] = analyser[q]
	w_orig_polar[2][] = 0
	polar2cart_wave(w_orig_polar, w_orig_cart)
	// if the angle-dispersive axis was horizontal, we'd need to rotate the detector
	//rotate_z_wave(w_orig_cart, 90)

	variable ii
	for (ii = 0; ii < nn; ii += 1)
		// step 2: rotate the detection vectors according to the manipulator angles
		// the order of rotations is important because we rotate about fixed axes
		// y-axis = tilt rotation axis
		// x-axis = polar rotation axis
		// z-axis = normal emission = azimuthal rotation axis
		w_rot_cart = w_orig_cart
		rotate_y_wave(w_rot_cart, -tilt[ii])
		rotate_x_wave(w_rot_cart, -theta[ii])
		rotate_z_wave(w_rot_cart, -phi[ii] - 90)
		// map the vectors back to the sample coordinate system
		cart2polar_wave(w_rot_cart, w_rot_polar)
		// copy to output
		polar[][ii] = w_rot_polar[1][p]
		azi[][ii] = w_rot_polar[2][p]
	endfor
end

static function line_average(source, dest)
	// is this function used?
	wave source
	wave dest
	
	variable ii
	variable nn = dimsize(source, 1)
	make /n=(dimsize(source, 0))/d/free line
	for (ii = 0; ii < nn; ii += 1)
		line = source[p][ii]
		wavestats /q line
		dest[][ii] = line[p] / v_max
	endfor
end

/// calculate the number of phis for a given theta
///
/// adapted from XPDplot 8.03
static function calc_nth(Theta_st, Theta_in, th, Phi_ran, Phi_ref, Holomode)
	Variable Theta_st, Theta_in, th, Phi_ran, Phi_ref
       String Holomode
	Variable The_step
	Variable deg2rad=0.01745329
	 	
	if ( cmpstr(Holomode, "Stereographic") == 0)
	   The_step =trunc( Phi_ran*sin(th*deg2rad)*Phi_ref/Theta_st )
	   if(th==90)
	   	The_step =trunc( Phi_ran*sin(th*pi/180)*Phi_ref/Theta_st )
	   endif
	else
		if (cmpstr(Holomode, "Parallel") == 0)
			The_step=trunc( Phi_ran*tan(th*deg2rad)*Phi_ref/Theta_st )
		else
			if ( cmpstr(Holomode, "h") == 0)
          			The_step=trunc( th/Theta_in*Phi_ran/Theta_st )
			else
				//altro
			endif
		endif
	endif
	
	return(The_step)
end

/// calculate delta-phi for a given theta
///
/// adapted from XPDplot 8.03
static function calc_phi_step(Theta_in, th, Theta_st, Phi_ran, Phi_ref, Holomode)
	Variable Theta_in, th, Theta_st, Phi_ran, Phi_ref
	String Holomode
	
	Variable Phi_st
	Variable deg2rad=0.01745329

	if ( cmpstr(Holomode, "Stereographic") == 0 )
		if ((th < 0.5) || (trunc(Phi_ran*sin(th*deg2rad)*Phi_ref/Theta_st) == 0))
			Phi_st=0.0
		else
			Phi_st=Phi_ran/(trunc(Phi_ran*sin(th*deg2rad)*Phi_ref/Theta_st))
		endif
		if(th==90)
			Phi_st=2.0
		endif
	endif

	if ( cmpstr(Holomode, "Parallel") == 0 )
		if((th < 0.5) || (trunc(Phi_ran*tan(th*deg2rad)*Phi_ref/Theta_st) == 0))
			Phi_st=0.0
		else
			Phi_st=Phi_ran/(trunc(Phi_ran*tan(th*deg2rad)*Phi_ref/Theta_st))
		endif
	endif

	if ( cmpstr(Holomode, "h") == 0 )
		if((th < 0.5) || (trunc(Phi_ran*sin(th*deg2rad)*Phi_ref/Theta_st) == 0))
			Phi_st=0.0
		else
		 	Phi_st=Phi_ran/trunc(th/Theta_in*Phi_ran/Theta_st)
		endif
	endif
	
	if (Phi_st==0)
		Phi_st=360
	endif
	
	return(Phi_st)
end

/// calculate delta-theta for a given theta
///
/// adapted from XPDplot 8.03
static function Calc_The_step(th, Theta_st, Holomode)
	String Holomode
	Variable th, Theta_st
	
	Variable deg2rad=0.01745329, dt_loc,The_step

	if ( (cmpstr(Holomode, "Stereographic")) ==0 )
		The_step=Theta_st 
	endif
	
	if ( (cmpstr(Holomode, "h")) ==0 )
		The_step=Theta_st 
	endif
 
 	if ( cmpstr(Holomode, "Parallel") == 0 )
		if(th < 89.5)
			dt_loc = Theta_st/cos(th*deg2rad)
	   		if(dt_loc > 10)
	   			dt_loc=10
	   		endif
	   		The_step=dt_loc
	  	else
	   		The_step=10
	  	endif
	endif
	return(The_step)
end

/// calculate the number of thetas for a pattern
///
/// adapted from XPDplot 8.03
static function CalcN_Theta(HoloMode,Theta_in,Theta_ran,Theta_st)
	String HoloMode
	Variable Theta_in,Theta_ran,Theta_st
	Variable n_theta, aux, aux1,ii

	aux = Theta_in
	aux1= Theta_in - Theta_ran
	ii = 0
	do
		aux = aux - Calc_The_step(aux, Theta_st, HoloMode)
		if(aux<=Theta_in-Theta_ran)
			aux=Theta_in-Theta_ran
		endif
		ii = ii+1
	while((aux>aux1)%&(Theta_in-aux<=Theta_ran))																//
	n_theta=ii+1
	Return(n_theta)
end

/// create a hemispherical, constant solid angle grid
///
/// all necessary waves are created in the current data folder
/// with step size 90 / (npol - 1)
///
/// adapted from XPDplot 8.03
///
/// @param	npol			number of polar angles, determines polar and azimuthal step size.
///						recommended 91 for 1-degree steps.
///
/// @param	nickname	name prefix for waves.
///						nick name must be unique in the current data folder.
///						otherwise existing waves get overwritten.
///						may be empty.
///
/// @param	xpdplot		XPDplot compatibility
///	@arg	0			(default) create the data structures required by this module
///	@arg	1			create additional waves and notebook required by XPDplot
///
function make_hemi_grid(npol, nickname, [xpdplot])
	variable npol
	string nickname
	variable xpdplot
	
	if (ParamIsDefault(xpdplot))
		xpdplot = 0
	endif
	
	string HoloMode = "h"
	variable Theta_in = 90
	variable Theta_ran = 90
	variable Theta_st = 90 / (npol - 1)
	variable Phi_ran = 360
	variable Phi_ref = 1
	variable Phi_in = 0
	
	variable n_theta = CalcN_Theta(HoloMode, Theta_in, Theta_ran, Theta_st)

	// wave names
	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"         // Intensity wave (counts/sec)
	else
		s_prefix = ""
		s_int = "values"	// "i" is not a valid wave name
	endif
	string s_polar = s_prefix + "pol"       // thetas for each int-point of the holo
	string s_azim = s_prefix + "az"        // phis for each int-point of the holo

	string s_index = s_prefix + "index"       // starting index for each theta
	string s_theta = s_prefix + "th"        // theta values
	string s_dphi = s_prefix + "dphi"        // delta phis at each theta
	string s_nphis = s_prefix + "nphis"     // number of phis at each theta

	string s_HoloData = s_prefix + "data"          // All holo exp.- parameter information
	string s_HoloInfo = s_prefix + "info"
	
	// the following two waves are used by the pearl-anglescan procedures but not by XPDplot
	string s_tot = s_prefix + "tot"	// accumulated counts at each point
	string s_weight = s_prefix + "wt"	// total accumulation time at each point (arb. units)

	make /O/D/n=(n_theta) $s_index /wave=index
	make /O/D/n=(n_theta) $s_theta /wave=theta
	make /O/D/n=(n_theta) $s_dphi /wave=dphi
	make /O/D/n=(n_theta) $s_nphis /wave=nphis
	
	//---------- calculate phi-step-size for this theta:
	variable aux = Calc_The_step(Theta_in, Theta_st, HoloMode)
	dphi[0] =  calc_phi_step(Theta_in, Theta_in, aux, Phi_ran, Phi_ref, HoloMode)
	Theta[0] = Theta_in
	nphis[0] = calc_nth(aux, Theta_in, Theta_in, Phi_ran, Phi_ref, HoloMode) 
	Index[0] = nphis[0]
	
	//---------- calculate number of phis, phi-step, and starting-index for each theta:
	variable ii = 1
	do
		Theta[ii] = Theta[ii-1] - aux
		if(Theta[ii] <= Theta_in-Theta_ran)
			Theta[ii] = Theta_in-Theta_ran
		endif
		aux =  Calc_The_step(Theta[ii], Theta_st, HoloMode)
		dphi[ii] = calc_phi_step(Theta_in, Theta[ii],  aux, Phi_ran, Phi_ref, HoloMode)
		nphis[ii] = calc_nth(aux, Theta_in, Theta[ii], Phi_ran, Phi_ref, HoloMode)
		Index[ii] = Index[ii-1] + nphis[ii]
		ii=ii+1
	while(ii < n_theta)

	if (Index[n_theta-1]==Index[n_theta-2])
		Index[n_theta-1]=Index[n_theta-2]+1
		nphis[n_theta-1]=1
	endif

	variable NumPoints = sum(nphis, 0, numpnts(nphis))
	
	//---------- calculate theta and phi for each data point:
	make /O/D/N=(NumPoints) $s_polar /wave=polar, $s_azim /wave=azim
	note azim, "version=1.6"
	
	ii = 0
	variable StartIndex = 0
	variable EndIndex
	do
		EndIndex=Index[ii]
		Polar[StartIndex, EndIndex-1]=Theta[ii]
		Azim[StartIndex, EndIndex-1]= mod(Phi_ran+(x-StartIndex)*dphi[ii]+Phi_in,Phi_ran)
		ii = ii + 1	
		StartIndex = EndIndex	
	while(ii < n_theta)

	duplicate /o azim, $s_int /wave=values
	duplicate /o azim, $s_tot /wave=totals
	duplicate /o azim, $s_weight /wave=weights
	values = nan
	totals = 0
	weights = 0

	// XPDplot metadata
	if (xpdplot)
		string s_FileName = ""
		string s_Comment = "created by pearl-anglescan-process.ipf"
		string s_HoloMode = "Stereographic"
		variable /g gb_SpectraFile = 0
	
		Make/O/D/n=22 $s_HoloData /wave=HoloData
		HoloData[0] = NaN     // v_StartKE
		HoloData[1] = NaN     // v_StoppKE
		HoloData[6] = NumPoints
		HoloData[7] = Theta_in
		HoloData[8] = Theta_ran
		HoloData[9] = Theta_st
		HoloData[11] = Phi_in
		HoloData[12] = Phi_ran
		HoloData[13] = Theta_st
		HoloData[15] = Phi_ref
		HoloData[16] = Phi_ran
		HoloData[17] = 0 // v_HoloBit (stereographic)

		Make/O/T/n=22 $s_HoloInfo /wave=HoloInfo
		HoloInfo[0] = s_FileName
		HoloInfo[1] = s_Comment
		HoloInfo[10] = s_HoloMode
		HoloInfo[11] = "" // s_MeasuringMode

		// notebook for XPDplot
		if (WinType(NickName) == 5)
			Notebook $NickName selection={startOfFile, endOfFile}	
			Notebook $NickName text=""
		else
			NewNotebook /F=0 /K=1 /N=$NickName /W=(5,40,341,260)
			Notebook $NickName defaultTab=140
			Notebook $NickName statusWidth=300
			Notebook $NickName backRGB=(56797,56797,56797)
			Notebook $NickName pageMargins={80,80,80,80}
			Notebook $NickName fSize=10
			Notebook $NickName fStyle=0,textRGB=(65535,0,26214) 
			Notebook $NickName textRGB=(65535,0,26214) 
		endif
		Notebook $NickName text = "File:\t" + s_FileName + "\r"
		Notebook $NickName text = "*** " + s_Comment + " ***\r\r"
		Notebook $NickName text = "Angle-Mode:\t" + s_HoloMode + "\r"
		Notebook $NickName text = "XPDplot Nickname:\t" + NickName + "\r"
	endif
end

/// finds the nick name given any hemi wave
///
/// the nick name is either the name of a child folder in the current data folder (PEARL specification),
/// or a prefix of the hemi wave names (XPDplot specification).
///
/// @returns the nick name
function /s get_hemi_nickname(w)
	wave w
	
	string prefix = get_hemi_prefix(w)
	string wname = nameofwave(w)
	string nickname
	
	if (strlen(prefix))
		nickname = prefix
	else
		string s_wave_df = GetWavesDataFolder(w, 1)
		dfref	parent_df = $(s_wave_df + "::")
		nickname = GetDataFolder(0, parent_df)
	endif
	
	return nickname
end

/// finds the prefix given any hemi wave
///
/// the prefix is the part of the wave name before the first underscore.
/// the prefix is used by XPDplot where it is identical to the nick name.
/// the prefix is empty in the PEARL specification.
///
/// @returns the prefix
function /s get_hemi_prefix(w)
	wave w

	string wname = nameofwave(w)
	string prefix
	if (ItemsInList(wname, "_") >= 2)
		prefix = StringFromList(0, wname, "_")
	else
		prefix = ""
	endif
	
	return prefix
end

/// finds the folder, prefix and name of holo waves given their nick name
///
/// the function looks for holo waves in the following order:
/// 1. if nickname is empty, check for prefix-less waves in current folder.
/// 2. if nickname is the name of a child folder in the current data folder, 
/// clear the (prefix-less) waves in the child folder.
/// 3. nickname is prefix of waves in current folder.
/// 4. nickname is prefix of waves in root folder.
///
/// @param[in]	nickname	folder name or name prefix of holo waves. may be empty.
///
/// @param[out]	prefix		name prefix of waves. may be empty.
///
/// @param[out]	intwave		name of intensity/values wave
///
/// @returns					reference of the data folder which contains the waves
///
function /df find_hemi_data(nickname, prefix, intwave)
	string nickname
	string &prefix
	string &intwave

	dfref datadf
	prefix = ""
	intwave = "values"
	if (strlen(nickname))
		if (DataFolderExists(nickname))
			datadf = $nickname
		else
			datadf = getdatafolderdfr()
			prefix = nickname + "_"
			intwave = prefix + "i"
			if (exists(intwave) != 1)
				datadf = root:
			endif
		endif
	else
		datadf = getdatafolderdfr()
		prefix = ""
		intwave = "values"
	endif
	return datadf
end

/// clear a hemispherical scan grid
///
/// values and weights waves are set to zero.
/// the intensity wave is set to NaN.
///
/// @param	nickname	folder name or name prefix of holo waves. may be empty.
///
function clear_hemi_grid(nickname)
	string nickname

	dfref datadf
	string s_prefix
	string s_int
	datadf = find_hemi_data(nickname, s_prefix, s_int)

	string s_totals = s_prefix + "tot"
	string s_weights = s_prefix + "wt"
	
	wave /sdfr=datadf /z w_values = $s_int
	wave /sdfr=datadf /z w_totals = $s_totals
	wave /sdfr=datadf /z w_weights = $s_weights

	if (waveexists(w_totals))
		w_totals = 0
	endif
	if (waveexists(w_weights))
		w_weights = 0
	endif
	if (waveexists(w_values))
		w_values = nan	
	endif
end

/// duplicate a hemispherical scan dataset.
///
/// this function works only for hemi scans created by make_hemi_grid() (or compatible functions).
/// the angle grid is recreated rather than copied point-by-point.
/// the new dataset is independent from the original one.
///
/// if the version of the source dataset is pre 1.6, it is converted to version 1.6.
///
/// @param	source_nickname		name prefix for waves. source data must be in current data folder.
///
/// @param	dest_folder			destination folder. folder must exist.
///
/// @param	dest_nickname		name prefix for destination waves.
///								must be unique in the current data folder.
///								otherwise existing waves get overwritten.
///								may be empty.
///
/// @param	xpdplot		XPDplot compatibility
///	@arg	0			(default) create the data structures required by this module
///	@arg	1			create additional waves and notebook required by XPDplot
///
function duplicate_hemi_scan(source_nickname, dest_folder, dest_nickname, [xpdplot])
	string source_nickname
	dfref dest_folder
	string dest_nickname
	variable xpdplot
	
	if (ParamIsDefault(xpdplot))
		xpdplot = 0
	endif
	
	dfref savedf = getdatafolderdfr()

	// source data
	string s_prefix = ""
	string s_int = "values"
	dfref source_df = find_hemi_data(source_nickname, s_prefix, s_int)
	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	string s_theta = s_prefix + "th"
	string s_tot = s_prefix + "tot"
	string s_weight = s_prefix + "wt"
	string s_matrix = s_prefix + "matrix"

	wave /sdfr=source_df theta1 = $s_theta
	wave /sdfr=source_df polar1 = $s_polar
	wave /sdfr=source_df azim1 = $s_azim
	wave /sdfr=source_df tot1 = $s_tot
	wave /sdfr=source_df weight1 = $s_weight
	wave /sdfr=source_df values1 = $s_int
	wave /sdfr=source_df /z matrix1 = $s_matrix
	
	variable npol = numpnts(theta1)
	
	setdatafolder dest_folder
	make_hemi_grid(npol, dest_nickname, xpdplot=xpdplot)
	
	// dest data
	dfref dest_df = find_hemi_data(dest_nickname, s_prefix, s_int)
	s_polar = s_prefix + "pol"
	s_azim = s_prefix + "az"
	s_theta = s_prefix + "th"
	s_tot = s_prefix + "tot"
	s_weight = s_prefix + "wt"
	s_matrix = s_prefix + "matrix"

	wave /sdfr=dest_df theta2 = $s_theta
	wave /sdfr=dest_df polar2 = $s_polar
	wave /sdfr=dest_df azim2 = $s_azim
	wave /sdfr=dest_df tot2 = $s_tot
	wave /sdfr=dest_df weight2 = $s_weight
	wave /sdfr=dest_df values2 = $s_int

	tot2 = tot1
	weight2 = weight1
	values2 = values1
	if (waveexists(matrix1))
		setdatafolder dest_df
		duplicate /o matrix1, $s_matrix
	endif
	
	if (!(NumberByKey("version", note(azim1), "=", "\r") >= 1.6))
		azim2 += 180	// changed 151030 (v1.6)
		azim2 = azim2 >= 360 ? azim2 - 360 : azim2
	endif
	
	setdatafolder saveDF
end

/// azimuthally rotate a hemispherical scan dataset.
///
/// this function works only for hemi scans created by make_hemi_grid() (or compatible functions).
///
/// @param	nickname		name prefix for waves. source data must be in current data folder.
/// @param	angle			azimuthal rotation angle in degrees.
///
function rotate_hemi_scan(nickname, angle)
	string nickname
	variable angle
	
	dfref savedf = getdatafolderdfr()

	string s_prefix = ""
	string s_int = "values"
	dfref df = find_hemi_data(nickname, s_prefix, s_int)

	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	string s_tot = s_prefix + "tot"
	string s_weight = s_prefix + "wt"
	
	wave /sdfr=df polar = $s_polar
	wave /sdfr=df azim = $s_azim
	wave /sdfr=df tot = $s_tot
	wave /sdfr=df weight = $s_weight
	wave /sdfr=df values = $s_int

	azim += angle
	azim = azim < 0 ? azim + 360 : azim
	azim = azim >= 360 ? azim - 360 : azim

	duplicate /free polar, neg_polar
	neg_polar = -polar
	sort {neg_polar, azim}, polar, azim, tot, weight, values
	
	setdatafolder saveDF
end

/// create waves for plotting a hemispherical angle scan.
///
/// the scan data must exist in the current data folder.
///
/// @param	nickname		name prefix of holo waves.
///							may be empty.
///
/// @param	projection		mapping function from polar to cartesian coordinates.
///							see @ref PageProjections for details.
///	@arg	kProjDist = 0		azimuthal equidistant
///	@arg	kProjStereo = 1	stereographic (default)
///	@arg	kProjArea = 2	azimuthal equal-area
///	@arg	kProjGnom = 3	gnomonic (0 <= polar < 90)
///	@arg	kProjOrtho = 4	orthographic
///
///
function /s prepare_hemi_scan_display(nickname, [projection])
	string nickname
	variable projection
	
	dfref savedf = getdatafolderdfr()

	if (ParamIsDefault(projection))
		projection = 1
	endif
	
	// hemi grid waves
	string s_prefix = ""
	string s_int = "values"
	dfref df = find_hemi_data(nickname, s_prefix, s_int)

	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	
	wave /sdfr=df /z values = $s_int
	wave /sdfr=df /z azim = $s_azim
	wave /sdfr=df /z polar = $s_polar

	setdatafolder df
	string s_ster_rad = s_prefix + "ster_rad"
	duplicate /o polar, $s_ster_rad /wave=ster_rad
	ster_rad = calc_graph_radius(polar, projection=projection)
	
	string s_ster_x = s_prefix + "ster_x"
	string s_ster_y = s_prefix + "ster_y"
	duplicate /o azim, $s_ster_x /wave=ster_x, $s_ster_y /wave=ster_y
	ster_x = ster_rad * cos(azim * pi / 180)
	ster_y = ster_rad * sin(azim * pi / 180)

	setdatafolder savedf
end

/// display a plot of a hemispherical angle scan.
///
/// the scan data must exist in the current data folder.
/// azimuth = 0 should be at 9 o'clock.
/// then the orientation is the same as the sample at normal emission and phi = 0,
/// the handle of the sample plate pointing to the left.
///
/// @param	nickname		name prefix of holo waves.
///							may be empty.
///
/// @param	projection		mapping function from polar to cartesian coordinates.
///							see @ref PageProjections for details.
///	@arg	kProjDist = 0		azimuthal equidistant
///	@arg	kProjStereo = 1	stereographic (default)
///	@arg	kProjArea = 2	azimuthal equal-area
///	@arg	kProjGnom = 3	gnomonic (0 <= polar < 90)
///	@arg	kProjOrtho = 4	orthographic
///
/// @param	graphtype		type of graph
///	@arg	1				(pol, az) trace in Igor "New Polar" (default).
///	@arg	2				XPDplot (reserved, not implemented).
///	@arg	3				matrix in Igor "New Polar".
///							the matrix wave is a 2D wave with X and Y scaling corresponding to the selected projection.
///							matrix waves can be created by interpolate_hemi_scan().
///							note: the pol and az waves are required as well.
///
/// @param	do_ticks			select which ticks to draw. 
///							value must be the arithmetic OR of all selected items.
///							default: 3
///	@arg	0				none
///	@arg	1				major azimuthal
///	@arg	2				minor azimuthal
///S
/// @param	do_grids			select which grids to draw. 
///							value must be the arithmetic OR of all selected items.
///							default: 3
///	@arg	0				none
///	@arg	1				radius at 0 and 90 degree azimuth
///	@arg	2				circle at 30 and 60 degree polar
///
/// @param	graphname		name of graph window. default: nickname
///							if empty, a default name is assigned.
///							if a window with this name is existing, the function brings it to the front, and does nothing else.
///
/// @returns the name of the graph window
///
function /s display_hemi_scan(nickname, [projection, graphtype, do_ticks, do_grids, graphname])
	string nickname
	variable projection
	variable graphtype
	variable do_ticks
	variable do_grids
	string graphname
	
	dfref savedf = getdatafolderdfr()

	if (ParamIsDefault(projection))
		projection = 1
	endif
	if (ParamIsDefault(graphtype))
		graphtype = 1
	endif
	if (ParamIsDefault(do_ticks))
		do_ticks = 3
	endif
	if (ParamIsDefault(do_grids))
		do_grids = 3
	endif
	if (ParamIsDefault(graphname))
		if (strlen(nickname) > 0)
			graphname = nickname
		else
			graphname = GetDataFolder(0)
		endif
	endif
	
	prepare_hemi_scan_display(nickname, projection=projection)
	
	// hemi grid waves
	string s_prefix = ""
	string s_int = "values"
	dfref df = find_hemi_data(nickname, s_prefix, s_int)

	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	string s_matrix = s_prefix + "matrix"
	string s_ster_rad = s_prefix + "ster_rad"
	
	wave /sdfr=df /z values = $s_int
	wave /sdfr=df /z azim = $s_azim
	wave /sdfr=df /z polar = $s_polar
	wave /sdfr=df /z ster_rad = $s_ster_rad
	wave /sdfr=df /z matrix = $s_matrix

	setdatafolder df
	variable azim_offset = 0
	if (!(NumberByKey("version", note(azim), "=", "\r") >= 1.6))
		DoAlert /T="display hemi scan" 0, "your dataset doesn't include the version 1.6 flag. if it was created with an earlier version that might be okay. please check that the orientation is correct!"
		azim_offset = 180	// changed 151030 (v1.6)
	endif

	string s_trace
	DoWindow $graphname
	if (v_flag)
		if (str2num(GetUserData(graphname, "", "graphtype")) == graphtype)
			// graph exists and will update automatically - do not recreate
			graphtype = 0
		else
			// graph exists - but needs recreating
			killwindow $graphname
		endif
	endif
	
	switch(graphtype)
		case 1:
			graphname = display_polar_graph(graphname, angle_offset=azim_offset, do_ticks=do_ticks)

			s_trace = WMPolarAppendTrace(graphname, ster_rad, azim, 360)
			ModifyGraph /W=$graphname mode($s_trace)=2, lsize($s_trace)=2
			ModifyGraph /W=$graphname zColor($s_trace)={values,*,*,BlueGreenOrange,0}
			
			ColorScale /W=$graphname /C /N=text0 /E=2 /F=0 /B=1 /A=RB /X=0.00 /Y=0.00 trace=polarY0
			ColorScale /W=$graphname /C /N=text0 side=2, width=5, heightPct=40, frame=0.50, lblMargin=0
			ColorScale /W=$graphname /C /N=text0 nticks=2, minor=1, tickLen=4.00, tickThick=0.50

			SetWindow $graphname, userdata(projection)=num2str(projection)
			SetWindow $graphname, userdata(graphtype)=num2str(graphtype)
			draw_hemi_axes(graphname, do_grids=do_grids)
			break
		case 3:
			graphname = display_polar_graph(graphname, angle_offset=azim_offset, do_ticks=do_ticks)

			s_trace = WMPolarAppendTrace(graphname, ster_rad, azim, 360)
			ModifyGraph /W=$graphname mode($s_trace)=0, lsize($s_trace)=0
			AppendImage /L=VertCrossing /B=HorizCrossing matrix

			ColorScale /W=$graphname /C /N=text0 /E=2 /F=0 /B=1 /A=RB /X=0.00 /Y=0.00 image=$s_matrix
			ColorScale /W=$graphname /C /N=text0 side=2, width=5, heightPct=40, frame=0.50, lblMargin=0
			ColorScale /W=$graphname /C /N=text0 nticks=2, minor=1, tickLen=4.00, tickThick=0.50
		
			SetWindow $graphname, userdata(projection)=num2str(projection)
			SetWindow $graphname, userdata(graphtype)=num2str(graphtype)
			draw_hemi_axes(graphname, do_grids=do_grids)
			break
	endswitch
	
	setdatafolder savedf
	return graphname
end

/// displays an empty polar graph
///
/// the graph is drawn using Wavemetrics "New Polar Graphs.ipf".
///
/// initially the graph is empty.
/// hemispherical scans are displayed by adding a trace that coveres the whole plot area,
/// and setting the trace color to a function of the intensity.
/// traces are added by calling WMPolarAppendTrace.
///
/// the following items of the graph > packages menu might be useful:
/// * modify polar graph
/// * color table control
/// * show polar cursors
/// * polar graph legend
///	
/// parameters can be changed programmatically as shown in the code of this function.
/// after programmatic parameter changes, call WMPolarAxesRedrawGraphNow(graphname).
///
/// @param	graphname	requested name of new graph window.
///						if empty, a default name is assigned.
///						if a window with this name is existing, the function brings it to the front, and does nothing else.
///
/// @param	angle_offset	azimuth (on screen) where angle 0 is plotted
///						(zeroAngleWhere parameter of polar graphs).
///						starting with version 1.6, the default is 0.
///						for hemi grids created with earlier versions, 
///						it should be set to 180 for correct orientation.
///
/// @param	do_ticks		select which ticks to draw. 
///						value must be the arithmetic OR of all selected items.
///						default: 3
///	@arg	0			none
///	@arg	1			major azimuthal
///	@arg	2			minor azimuthal
///
/// @returns the name of the graph window.
///
/// @version 1.7
/// interface change: the trace drawing code is moved to display_hemi_scan,
/// so that this function can be reused by other graph types, e.g. display_scanlines.
///
static function /s display_polar_graph(graphname, [angle_offset, do_ticks])
	
	string graphname
	variable angle_offset
	variable do_ticks
	
	dfref savedf = GetDataFolderDFR()
	
	if (ParamIsDefault(angle_offset))
		angle_offset = 0
	endif
	if (ParamIsDefault(do_ticks))
		do_ticks = 3
	endif
	
	if ((strlen(graphname) == 0) || (wintype(graphname) == 0))
		Display /k=1 /W=(10,45,360,345)
		DoWindow /C $graphname
		graphname = WMNewPolarGraph("", graphname)
		WMPolarGraphSetVar(graphname, "zeroAngleWhere", angle_offset)
		
		WMPolarGraphSetVar(graphname, "angleAxisThick", 0.5)
		WMPolarGraphSetStr(graphname, "doMajorAngleTicks", "manual")
		WMPolarGraphSetVar(graphname, "majorAngleInc", 30) // major ticks in 30 deg steps
		WMPolarGraphSetVar(graphname, "minorAngleTicks", 2) // minor ticks in 10 deg steps
		WMPolarGraphSetStr(graphname, "angleTicksLocation", "Outside")
		WMPolarGraphSetVar(graphname, "doAngleTickLabelSubRange", 1)
		WMPolarGraphSetVar(graphname, "angleTickLabelRangeStart", 0)
		WMPolarGraphSetVar(graphname, "angleTickLabelRangeExtent", 90)
		WMPolarGraphSetStr(graphname, "angleTickLabelNotation", "%g°")
  
		WMPolarGraphSetVar(graphname, "doPolarGrids", 0)
		WMPolarGraphSetVar(graphname, "doRadiusTickLabels", 0)
		WMPolarGraphSetStr(graphname, "radiusAxesWhere", "  Off")	// note the leading spaces, cf. WMPolarAnglesForRadiusAxes
		WMPolarGraphSetStr(graphname, "radiusTicksLocation", "Off")

		WMPolarGraphSetVar(graphname, "majorTickLength", 2)
		WMPolarGraphSetVar(graphname, "majorTickThick", 0.5)
		WMPolarGraphSetVar(graphname, "minorTickLength", 1)
		WMPolarGraphSetVar(graphname, "minorTickThick", 0.5)
		WMPolarGraphSetVar(graphname, "tickLabelOpaque", 0)
		WMPolarGraphSetVar(graphname, "tickLabelFontSize", 7)

		// changes
		if (do_ticks & 1)
			WMPolarGraphSetStr(graphname, "angleTicksLocation", "Outside")
		else
			WMPolarGraphSetStr(graphname, "angleTicksLocation", "Off")
		endif
		if (do_ticks & 2)
			WMPolarGraphSetVar(graphname, "doMinorAngleTicks", 1)
		else
			WMPolarGraphSetVar(graphname, "doMinorAngleTicks", 0)
		endif
		
		DoWindow /T $graphname, graphname

		// cursor info in angles
		string graphdf = "root:packages:WMPolarGraphs:" + graphname
		setdatafolder graphdf
		// current theta, phi coordinates are stored in global variables in the package folder of the graph
		variable /g csrA_theta
		variable /g csrA_phi
		variable /g csrB_theta
		variable /g csrB_phi
		// the text box is hidden initially. it shows up and hides with the cursor info box.
		string tb 
		tb = "\\{"
		tb = tb + "\"A = (%.1f, %.1f)\","
		tb = tb + graphdf + ":csrA_theta,"
		tb = tb + graphdf + ":csrA_phi"
		tb = tb + "}"
		TextBox /W=$graphname /A=LT /B=1 /E=2 /F=0 /N=tb_angles /X=1 /Y=1 /V=0 tb
		tb = "\\{"
		tb = tb + "\"B = (%.1f, %.1f)\","
		tb = tb + graphdf + ":csrB_theta,"
		tb = tb + graphdf + ":csrB_phi"
		tb = tb + "}"
		AppendText /W=$graphname /N=tb_angles tb
		// updates are triggered by a window hook
		SetWindow $graphname, hook(polar_graph_hook)=PearlAnglescanProcess#polar_graph_hook
	else
		// graph window exists
		DoWindow /F $graphname
	endif

	setdatafolder savedf
	return graphname
end

/// draw polar and azimuthal grids in an existing polar graph.
///
/// the function adds the following draw objects to a polar graph:
/// * concentric circles at polar angles 0, 30, and 60 degrees with labels.
/// * radial axes at 0 and 90 degree azimuth.
///
/// the objects are added to the ProgFront drawing layer and will appear in front of the data trace.
/// in interactive drawing mode, you can select the active drawing layer by clicking the tree icon
/// while holding the Alt key.
///
/// the graph must have been created by display_polar_graph().
/// the function reads the projection mode from the window user data "projection".
///
/// @param	graphname	name of graph window.
///
/// @param	do_grids			select which optional grids to draw. 
///							value must be the arithmetic OR of all selected items.
///							default: 3
///	@arg	0				none
///	@arg	1				radius at 0 and 90 degree azimuth
///	@arg	2				circle at 30 and 60 degree polar
///
/// @warning EXPERIMENTAL!
///			this function is under development.
///			the interface and behaviour of this function may change significantly in future versions.
static function /s draw_hemi_axes(graphname, [do_grids])
	string graphname
	variable do_grids
	
	if (ParamIsDefault(do_grids))
		do_grids = 3
	endif

	dfref savedf = GetDataFolderDFR()

	string sproj = GetUserData(graphname, "", "projection")
	variable projection = str2num("0" + sproj)

	SetDrawLayer /W=$graphname ProgFront

	// polar axis
	SetDrawEnv /W=$graphname xcoord=HorizCrossing, ycoord=VertCrossing
	SetDrawEnv /W=$graphname linethick= 0.5
	SetDrawEnv /W=$graphname dash=2
	SetDrawEnv /W=$graphname fillpat=0
	SetDrawEnv /W=$graphname fname="default", fsize=7
	SetDrawEnv /W=$graphname textxjust=1, textyjust=1
	//SetDrawEnv /W=$graphname linefgc=(65535,65535,65535)
	SetDrawEnv /W=$graphname save

	if (do_grids & 1)
		DrawLine /W=$graphname 0, -2, 0, 2
		DrawLine /W=$graphname -2, 0, 2, 0
	endif

	variable radi
	if (do_grids & 2)
		radi = calc_graph_radius(0.5, projection=projection)
		DrawOval /W=$graphname -radi, radi, radi, -radi
		radi = calc_graph_radius(30, projection=projection)
		DrawOval /W=$graphname -radi, radi, radi, -radi
		radi = calc_graph_radius(60, projection=projection)
		DrawOval /W=$graphname -radi, radi, radi, -radi

		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°"
		radi = calc_graph_radius(60, projection=projection)
		DrawText /W=$graphname radi, -0.1, "60°"
	endif
	
	setdatafolder savedf
end

/// draw the circle of a diffraction cone in a stereographic polar graph.
///
/// the diffraction cone consists of a circle marking the diffraction ring, and a dot marking the axis.
/// the cone is drawn as a group of draw objects on the UserFront layer.
/// the objects can be edited interactively.
///
/// @param graphname		name of graph window (not implemented yet).
///
/// @param groupname		name of a drawing group.
///							if the group exists (from a previous cone) it is replaced.
///							if the group doesn't exist, a new one is created.
///
/// @param theta_axis			polar angle of the cone axis in degrees.
///
/// @param theta_inner		polar angle of the innermost point of the circle in degrees.
///
/// @param phi				azimuthal angle of the cone axis in degrees.
///
/// @warning EXPERIMENTAL!
///			this function is under development.
///			the interface and behaviour of this function may change significantly in future versions.
///
function draw_diffraction_cone(graphname, groupname, theta_axis, theta_inner, phi)
	string graphname
	string groupname
	
	variable theta_axis
	variable theta_inner
	variable phi
	
	variable r_axis = calc_graph_radius(theta_axis)
	variable r_inner = calc_graph_radius(theta_inner)
	variable r_outer = calc_graph_radius(2 * theta_axis - theta_inner)
	
	SetDrawEnv push
	SetDrawLayer UserFront
	DrawAction getgroup=$groupname, delete
	SetDrawEnv gstart, gname=$groupname
	variable xc, yc, xr, yr

	// cone periphery
	variable r_center = (r_outer + r_inner) / 2
	variable r_radius = (r_outer - r_inner) / 2
	xc = r_center * cos(phi * pi / 180)
	yc = r_center * sin(phi * pi / 180)
	xr = r_radius
	yr = r_radius
	SetDrawEnv xcoord=HorizCrossing, ycoord=VertCrossing
	SetDrawEnv dash=11, fillpat=0
	DrawOval xc - xr, yc - yr, xc + xr, yc + yr

	// cone axis	
	xc = r_axis * cos(phi * pi / 180)
	yc = r_axis * sin(phi * pi / 180)
	r_radius = calc_graph_radius(2)
	xr = r_radius
	yr = r_radius
	SetDrawEnv xcoord=HorizCrossing, ycoord=VertCrossing
	SetDrawEnv fillfgc=(0,0,0)
	DrawOval xc - xr, yc - yr, xc + xr, yc + yr
	
	SetDrawEnv gstop
	SetDrawEnv pop
end

/// display a polar graph with lines indicating the angles covered by an angle scan.
///
/// @param	nickname	nick name for output data.
///						this will become the name of a child folder containing the output.
///
/// @param	alpha_lo		low limit of the analyser angle.
///
/// @param	alpha_hi		high limit of the analyser angle.
///
/// @param	m_theta		manipulator theta angles, 0 = normal emission. size = dimsize(data, 1)
///
/// @param	m_tilt		manipulator tilt angles, 0 = normal emission. size = dimsize(data, 1)
///
/// @param	m_phi		manipulator phi angles, 0 = azimuthal origin. size = dimsize(data, 1)
///
/// @param	folding		rotational averaging, default = 1
///
/// @param	projection	mapping function from polar to cartesian coordinates. see calc_graph_radius().
///
/// @remark	this function is extremely slow.
///
function /s display_scanlines(nickname, alpha_lo, alpha_hi, m_theta, m_tilt, m_phi, [folding, projection])
	string nickname
	variable alpha_lo
	variable alpha_hi
	wave m_theta
	wave m_tilt
	wave m_phi
	variable folding
	variable projection
	
	if (ParamIsDefault(folding))
		folding = 1
	endif
	if (ParamIsDefault(projection))
		projection = 1
	endif
	
	// sort out data folder structure
	dfref saveDF = GetDataFolderDFR()
	newdatafolder /s/o $nickname
	string graphname = "graph_" + nickname
	
	duplicate /free m_tilt, loc_m_tilt
	loc_m_tilt = -m_tilt
	
	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
	analyser = x
	
	convert_angles_ttpa2polar(m_theta, loc_m_tilt, m_phi, analyser, d_polar, d_azi)
	duplicate /free d_polar, d_radius
	d_radius = calc_graph_radius(d_polar, projection=projection)
	d_azi += 180	// changed 151030 (v1.6)
	
	graphname = display_polar_graph(graphname)
	SetWindow $graphname, userdata(projection)=num2str(projection)
	
	variable ifold
	variable iang
	variable nang = numpnts(m_theta)
	string s_rad
	string s_azi
	string s_trace
	for (ifold = 0; ifold < folding; ifold += 1)
		d_azi = d_azi >= 360 ? d_azi - 360 : d_azi
		for (iang = 0; iang < nang; iang += 1)
			sprintf s_rad, "rad_%d_%d", ifold, iang
			duplicate /o analyser, $s_rad
			wave w_rad = $s_rad
			w_rad = d_radius[p][iang]
			
			sprintf s_azi, "azi_%d_%d", ifold, iang
			duplicate /o analyser, $s_azi
			wave w_azi = $s_azi
			w_azi = d_azi[p][iang]
			
			if (numtype(sum(w_rad)) == 0)
				s_trace = WMPolarAppendTrace(graphname, w_rad, w_azi, 360)
				ModifyGraph /w=$graphname mode($s_trace)=0, lsize($s_trace)=0.5
			endif
		endfor
		d_azi += 360 / folding
	endfor

	draw_hemi_axes(graphname)

	setdatafolder saveDF
	return graphname
end

/// @page PageProjections Projections
///
/// the functions of the anglescan package support the following map projections.
/// for a description of the different projections, see, for example,
/// https://en.wikipedia.org/wiki/Map_projection
///
/// | Selector | Projection | Function | Properties |
/// | :----: | :----: | :----: | :---- |
/// | kProjDist = 0 | azimuthal equidistant | r = c * theta | radius is proportional to polar angle. |
/// | kProjStereo = 1 | stereographic | r = c * tan theta/2 | circles on sphere map to circles. |
/// | kProjArea = 2 | azimuthal equal-area | r = c * sin theta/2 | preserves area measure. |
/// | kProjGnom = 3 | gnomonic | r = c * tan theta | great circles map to straight lines. |
/// | kProjOrtho = 4 | orthographic | r = c * sin theta | k-space mapping in ARPES and LEED. |
///
/// the projections in this package are defined for 0 <= theta < 90.
///
constant kProjDist = 0
constant kProjStereo = 1
constant kProjArea = 2
constant kProjGnom = 3
constant kProjOrtho = 4

static constant kProjScaleDist = 2
static constant kProjScaleStereo = 2
static constant kProjScaleArea = 2
// scaled so that radius(gnom) = radius(stereo) for polar = 88
static constant kProjScaleGnom = 0.06744519021
static constant kProjScaleOrtho = 2

/// calculate the projected polar angle
///
/// @param	polar			polar angle in degrees
///
/// @param	projection		mapping function from polar to cartesian coordinates.
///							see @ref PageProjections for details.
///	@arg	kProjDist = 0		azimuthal equidistant
///	@arg	kProjStereo = 1	stereographic (default)
///	@arg	kProjArea = 2	azimuthal equal-area
///	@arg	kProjGnom = 3	gnomonic (0 <= polar < 90)
///	@arg	kProjOrtho = 4	orthographic
///
/// @return projected radius.
///	the radius is scaled such that grazing emission maps to 2.
threadsafe function calc_graph_radius(polar, [projection])
	variable polar
	variable projection
	
	if (ParamIsDefault(projection))
		projection = 1
	endif

	variable radius
	switch(projection)
		case kProjStereo: // stereographic
			radius = kProjScaleStereo * tan(polar / 2 * pi / 180)
			break
		case kProjArea: // equal area
			radius = kProjScaleArea * sin(polar / 2 * pi / 180)
			break
		case kProjGnom: // gnomonic
			radius = polar < 90 ? kProjScaleGnom * tan(polar * pi / 180) : inf
			break
		case kProjOrtho: // orthographic
			radius = kProjScaleOrtho * sin(polar * pi / 180)
			break
		default: // equidistant
			radius = kProjScaleDist * polar / 90
	endswitch
	
	return radius
end

/// calculate polar angle from Cartesian coordinate
///
/// this is the reverse mapping to calc_graph_radius()
///
/// @param	x, y				projected Cartesian coordinate
///
/// @param	projection		mapping function from polar to cartesian coordinates.
///							see @ref PageProjections for details.
///	@arg	kProjDist = 0		azimuthal equidistant
///	@arg	kProjStereo = 1	stereographic (default)
///	@arg	kProjArea = 2	azimuthal equal-area
///	@arg	kProjGnom = 3	gnomonic (0 <= polar < 90)
///	@arg	kProjOrtho = 4	orthographic
///
/// @returns polar angle in degrees
///
threadsafe function calc_graph_polar(x, y, [projection])
	variable x
	variable y
	variable projection	
	
	if (ParamIsDefault(projection))
		projection = 1
	endif

	variable radius
	variable polar
	
	radius = sqrt(x^2 + y^2)
	switch(projection)
		case kProjStereo: // stereographic
			polar = 2 * atan(radius / kProjScaleStereo) * 180 / pi
			break
		case kProjArea: // equal area
			polar = 2 * asin(radius / kProjScaleArea) * 180 / pi
			break
		case kProjGnom: // gnomonic
			polar = atan(radius / kProjScaleGnom) * 180 / pi
			break
		case kProjOrtho: // orthographic
			polar = asin(radius / kProjScaleOrtho) * 180 / pi
			break
		default: // equidistant
			polar = 90 * radius / kProjScaleDist
	endswitch

	return polar
end

/// calculate azimuthal angle from Cartesian coordinate
///
/// @param	x, y				projected Cartesian coordinate
///
/// @param	projection		mapping function from polar to cartesian coordinates.
///							all supported projections are azimuthal, they have no effect on the azimuthal coordinate.
///							see @ref PageProjections for details.
///	@arg	kProjDist = 0		azimuthal equidistant
///	@arg	kProjStereo = 1	stereographic (default)
///	@arg	kProjArea = 2	azimuthal equal-area
///	@arg	kProjGnom = 3	gnomonic (0 <= polar < 90)
///	@arg	kProjOrtho = 4	orthographic
///
/// @param	zeroAngle		zeroAngleWhere parameter of polar graphs
///	@arg	0 (default)		zero is at the 3 o'clock position
///	@arg	180				zero is at the 9 o'clock position
///	@arg	other values		not tested
///		
/// @returns polar angle in degrees
///
threadsafe function calc_graph_azi(x, y, [projection,zeroAngle])
	variable x
	variable y
	variable projection	
	variable zeroAngle
	
	if (ParamIsDefault(projection))
		projection = 1
	endif
	if (ParamIsDefault(zeroAngle))
		zeroAngle = 0
	endif

	variable azi
	if (x > 0)
		azi = atan(y / x) * 180 / pi
	else
		azi = atan(y / x) * 180 / pi + 180
	endif
	
	azi += zeroAngle
	if (azi < 0)
		azi += 360
	endif
	if (azi >= 360)
		azi -= 360
	endif
	if (numtype(azi) != 0)
		azi = 0
	endif
	
	return azi
end

/// update the angles info based on cursors A and B of a given polar graph window
///	
/// the function reads the projection mode from the user data of the graph window
/// and the zeroAngleWhere variable from the associated WMPolarGraph data folder.
///
/// the calculated angles are written to the csrA_theta, csrA_phi, csrB_theta, and csrB_phi
/// global variables in the polar graph data folder.
/// the angles text box of the graph updates from to these variables dynamically.
///
/// @param	graphname		name of polar graph window
///
static function update_polar_info(graphname)
	string graphname
	
	dfref savedf = GetDataFolderDFR()

	string graphdf = "root:packages:WMPolarGraphs:" + graphname
	setdatafolder graphdf
	
	nvar csrA_theta
	nvar csrA_phi
	nvar csrB_theta
	nvar csrB_phi

	string sproj = GetUserData(graphname, "", "projection")
	variable projection = str2num("0" + sproj)
	nvar zeroAngleWhere
	
	variable x = hcsr(A, graphname)
	variable y = vcsr(A, graphname)
	csrA_theta = calc_graph_polar(x, y, projection=projection)
	csrA_phi = calc_graph_azi(x, y, projection=projection, zeroAngle=zeroAngleWhere)

	x = hcsr(B, graphname)
	y = vcsr(B, graphname)
	csrB_theta = calc_graph_polar(x, y, projection=projection)
	csrB_phi = calc_graph_azi(x, y, projection=projection, zeroAngle=zeroAngleWhere)
	
	setdatafolder savedf
end

/// polar graph window hook
///
/// this hook converts the cursor positions to polar coordinates
/// and displays them in a text box on the graph.
/// the text box is visible while the cursor info box is visible.
static function polar_graph_hook(s)
	STRUCT WMWinHookStruct &s

	Variable hookResult = 0

	switch(s.eventCode)
		case 7:	// cursor moved
			update_polar_info(s.winname)
			break
		case 20: // show info
			TextBox /W=$s.winname /N=tb_angles /C /V=1
			break
		case 21: // hide info
			TextBox /W=$s.winname /N=tb_angles /C /V=0
			break
	endswitch

	return hookResult		// 0 if nothing done, else 1
end

function set_polar_graph_cursor(nickname, cursorname, polar_angle, azim_angle, [graphname])
	string nickname
	string cursorname
	variable polar_angle
	variable azim_angle
	string graphname
	
	if (ParamIsDefault(graphname))
		if (strlen(nickname) > 0)
			graphname = nickname
		else
			graphname = GetDataFolder(0)
		endif
	endif
	
	string s_prefix = ""
	string s_int = "values"
	dfref df = find_hemi_data(nickname, s_prefix, s_int)

	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	wave /sdfr=df /z azim = $s_azim
	wave /sdfr=df /z polar = $s_polar

	FindLevel /P /Q polar, polar_angle
	if (v_flag == 0)
		variable polar_level = floor(v_levelx)
		FindLevel /P /Q /R=[polar_level] azim, azim_angle
		if (v_flag == 0)
			variable azim_level = round(v_levelx)
			string tracename = "polarY0"
			Cursor /W=$graphname /P $cursorname $traceName azim_level
		endif
	endif
end


/// add arbitrary angle scan data to a hemispherical scan grid.
///
/// the function fills the input data into bins defined by the hemi scan grid.
/// it sums up the values and weights of the data points which fall into each 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).
///
/// 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/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
	wave values
	wave polar
	wave azi
	wave weights

	dfref savedf = GetDataFolderDFR()

	if (ParamIsDefault(weights))
		duplicate /free values, weights
		weights = 1
	endif

	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"
	string s_nphis = s_prefix + "nphis"
	
	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
	duplicate /free weights, weights_copy
	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
	
	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 = 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
	
	// destination slice coordinates
	variable ipol = BinarySearch(w_theta, polar)
	if (ipol < 0)
		return -1
	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]
	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 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 function determines the bin size at the given polar angle,
/// sums up the values and weights of the data points which fall into each 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
	wave values
	variable polar
	wave azi
	wave weights

	dfref savedf = GetDataFolderDFR()

	duplicate /free values, values_copy
	duplicate /free azi, azi_copy
	if (ParamIsDefault(weights))
		duplicate /free values, weights_copy
		weights_copy = 1
	else
		duplicate /free weights, weights_copy
	endif

	// hemi grid waves
	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"
	string s_nphis = s_prefix + "nphis"
	
	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
	
	add_aziscan_core(values_copy, weights_copy, polar, azi_copy, w_theta, w_azim, w_index, w_dphi, w_totals, w_weights)

	w_values = w_weights > 0 ? w_totals / w_weights : nan

	SetDataFolder savedf
end

/// interpolate a hemispherical scan onto a rectangular grid
///
/// the scan data must exist in the current data folder
/// or in the sub-folder given by the nickname parameter.
///
/// the interpolated data is written to a new two-dimensional wave "matrix".
/// the wave has a fixed size of 181 x 181 points optimized for 1-degree polar steps.
///
/// missing values (nan) are interpolated.
/// this works well only if the missing values are reasonable sparse.
/// the function also applies a gaussian filter to smooth the image.
/// empty rings at high polar angles map are preserved.
///
/// to display the result call display_hemi_scan() with graphtype=3.
///
/// @param	nickname		name prefix of holo waves.
///							may be empty.
///
/// @param	projection		mapping function from polar to cartesian coordinates.
///							see @ref PageProjections for details.
///	@arg	kProjDist = 0		azimuthal equidistant
///	@arg	kProjStereo = 1	stereographic (default)
///	@arg	kProjArea = 2	azimuthal equal-area
///	@arg	kProjGnom = 3	gnomonic (0 <= polar < 90)
///	@arg	kProjOrtho = 4	orthographic
///
///
function interpolate_hemi_scan(nickname, [projection])
	string nickname
	variable projection

	dfref savedf = GetDataFolderDFR()

	if (ParamIsDefault(projection))
		projection = 1
	endif
	
	string s_prefix = ""
	string s_int = "values"
	dfref df = find_hemi_data(nickname, s_prefix, s_int)
	prepare_hemi_scan_display(nickname, projection=projection)
	
	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	string s_matrix = s_prefix + "matrix"
	string s_ster_rad = s_prefix + "ster_rad"
	string s_ster_x = s_prefix + "ster_x"
	string s_ster_y = s_prefix + "ster_y"
	
	wave /sdfr=df /z values = $s_int
	wave /sdfr=df /z azim = $s_azim
	wave /sdfr=df /z polar = $s_polar
	wave /sdfr=df /z ster_rad = $s_ster_rad
	wave /sdfr=df /z ster_x = $s_ster_x
	wave /sdfr=df /z ster_y = $s_ster_y

	variable min_ster_x = wavemin(ster_x)
	variable max_ster_x = wavemax(ster_x)
	variable x0 = min_ster_x
	variable xn = 181
	variable dx = (max_ster_x - min_ster_x) / (xn - 1)
	make /n=(numpnts(ster_x), 3) /free triplet
	triplet[][0] = ster_x[p]
	triplet[][1] = ster_y[p]
	triplet[][2] = values[p]

	variable size = 181
	setdatafolder df
	make /n=(size, size) /d /o $(s_prefix + "matrix") /wave=matrix
	make /n=(size, size) /free mnorm
	ImageFromXYZ /as {ster_x, ster_y, values}, matrix, mnorm
	matrix /= mnorm
	matrixfilter NanZapMedian, matrix
	matrixfilter gauss, matrix
	
	duplicate /free values, ster_finite
	ster_finite = (numtype(values) == 0) * (ster_x^2 + ster_y^2)
	variable ster_max = wavemax(ster_finite)
	matrix = (x^2 + y^2) <= ster_max ? matrix : nan
	
	setdatafolder savedf
end

/// map angle scan data onto a rectangular grid in stereographic projection
///
/// accepts angle-scan data as returned by adh5_load_reduced,
/// maps them onto a rectangular grid in stereographic projection
///
/// @param data 2D data wave, X-axis = analyser angle, Y-axis = manipulator scan (no specific ordering required)
///
/// @pre manipulator angles as attributes in attr folder next to the data wave
///
/// @warning EXPERIMENTAL
function quick_pizza_image(data, nickname, theta_offset, tilt_offset, phi_offset, [npolar, nograph, folding])
	wave data // 2D intensity wave, see requirements above
	string nickname // nick name for output data
		// in default mode, this will be the name of a child folder containing the output
		// in XPDplot mode, this will be a prefix of the generated data in the root folder
	variable theta_offset // manipulator theta angle corresponding to normal emission
	variable tilt_offset // manipulator tilt angle corresponding to normal emission
	variable phi_offset // manipulator phi angle corresponding to phi_result = 0
	variable npolar // number of polar angles, determines polar and azimuthal step size
		// default = 91 (1 degree steps)
	variable nograph // 0 (default) = display a new polar graph
		// 1 = don't display a new graph (if a graph is existing from a previous call, it will update)
	variable folding // rotational averaging, default = 1
		
	if (ParamIsDefault(npolar))
		npolar = 91
	endif
	if (ParamIsDefault(nograph))
		nograph = 0
	endif
	if (ParamIsDefault(folding))
		folding = 1
	endif
	string graphname = "graph_" + nickname
	string s_prefix = ""
	
	// sort out data folder structure
	dfref saveDF = GetDataFolderDFR()
	dfref dataDF = GetWavesDataFolderDFR(data)
	setdatafolder dataDF
	if (DataFolderExists(":attr"))
		setdatafolder :attr
	endif
	dfref attrDF = GetDataFolderDFR()
	setdatafolder dataDF
	newdatafolder /s/o $nickname
	dfref destDF = GetDataFolderDFR()
	
	// performance monitoring
	variable timerRefNum
	variable /g xyz_perf_secs
	timerRefNum = startMSTimer

	wave /sdfr=attrDF ManipulatorTheta
	wave /sdfr=attrDF ManipulatorTilt
	wave /sdfr=attrDF ManipulatorPhi
	duplicate /free ManipulatorTheta, m_theta
	duplicate /free ManipulatorTilt, m_tilt
	duplicate /free ManipulatorPhi, m_phi
	m_theta -= theta_offset
	m_tilt -= tilt_offset
	m_tilt *= -1 // checked 140702
	m_phi -= phi_offset
	//m_phi *= -1 // checked 140702
	
	make /n=1/d/free d_polar, d_azi
	convert_angles_ttpd2polar(m_theta, m_tilt, m_phi, data, d_polar, d_azi)
	d_azi += 180	// changed 151030 (v1.6)
	d_azi = d_azi >= 360 ? d_azi - 360 : d_azi

	duplicate /free data, values
	variable nn = dimsize(values, 0) * max(dimsize(values, 1), 1)
	redimension /n=(nn) values, d_polar, d_azi
	duplicate /o d_polar, ster_rad, ster_x, ster_y

	variable projection = 1
	switch(projection)
		case 1: // stereographic
			ster_rad =  2 * tan(d_polar / 2 * pi / 180)
			break
		case 2: // azimuthal
			ster_rad = 2 * cos((180 - d_polar) / 2 * pi / 180)
			break
	endswitch
	string s_ster_x = s_prefix + "ster_x"
	string s_ster_y = s_prefix + "ster_y"
	
	nn = 401
	make /n=(nn, nn) /d /o matrix
	make /n=(nn, nn) /free mnorm
	setscale /i x -2, +2, matrix, mnorm
	setscale /i y -2, +2, matrix, mnorm
	matrix = 0
	mnorm = 0
	
	variable ifold
	for (ifold = 0; ifold < folding; ifold += 1)
		ster_x = ster_rad * cos(d_azi * pi / 180)
		ster_y = ster_rad * sin(d_azi * pi / 180)
		ImageFromXYZ {ster_x, ster_y, values}, matrix, mnorm
		d_azi = d_azi >= 180 ? d_azi + 360 / folding - 180 : d_azi + 360 / folding
	endfor
	
	matrix /= mnorm
	matrixfilter /n=5 NanZapMedian matrix
	matrixfilter /n=3 gauss matrix

	if (!nograph)
		display /k=1
		appendimage matrix
		modifygraph width={Plan,1,bottom,left}
	endif
	
	if (timerRefNum >= 0)
		xyz_perf_secs = stopMSTimer(timerRefNum) / 1e6
	endif

	setdatafolder saveDF
end

/// save a hemispherical scan to an Igor text file
function save_hemi_scan(nickname, pathname, filename)
	string nickname
	string pathname
	string filename
	
	dfref savedf = getdatafolderdfr()

	// source data
	string s_prefix = ""
	string s_int = "values"
	dfref df = find_hemi_data(nickname, s_prefix, s_int)

	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	string s_theta = s_prefix + "th"
	string s_tot = s_prefix + "tot"
	string s_weight = s_prefix + "wt"

	wave /sdfr=df theta1 = $s_theta
	wave /sdfr=df polar1 = $s_polar
	wave /sdfr=df azim1 = $s_azim
	wave /sdfr=df tot1 = $s_tot
	wave /sdfr=df weight1 = $s_weight
	wave /sdfr=df values1 = $s_int
	
	save /m="\r\n" /o /p=$pathname /t theta1, polar1, azim1, tot1, weight1, values1 as filename
	
	setdatafolder saveDF
end

/// load a hemispherical scan from an Igor text file
///
/// @todo function not implemented
function load_hemi_scan(nickname, pathname, filename)
	string nickname
	string pathname
	string filename
	
	dfref savedf = getdatafolderdfr()

	//loadwave /p=$pathname /t theta1, polar1, azim1, tot1, weight1, values1 as filename
	//LoadWave /t/p=pearl_explorer_filepath/q filename
	//svar waves = s_wavenames
	//if (v_flag > 0)
	//	string /g pearl_explorer_import = "load_itx_file"
	//endif
	
	setdatafolder saveDF
end

/// import a hemispherical scan from theta-phi-intensity waves and display it
///
/// in the tpi format, the hemi scan data is represented
/// by a triple of flat one-dimensional waves
/// corresponding to the polar angle (theta), azimuthal angle (phi) and intensity.
/// no specific sort order is required.
///
/// @param	nickname	nick name for output data
///	@arg in default mode, this will become the name of a child folder containing the output.
///	@arg in XPDplot mode, this will become a prefix of the generated data in the root folder.
///
/// @param	theta		theta angles, 0 = normal emission.
///
/// @param	phi			phi angles, 0 = azimuthal origin. size = dimsize(data, 1)
///
/// @param	intensity		intensity wave, see requirements above.
///
/// @param	npolar		number of polar angles, determines polar and azimuthal step size.
///						default = 91 (1 degree steps)
///
/// @param	folding		rotational averaging.
///						example: 3 = average to 3-fold symmetry.
///						default = 1.
///
/// @param	nograph		display a new graph window?
///	@arg	0			(default) display a new polar graph
///	@arg	1			don't display a new graph
///
/// @param	xpdplot		XPDplot compatibility
///	@arg	0			(default) create waves in child folder $nickname
///	@arg	1			create waves in root folder (compatible with XPDplot)
///
function import_tpi_scan(nickname, theta, phi, intensity, [folding, npolar, nograph, xpdplot])
	string nickname
	wave theta
	wave phi
	wave intensity
	
	variable folding
	variable npolar
	variable nograph
	variable xpdplot
	
	if (ParamIsDefault(npolar))
		npolar = 91
	endif
	if (ParamIsDefault(nograph))
		nograph = 0
	endif
	if (ParamIsDefault(folding))
		folding = 1
	endif
	if (ParamIsDefault(xpdplot))
		xpdplot = 0
	endif
	
	make_hemi_grid(npolar, nickname, xpdplot=xpdplot)

	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 += 360 / folding
	endfor

	if (nograph==0)
		display_hemi_scan(nickname)
	endif
end

/// trim a hemispherical scan at grazing angle
///
/// the function recalaculates the values wave from totals and weights
/// but sets elements above a given polar angle to nan.
///
/// @param nickname		name of the scan dataset.
///						can be empty if no prefix is used.
///						the dataset must be in the current datafolder.
///
/// @param theta_max	highest polar angle to keep (0...90 degrees).
///
function trim_hemi_scan(nickname, theta_max)
	string nickname
	variable theta_max

	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"

	wave /sdfr=df w_polar = $s_polar
	wave /sdfr=df w_values = $s_int
	wave /sdfr=df w_totals = $s_totals
	wave /sdfr=df w_weights = $s_weights

	w_values = w_polar <= theta_max ? w_totals / w_weights : nan
end

/// extract a polar cut from a hemispherical scan.
///
/// for each polar angle, the function first extracts all azimuthal angles.
/// the intensity is then interpolated between the nearest neighbours of the given azimuth.
///
/// the hemi grid must have been created in the current data folder by the make_hemi_grid function.
/// correct ordering is required.
///
/// @param nickname		name of the scan dataset.
///						can be empty if no prefix is used.
///						the dataset must be in the current datafolder.
///
/// @param azim			azimuthal angle in degrees
///
/// @return				reference of the created wave.
///						the wave has the same name as the intensity wave of the dataset
///						with the suffix "_azi" and the azimuthal angle rounded to integer.
///						it is created in the same datafolder as the original data.
///
function /wave hemi_polar_cut(nickname, azim)
	string nickname
	variable azim

	dfref savedf = getdatafolderdfr()
	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"
	string s_nphis = s_prefix + "nphis"
	string s_cut
	sprintf s_cut, "%s_azi%03u", s_int, round(azim)

	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

	variable npol = numpnts(w_theta)
	variable ipol
	variable pol_st = abs(w_theta[1] - w_theta[0])
	variable pol
	variable pol1, pol2
	variable nsel
	
	setdatafolder df
	make /n=(npol) /o $s_cut
	wave w_cut = $s_cut
	setscale /i x w_theta[0], w_theta[numpnts(w_theta)-1], "deg", w_cut
	make /n=1 /free azi_slice
	make /n=1 /free values_slice
	
	for (ipol = 0; ipol < npol; ipol += 1)
		pol = w_theta[ipol]
		pol1 = pol - pol_st / 2
		pol2 = pol + pol_st / 2
		extract /free /indx w_polar, sel, (pol1 < w_polar) && (w_polar <= pol2)
		nsel = numpnts(sel)
		if (nsel > 0)
			redimension /n=(nsel+2) azi_slice, values_slice
			azi_slice[1, nsel] = w_azim[sel[p-1]]
			azi_slice[0] = azi_slice[nsel] - 360
			azi_slice[nsel+1] = azi_slice[1] + 360
			values_slice[1, nsel] = w_values[sel[p-1]]
			values_slice[0] = values_slice[nsel]
			values_slice[nsel+1] = values_slice[1]
			w_cut[ipol] = interp(azim, azi_slice, values_slice)
		else
			w_cut[ipol] = nan
		endif
	endfor
	
	setdatafolder savedf
	return w_cut
end

/// extract an azimuthal cut from a hemispherical scan
///
/// the function extracts all azimuthal angles that are present for the given polar angle.
///
/// the hemi grid must have been created in the current data folder by the make_hemi_grid function.
/// correct ordering is required.
///
/// @param nickname		name of the scan dataset.
///						can be empty if no prefix is used.
///						the dataset must be in the current datafolder.
///
/// @param pol			polar angle in degrees
///
/// @return				reference of the created wave.
///						the wave has the same name as the intensity wave of the dataset
///						with the suffix "_azi" and the azimuthal angle rounded to integer.
///						it is created in the same datafolder as the original data.
///
function /wave hemi_azi_cut(nickname, pol)
	string nickname
	variable pol

	dfref savedf = getdatafolderdfr()
	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"
	string s_nphis = s_prefix + "nphis"
	string s_cut
	sprintf s_cut, "%s_pol%03u", s_int, round(pol)

	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

	variable pol_st = abs(w_theta[1] - w_theta[0])
	variable pol1, pol2
	variable nsel

	pol1 = pol - pol_st / 2
	pol2 = pol + pol_st / 2
	extract /free /indx w_polar, sel, (pol1 < w_polar) && (w_polar <= pol2)
	nsel = numpnts(sel)
	if (nsel > 0)
		setdatafolder df
		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
		setdatafolder savedf
		return w_cut
	else
		setdatafolder savedf
		return $""
	endif
	setdatafolder savedf
end

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 /z index = w_quantilesindex
	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.
///
/// set the minimum and maximum values of the pseudocolor scale
/// such that a specified percentile of the distribution lies outside the limits.
///
/// the new contrast is applied to traces and images of the selected graph
/// that have pseudocolor tables.
///
/// the function is not specific to angle scans.
/// it can be used for any pseudocolor trace or image plots except contour plots.
///
/// @param	pcmin		percentile below the minimum color (0-100).
/// @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, reversecolors, symmetric])
	variable pcmin
	variable pcmax
	string graphname
	string colortable
	variable reversecolors
	variable symmetric
	
	if (ParamIsDefault(graphname))
		graphname = ""
	endif
	if (ParamIsDefault(colortable))
		colortable = ""
	endif
	if (ParamIsDefault(reversecolors))
		reversecolors = 0
	endif
	if (ParamIsDefault(symmetric))
		symmetric = 0
	endif
	
	dfref save_df = GetDataFolderDFR()
	
	string objname
	string info
	string wname
	string ctab
	variable rev
	variable n
	variable i
	variable vmin
	variable vmax

	string traces = TraceNameList(graphname, ";", 1+4)
	n = ItemsInList(traces, ";")
	for (i = 0; i < n; i += 1)
		objname = StringFromList(i, traces, ";")
		info = TraceInfo(graphname, objname, 0)
		if (strlen(info) > 0)
			info = StringByKey("RECREATION", info, ":", ";")
			info = StringByKey("zColor(x)", info, "=", ";")
			if (strlen(info) > 2)
				info = info[1,strlen(info)-2]
				wname = StringFromList(0, info, ",")
				wave w = $wname
				ctab = StringFromList(3, info, ",")
				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
			endif
		endif
	endfor
	
	string images = ImageNameList(graphname, ";")
	n = ItemsInList(images, ";")
	for (i = 0; i < n; i += 1)
		objname = StringFromList(i, images, ";")
		wave w = ImageNameToWaveRef(graphname, objname)
		info = ImageInfo(graphname, objname, 0)
		if (strlen(info) > 0)
			info = StringByKey("RECREATION", info, ":", ";")
			info = StringByKey("ctab", info, "=", ";")
			if (strlen(info) > 2)
				info = info[1,strlen(info)-2]
				ctab = StringFromList(2, info, ",")
				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
			endif
		endif
	endfor

	setdatafolder save_df
end

/// k-space mapping of 2D angle-energy distribution (scienta image)
///
/// courtesy of F. Matsui
///
/// @param inwave	2D wave, x = kinetic energy (eV), y = polar angle (deg)
///						note: the kinetic energy is with reference to the vacuum level at the sample.
///						if the work functions of the analyser and the sample differ:
///						Ekin,sample = Ekin,analyser + WFanalyser - WFsample
///						where WFanalyser = Ephot - EFermi
///
/// @return the output wave has the name of the input wave with the suffix "_k".
///
Function AngleToK(inwave)
	Wave inwave
	String newname = NameofWave(inwave)+"_k"
	Duplicate/O inwave, $newname
	Wave outwave = $newname
	Variable rows,columns,xdelta,xoffset,ydelta,yoffset,kmin,kmax,Emax
	// inwave parameters
	rows = DimSize(inwave,0)
	columns = DimSize(inwave,1)
	xdelta = DimDelta(inwave,0)
	xoffset = DimOffset(inwave,0)
	ydelta = DimDelta(inwave,1)
	yoffset = DimOffset(inwave,1)
	Emax= xoffset + xdelta*(rows-1)
	kmin = 0.5123*sqrt(Emax)*sin(pi/180*(yoffset)) // calculate the k boundaries
	kmax = 0.5123*sqrt(Emax)*sin(pi/180*(yoffset+(columns-1)*ydelta))
	SetScale/I y kmin,kmax,"Ang^-1", outwave
	// scale the y axis
	outwave = interp2D(inwave, x, 180/pi*asin(y/ (0.5123*sqrt(x)))) // recalculate to k
	outwave = (NumType(outwave)==2) ? 0 : outwave // replace NaNs (optional)
End