#pragma rtGlobals=3		// Use modern global access method and strict wave access.
#pragma version = 1.8
#pragma IgorVersion = 6.2
#pragma ModuleName = PearlAnglescanProcess
#include "pearl-vector-operations"
#include "pearl-polar-coordinates"
#include <New Polar Graphs>

// copyright (c) 2013-17 Paul Scherrer Institut
//
// 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-17 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.


/// 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.
///
/// @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 (default) binomial, see Igor's Smooth operation
/// @arg 2 boxcar, see Igor's Smooth operation
/// @arg 3 scienta_ang_transm() function fit
/// @arg 4 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 = 1
	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 dist = ad_profile_x(strip, -inf, inf, "")
	variable div = mean(dist)
	dist /= div

	if (check)
		duplicate /o dist, check_dist
	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(dist, fit_params)
			dist = scienta_ang_transm(fit_params, x)
			break
		case 4:
			loess /smth=(smooth_factor) srcWave=dist
			break
	endswitch
	
	if (check)
		duplicate /o dist, check_smoo
	endif

	// divide
	if (check != 2)
		strip /= dist[p]
	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&degree;.
///
/// this normalization may be useful if the intensity varies with a 360&degree; 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
	wave dist = ad_profile_y(strip, -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
	wave dist = ad_profile_y(strip, -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 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
	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

/// 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
	if (strlen(source_nickname))
		string s_prefix = source_nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	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 theta1 = $s_theta
	wave polar1 = $s_polar
	wave azim1 = $s_azim
	wave tot1 = $s_tot
	wave weight1 = $s_weight
	wave values1 = $s_int
	wave /z matrix1 = $s_matrix
	
	variable npol = numpnts(theta1)
	
	setdatafolder dest_folder
	make_hemi_grid(npol, dest_nickname, xpdplot=xpdplot)
	
	// dest data
	if (strlen(dest_nickname))
		s_prefix = dest_nickname + "_"
		s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	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 theta2 = $s_theta
	wave polar2 = $s_polar
	wave azim2 = $s_azim
	wave tot2 = $s_tot
	wave weight2 = $s_weight
	wave values2 = $s_int

	tot2 = tot1
	weight2 = weight1
	values2 = values1
	if (waveexists(matrix1))
		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()

	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	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 polar = $s_polar
	wave azim = $s_azim
	wave tot = $s_tot
	wave weight = $s_weight
	wave 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

/// 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
///	@arg	0				linear
///	@arg	1				stereographic (default)
///	@arg	2				azimuthal
///	@arg	3				gnomonic (0 <= polar < 90).
///
/// @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
	
	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
	
	// hemi grid waves
	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	string s_matrix = s_prefix + "matrix"
	
	wave /z values = $s_int
	wave /z azim = $s_azim
	wave /z polar = $s_polar
	wave /z matrix = $s_matrix

	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)

	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
	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)
			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)
			draw_hemi_axes(graphname, do_grids=do_grids)
			break
	endswitch
	
	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
	
	if (strlen(nickname))
		string s_prefix = nickname + "_"
	else
		s_prefix = ""
	endif
	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	wave /z azim = $s_azim
	wave /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 an arbitrary angle scan 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,
