#pragma rtGlobals=3		// Use modern global access method and strict wave access.
#include "pearl-area-profiles"

// $Id$

/// @file
/// @brief tools for analysing the Fermi edge measured by the Scienta EW4000 analyser.
/// @ingroup ArpesPackage
///
/// proposed procedure
///
/// * angular normalization
/// * fit curved fermi function
/// * calculate corrected energy coordinates and map to single independent variable
/// * fit normal fermi function
///
/// @author matthias muntwiler, matthias.muntwiler@psi.ch
/// @author thomas dienel
///
/// @copyright 2013-15 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


function analyse_curved_edge(data)
	wave data	// 2D counts data, x-scale = kinetic energy, y-scale = analyser angle
	variable guess_EF	// initial guess of the Fermi level
	
	dfref savedf = GetDataFolderDFR()
	dfref datadf = GetWavesDataFolderDFR(data)
	setdatafolder datadf
	
	make /n=5 /d /o w_coef_int
	make /n=7 /d /o w_coef_curved

	// 1) integrate data
	wave xint = ad_profile_x(data, -inf, inf, "", noavg=1)
	duplicate /free xint, xint_sig
	xint_sig = sqrt(xint * 4) / 4
	variable xmin = wavemin(xint)
	variable xmax = wavemax(xint)
	variable xmean = mean(xint)

	// 2) fit integrated data for reference
	w_coef_int[0] = xmin
	w_coef_int[1] = 0
	w_coef_int[2] = xmax - xmin
	w_coef_int[3] = dimoffset(xint, 0) + dimdelta(xint, 0) * dimsize(xint, 0) / 2
	w_coef_int[4] = 100
	FuncFit /NTHR=0 FermiFuncLinDOS w_coef_int xint /D /I=1 /W=xint_sig
	wave w_sigma
	duplicate /o w_sigma, w_sigma_int
	
	// 3) normalize data
	wave yavg = ad_profile_y(data, -inf, inf, "")
	variable ymean = mean(yavg)

	duplicate /o data, data_norm, data_sig
	data_sig = sqrt(data * 4) / 4
	data_norm = data_norm * ymean / yavg[q]
	data_sig = data_sig * ymean / yavg[q]
	
	// 4) fit normalized data
	wave xavg = ad_profile_x(data, -inf, inf, "")
	xmin = wavemin(xavg)
	xmax = wavemax(xavg)
	xmean = mean(xavg)
	
	w_coef_curved[0] = {0, 0, 1, 95.5, 100, 0, -0.0001}
	w_coef_curved[0] = xmin
	w_coef_curved[2] = xmax - xmin
	w_coef_curved[3] = w_coef_int[3]
	w_coef_curved[4] = w_coef_int[4]
	
	//variable box = min(11, numpnts(xavg) / 5)
	//FindLevel /B=(box) /EDGE=2 /Q xavg, (xmin + xmax) / 2
	//if (v_flag == 0)
	//	w_coef_curved[3] = v_levelx
	//else
	//	w_coef_curved[3] =  dimoffset(data, 0) + dimdelta(data, 0) * dimsize(data, 0) / 2
	//endif
	
	duplicate /o data_norm, fit_data_norm
	FuncFitMD /X=1 /NTHR=0 FermiFuncLinDOS_2Dcorr w_coef_curved  data_norm /D=fit_data_norm /I=1 /W=data_sig
	wave w_sigma
	duplicate /o w_sigma, w_sigma_curved

	display /k=1; appendimage data_norm
	ModifyImage data_norm ctab= {xmin,xmax,Grays,0}
	AppendMatrixContour fit_data_norm
	
	setdatafolder savedf
	return 0
end

