function [nux,nuy,ampx,ampy]=findfreq(fadx,fady)
% function [nux,nuy,ampx,ampy]=findfreq(fadx,fady)
%
% this function calculates the tunes and the amplitude of the peak
% for given transverse oscillation data (turn by turn - FAD) by
% using an FFT with sine window and interpolation.
%
% return values:
% nux, nuy 	the two transverse betatron tunes
% ampx, ampy	the peak amplitude
%
% input values:
% fadx		horizontal FAD data (number of BPMs and number of turns is variable)
% fady		vertical FAD data
%
% Christoph Steier, July 1999

sfx=size(fadx);
sfy=size(fady);

mfadx = fadx-ones(sfx(1),1)*mean(fadx);
mfady = fady-ones(sfy(1),1)*mean(fady);

fftx = fft(mfadx);
ffty = fft(mfady);

fftx = abs(fftx);
ffty = abs(ffty);

% minx=int32(sfx(1)/20);
maxx=int32(sfx(1)/2);
% miny=int32(sfy(1)/20);
maxy=int32(sfy(1)/2);
% dminx=double(minx);
dmaxx=double(maxx);
% dminy=double(miny);
dmaxy=double(maxy);

[dummy,indx] = max(fftx(1:dmaxx,:));
[dummy,indy] = max(ffty(1:dmaxy,:));

dmaxx = indx + 5;
dminx = indx - 5;
dmaxy = indy + 5;
dminy = indy - 5;

sfadx = mfadx.* (ones(sfx(2),1)*sin(pi*[0:1/(sfx(1)-1):1]))';
sfady = mfady.* (ones(sfy(2),1)*sin(pi*[0:1/(sfy(1)-1):1]))';

fftx = fft(sfadx);
ffty = fft(sfady);

fftx = abs(fftx);
ffty = abs(ffty);

[dummy,indx] = max(fftx(dminx:dmaxx,:));
[dummy,indy] = max(ffty(dminy:dmaxy,:));

indx = indx + dminx;
indy = indy + dminy;

idx = [];
ampx = []; ampx2 = [];

for n = 1:sfx(2)

	if indx(n) == 1
		indx1=2; indx3=2;
	elseif indx(n) == (sfx(1)/2)
		indx1=(sfx(1)/2-1); indx3=(sfx(1)/2-1);
	else
		indx1=indx(n)-1; indx3=indx(n)+1;
	end

	if (fftx(indx3,n)>fftx(indx1,n))
		ampx(n) = fftx(indx(n),n); ampx2(n) = fftx(indx3,n);
		idx(n) = indx(n);
	else
		ampx(n) = fftx(indx1,n); ampx2(n) = fftx(indx(n),n);
		if (indx(n) ~= 1)
			idx(n) = indx1;
		else
			idx(n) = 0;
		end
	end
end

nux = 1/(sfx(1)) * (idx-1 + (2*ampx2./(ampx+ampx2)) - 0.5);

idy = [];
ampy = []; ampy2 = [];

for n = 1:(sfy(2))

	if indy(n) == 1
		indy1=2; indy3=2;
	elseif indy(n) == (sfy(1)/2)
		indy1=(sfy(1)/2-1); indy3=(sfy(1)/2-1);
	else
		indy1=indy(n)-1; indy3=indy(n)+1;
	end

	if (ffty(indy3,n)>ffty(indy1,n))
		ampy(n) = ffty(indy(n),n); ampy2(n) = ffty(indy3,n);
		idy(n) = indy(n);
	else
		ampy(n) = ffty(indy1,n); ampy2(n) = ffty(indy(n),n);
		if (indy(n) ~= 1)
			idy(n) = indy1;
		else
			idy(n) = 0;
		end
	end
end

nuy = 1/sfy(1) * (idy-1 + (2*ampy2./(ampy+ampy2)) - 0.5);

ampx = 2 * diag(fftx(indx,:),0)' ...
./( sin(pi*(indx-1+0.5-nux*sfx(1)))./(pi*(indx-1+0.5-nux*sfx(1))) + ...
sin(pi*(indx-1-0.5-nux*sfx(1)))./(pi*(indx-1-0.5-nux*sfx(1))));
ampy = 2 * diag(ffty(indy,:),0)' ...
./( sin(pi*(indy-1+0.5-nuy*sfy(1)))./(pi*(indy-1+0.5-nuy*sfy(1))) + ...
sin(pi*(indy-1-0.5-nuy*sfy(1)))./(pi*(indy-1-0.5-nuy*sfy(1))));

% ampx = 2 * diag(fftx(indx,:),0)' ...
% ./( sinc(0.5+nux*1024) + sinc(-0.5+nux*1024));
% ampy = 2 * diag(ffty(indy,:),0)' ...
% ./( sinc(0.5+nuy*1024) + sinc(-0.5+nuy*1024));

% nux = 1-nux;
% nuy = 1-nuy;