/// and adds all data points which fall into a bin.
/// a point which lies exactly on the upper boundary falls into the next bin.
/// this function does not clear previous values before adding new data.
/// values are added to the _tot wave, weights to the _wt wave.
/// the intensity (_i) wave is calculated as _tot / _wt.
function hemi_add_anglescan(nickname, values, polar, azi, [weights])
	string nickname // name prefix of holo waves.
		// may be empty.
	wave values	// intensity values
		// the wave can be one- or two-dimensional.
		// no specific order required, the function sorts the arrays internally
	wave polar	// polar coordinates. allowed range 0 <= theta <= 90
		// dimensions corresponding to value.
	wave azi	// azimuthal coordinates. allowed range -360 <= phi < +360
		// dimensions corresponding to value.
	wave weights	// total accumulation time of each point of values. default = 1

	if (ParamIsDefault(weights))
		duplicate /free values, weights
		weights = 1
	endif
	
	// quick check whether hemi grid is existing
	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	string s_theta = s_prefix + "th"

	wave /z w_values = $s_int
	wave /z w_azim = $s_azim
	wave /z w_polar = $s_polar
	wave /z w_theta = $s_theta
	if (!waveexists(w_values) || !waveexists(w_azim) || !waveexists(w_polar))
		abort "Missing hemispherical scan grid. Please call make_hemi_grid() first."
	endif
	
	// make internal copies, one-dimensional, ordered in theta
	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
	
	variable pol
	variable pol_st = abs(w_theta[1] - w_theta[0])
	variable pol1, pol2
	
	duplicate /free azi_copy, azi_slice
	duplicate /free values_copy, values_slice
	duplicate /free weights_copy, weights_slice
	for (pol = 90; pol >= 0; pol -= pol_st)
		pol1 = pol - pol_st / 2
		pol2 = pol + pol_st / 2
		extract /free /indx polar_copy, sel, (pol1 < polar_copy) && (polar_copy <= pol2)
		if (numpnts(sel) > 0)
			redimension /n=(numpnts(sel)) azi_slice, values_slice, weights_slice
			azi_slice = azi_copy[sel]
			values_slice = values_copy[sel]
			weights_slice = weights_copy[sel]
			hemi_add_aziscan(nickname, values_slice, pol, azi_slice, weights=weights_slice)
		endif
	endfor
end

/// add an azimuthal scan 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,
/// and calculates the mean values of the data points which fall into a bin.
/// a point which lies exactly on the upper boundary falls into the next bin.
function hemi_add_aziscan(nickname, values, polar, azi, [weights])
	string nickname // name prefix of holo waves.
		// may be empty.
	wave values // intensity values of the azimuthal scan at the positions given in the azi parameter
	variable polar // polar angle where to add the azi scan
	wave azi // angle positions of the azimuthal scan
		// acceptable range: >= -360 and < +360
		// no specific order required, the function sorts the array internally
	wave weights	// total accumulation time of each point of values. default = 1

	if (ParamIsDefault(weights))
		duplicate /free values, weights
		weights = 1
	endif
		
	// hemi grid waves
	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	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 w_polar = $s_polar
	wave w_azim = $s_azim
	wave w_values = $s_int
	wave w_totals = $s_totals
	wave w_weights = $s_weights
	wave w_index = $s_index
	wave w_theta = $s_theta
	wave w_dphi = $s_dphi
	wave w_nphis = $s_nphis
	
	// destination slice coordinates
	//polar = round(polar)
	//variable ipol = 90 - polar
	variable ipol = BinarySearch(w_theta, polar)
	if (ipol < 0)
		abort "assertion failed in hemi_add_aziscan(): polar angle not found in grid."
	endif
	
	variable d1, d2
	if (ipol >= 1)
		d1 = w_index[ipol - 1]
	else
		d1 = 0
	endif
	d2 = w_index[ipol] - 1
	variable nd = d2 - d1 + 1
	variable dphi = w_dphi[ipol]
	variable az1, az2
	
	// source slice coordinates	
	// order the slice from -dphi/2 to 360-dphi/2
	azi = azi < 0 ? azi + 360 : azi
	azi = azi >= 360 - dphi/2 ? azi - 360 : azi
	duplicate /free values, sel_values
	duplicate /free weights, sel_weights
	
	// loop over destination
	variable id
	variable v1, v2, w1, w2
	for (id = 0; id < nd; id += 1)
		az1 = (id - 0.5) * dphi
		az2 = (id + 0.5) * dphi
		extract /free /indx azi, sel, (az1 <= azi) && (azi < az2)
		if (numpnts(sel) > 0)
			redimension /n=(numpnts(sel)) sel_values, sel_weights
			sel_values = values[sel]
			sel_weights = weights[sel]
			v1 = w_totals[d1 + id]
			w1 = w_weights[d1 + id]
			if ((numtype(v1) == 2) || (w1 <= 0))
				v1 = 0
				w1 = 0
			endif
			v2 = sum(sel_values)
			w2 = sum(sel_weights)
			w_totals[d1 + id] = v1 + v2
			w_weights[d1 + id] = w1 + w2
		endif
	endfor
	w_values[d1, d1 + nd - 1] = w_totals[p] / w_weights[p]