function record_results(index)
	variable index
	
	dfref savedf = GetDataFolderDFR()
	wave /sdfr=root: Tint, Tint_sig, Tcurv, Tcurv_sig, EFcurv, EFcurv_sig
	wave w_coef_int, w_coef_curved
	wave w_sigma_int, w_sigma_curved
	
	Tint[index] = w_coef_int[4]
	Tint_sig[index] = w_sigma_int[4]
	Tcurv[index] = w_coef_curved[4]
	Tcurv_sig[index] = w_sigma_curved[4]
	EFcurv[index] = w_coef_curved[3]
	EFcurv_sig[index] = w_sigma_curved[3]
	
	setdatafolder savedf
end

function integrate_curved_edge(data, data_sig)
	wave data
	wave /z data_sig
	//wave coef
	
	string name = nameofwave(data)
	
	duplicate /o data, data_out
	redimension /n=(dimsize(data,0)) data_out
	data_out = 0
	
	if (waveexists(data_sig))
		duplicate /o data_sig, sig_out
		redimension /n=(dimsize(data,0)) sig_out
		sig_out = 0
	endif

	wave ywgt = ad_profile_y(data, -inf, inf, "", noavg=1)
	ywgt = 1 / ywgt / 4
	//ywgt = 1 / 4
	variable sum_ywgt = sum(ywgt)
	
	variable nx = dimsize(data, 0)
	variable ny = dimsize(data, 1)
	variable iy
	variable yy
	variable dx
	variable dy
	variable dp
	variable dp_min = 0
	
	wave PassEnergy = :attr:PassEnergy
	wave NumSlices = :attr:NumSlices
	
	for (iy = 0; iy < ny; iy += 1)
		dy = dimoffset(data, 1) + dimdelta(data, 1) * iy
		dy = dy / dimdelta(data, 1)
		dx = slit_shift(dy * 902 / NumSlices[0], PassEnergy[0])
		dp = round(dx / dimdelta(data, 0))	// <= 0
		dp_min = min(dp_min, dp)
		
		data_out[] = p+dp >= 0 ? data_out + data[p+dp][iy] * ywgt[iy] : data_out

		if (waveexists(data_sig))
			sig_out = p+dp >= 0 ? sig_out + data_sig[p+dp][iy] * ywgt[iy] : sig_out
		endif
	endfor
	
	data_out /= sum_ywgt
	data_out[0, -dp_min][] = nan

	if (waveexists(sig_out))
		sig_out /= sum_ywgt
		//sig_out[0, -dp_min] = nan
	endif
end

function slit_correction(data, data_out, epass)
	wave data	// must be image with original dimensions (no cropping!)
	wave /z data_out	// 2D or 1D wave to receive the result
		// X dimension must be identical to the one of data
		// if 2D, Y dimension must be either identical to the one of data
		// if 1D, the result will be the sum of the corrected slices
	variable epass	// pass energy
	
	if (!WaveExists(data_out))
		string name = nameofwave(data) + "_corr"
		duplicate /o data, $name
		wave data_out = $name
		//redimension /n=(dimsize(data,0)) data_out
	endif
	data_out = 0
	
	variable nx = dimsize(data, 0)
	variable ny = dimsize(data, 1)
	variable iy
	variable yy
	variable dx
	variable dy
	variable dp
	variable dp_min = 0
	
	for (iy = 0; iy < ny; iy += 1)
		dy = dimoffset(data, 1) + dimdelta(data, 1) * iy
		dy = dy / dimdelta(data, 1)
		dx = slit_shift(dy * 902 / ny, epass)
		dp = round(dx / dimdelta(data, 0))	// <= 0
		dp_min = min(dp_min, dp)
		
		if (wavedims(data_out) >= 2)
			data_out[][iy] = p+dp >= 0 ? data_out + data[p+dp][iy] : data_out
		else
			data_out = p+dp >= 0 ? data_out + data[p+dp][iy] : data_out
		endif
	endfor
	
	data_out[0, -dp_min][] = nan
end

