%%%%%% load data from ac.bin
a=fopen('ac.bin');
b=fread(a,'double');

%%%%%%%%%%%%%%% parameters
tsteps=20001;
nmodes=2;
dt=5.89664e-12;

%reshape
y=reshape(b,length(b)/tsteps,tsteps);
%make it even
y=y(:,1:end-1);

%compute frequencies
dimt=tsteps-1;
invL=1./(dt*dimt);
dimd=ceil(dimt/2);
f=zeros(2*dimd,1);
for t=1:dimd
  f(t)=invL*(t-1);
end
for t=1:dimd
  f(t+dimd)=invL*(-dimd-2+t);
end

%fft
yf=fft(y,[],2)/dimt;

% 2D
am=yf(1:nmodes,:);
bm=yf(nmodes+1:nmodes*2,:);
py=2*abs(real(sum(am.*conj(bm))));

plot(f/3e8,py);axis tight;xlim([0 5]);