end

/// interpolate a hemispherical scan onto a rectangular grid
///
/// @warning experimental
/// this function has been tested for one specific set of scan parameters.
/// the interface and code may change at any time.
/// the function depends on the ster_x and ster_y waves that are created by display_hemi_scan.
///
function interpolate_hemi_scan(nickname)
	string nickname

	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif

	string s_polar = s_prefix + "pol"
	string s_azim = s_prefix + "az"
	string s_ster_x = s_prefix + "ster_x"
	string s_ster_y = s_prefix + "ster_y"

	wave values = $s_int
	wave azim = $s_azim
	wave polar = $s_polar
	wave ster_x = $s_ster_x
	wave 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]
	//ImageInterpolate /stw /s={x0, dx, xn, x0, dx, xn} voronoi triplet

	variable size = 181
	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
	
	matrix = (x^2 + y^2) < 4 ? matrix : nan
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
	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	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 theta1 = $s_theta
	wave polar1 = $s_polar
	wave azim1 = $s_azim
	wave tot1 = $s_tot
	wave weight1 = $s_weight
	wave 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)
		hemi_add_anglescan(nickname, intensity, theta, fold_phi)
		fold_phi = fold_phi >= 180 ? fold_phi + 360 / folding - 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

	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	string s_totals = s_prefix + "tot"
	string s_weights = s_prefix + "wt"
	string s_polar = s_prefix + "pol"

	wave w_polar = $s_polar
	wave w_values = $s_int
	wave w_totals = $s_totals
	wave 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

	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	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 w_polar = $s_polar
	wave w_azim = $s_azim
	wave w_values = $s_int
	wave w_totals = $s_totals
	wave w_weights = $s_weights
	wave w_index = $s_index
	wave w_theta = $s_theta
	wave w_dphi = $s_dphi
	wave 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
	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
	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

	if (strlen(nickname))
		string s_prefix = nickname + "_"
		string s_int = s_prefix + "i"
	else
		s_prefix = ""
		s_int = "values"
	endif
	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 w_polar = $s_polar
	wave w_azim = $s_azim
	wave w_values = $s_int
	wave w_totals = $s_totals
	wave w_weights = $s_weights
	wave w_index = $s_index
	wave w_theta = $s_theta
	wave w_dphi = $s_dphi
	wave 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)
		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
		return w_cut
	else
		return $""
	endif
end

static function check_contrast(values, pcmin, pcmax, vmin, vmax)
	wave values
	variable pcmin
	variable pcmax
	variable &vmin
	variable &vmax
	
	dfref save_df = GetDataFolderDFR()
	dfref dfr = NewFreeDataFolder()
	setdatafolder dfr
	StatsQuantiles /inan /iw /q /z values
	wave index = w_quantilesindex
	variable imin = round(numpnts(index) * pcmin / 100)
	variable imax = round(numpnts(index) * (100 - pcmax) / 100)
	vmin = values[index[imin]]
	vmax = values[index[imax]]
	KillDataFolder dfr
	setdatafolder save_df
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.
///
function set_contrast(pcmin, pcmax, [graphname, colortable])
	variable pcmin
	variable pcmax
	string graphname
	string colortable
	
	if (ParamIsDefault(graphname))
		graphname = ""
	endif
	if (ParamIsDefault(colortable))
		colortable = ""
	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
				endif
				check_contrast(w, pcmin, pcmax, vmin, vmax)
				ModifyGraph /w=$graphname zColor($objname)={w, vmin, vmax, $ctab, rev}
			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
				endif
				check_contrast(w, pcmin, pcmax, vmin, vmax)
				ModifyImage /w=$graphname $objname ctab={vmin, vmax, $ctab, rev}
			endif
		endif
	endfor

	setdatafolder save_df
end