//------------------------------------------------------------------------------
threadsafe Function FermiFuncLinDOS2D_corr(w,x,y) : FitFunc
// linear density of states below Fermi level
// 2D data with corrections:
// - straight slit (slit shift)
// - transmission function (polynomial)
//------------------------------------------------------------------------------
	Wave w; Variable x; variable y
	// w[0] = background far above the fermi level
	// w[1] = slope of the linear background
	// w[2] = amplitude
	// w[3] = fermi level in eV
	// w[4] = temperature in K
	// w[5] = transmission - linear term
	// w[6] = transmission - quadratic term
	// w[7] = transmission - cubic term
	// w[8] = pass energy (hold this value)
	// w[9] = y scale: pixels / unit of y

	variable pos = w[3] + slit_shift(y * w[9], w[8])
	variable transm = 1 + w[5] * y + w[6] * y^2 + w[7] * y^3
	variable fermi = (w[1] * min(x - pos, 0) + w[2]) / ( exp( (x - pos)  /  (kBoltzmann * w[4]) ) + 1.0  )
	return transm * (fermi + w[0])
end

//------------------------------------------------------------------------------
threadsafe Function FermiFuncLinDOS_2Dcorr_old(w,x,y) : FitFunc
// linear density of states below Fermi level
// 2D data with a polynomial shift of the Fermi level in the second dimension
//------------------------------------------------------------------------------
	Wave w; Variable x; variable y
	// w[0] = background far above the fermi level
	// w[1] = slope of the linear background
	// w[2] = amplitude
	// w[3] = fermi level in eV
	// w[4] = temperature in K
	// w[5] = shift - linear term
	// w[6] = shift - quadratic term

	variable pos = w[3] + w[5] * y + w[6] * y^2
	return w[0] + (w[1] * min(x - pos, 0) + w[2]) / ( exp( (x - pos)  /  (kBoltzmann * w[4]) ) + 1.0  )
end

/// MCP radius seen by the camera in pixels
static constant mcp_radius_pix = 555
/// physical size (radius) of the MCP in mm
static constant mcp_radius_mm = 20
/// physical size (radius) of the hemisphere in mm
static constant hemi_radius_mm = 200
/// energy range imaged on MCP relative to the pass energy
static constant mcp_radius_epass = 0.04

threadsafe function slit_shift(ypix, epass)
	// calculates the energy shift due to straight slit
	// the radius of the curve is 1/2 of the hemisphere radius
	variable ypix	// vertical (angle/position) pixels distance from center
		// = slice coordinate * 902 / number of slices
	variable epass // pass energy
	
	variable rpix = mcp_radius_pix * hemi_radius_mm / mcp_radius_mm
	//variable rpix = mcp_radius_pix * hemi_radius_mm / 2 / mcp_radius_mm
	variable rene = epass * mcp_radius_epass * hemi_radius_mm / 2 / mcp_radius_mm
	
	variable isin = asin(ypix / rpix)
	variable dene = rene * (cos(isin) - 1)
	
	return dene
end

function show_shift(data)
	wave data
	variable epass
	
	variable ny = dimsize(data, 1)
	make /o /n=(ny) shift_x
	setscale /i x -ny/2, ny/2, "", shift_x
	wave PassEnergy = :attr:PassEnergy
	wave NumSlices = :attr:NumSlices
	
	shift_x = slit_shift(x * 902 / NumSlices[0], PassEnergy[0])
	
	variable ecenter = dimoffset(data, 0) + dimdelta(data, 0) * dimsize(data, 0) / 2
	shift_x += ecenter
end

/// calculate the energy resolution of the analyser
///
/// @param	epass	pass energy in eV
/// @param	slit		analyser entrance slit in mm
///
/// @return energy resolution (FWHM)
function analyser_energy_resolution(epass, slit)
	variable epass
	variable slit

	variable respow
	if (epass < 4)
		respow = 1110
	elseif (epass < 8)
		respow = 1400
	else
		respow = 1750
	endif
	
	return epass * max(0.2, slit) / 0.2 / respow	
end