%Compute kinetic energy spectra from model data

fid = fopen('ng03-Kxyz500-simpleLaplace.dat','r','n');
chunk = fread(fid,'float');
chunk = reshape(chunk,305,305,30,3);
u = chunk(:,:,:,1);
v = chunk(:,:,:,2);
w = chunk(:,:,:,3);

nd = size(chunk,1);

field = u;

%Removing linear trend in field in both x and y directions...

for j=1:nd
    slope(j,:) = (field(nd,j,:) - field(1,j,:))./(nd-1);
end

slope = reshape(slope,1,305,30);
slope = repmat(slope,305,1);

for i=1:nd
    for j=1:nd
        fieldp(i,j,:) = field(i,j,:) - 0.5*(2*i - nd - 1)*slope(i,j,:);
    end
end

clear slope

for i=1:nd
    slope(i,:) = (fieldp(i,nd,:) - fieldp(i,1,:))/(nd-1);
end

slope = reshape(slope,1,305,30);
slope = repmat(slope,305,1);

for j=1:nd
    for i=1:nd
        fieldd(i,j,:) = fieldp(i,j,:) - 0.5*(2*j - nd - 1)*slope(i,j,:);
    end
end


%Computing 1D Fourier transform...x-direction for now
dx = 2;                             %grid spacing, km
N = size(fieldd,1);
%nfft = 2^nextpow2(N);               %number of samples in fft

%if mod(N,2) == 0
%    k=-N/2:N/2-1;                  %even number of samples
%else
%    k=-(N-1)/2:(N-1)/2;            %odd number of samples
%end

%half = floor(nfft/2);
half = floor(N/2);                  %truncate to integer

%amp = abs(fft(fieldd,nfft,1)/N);   %discrete 1D Fourier transform, points = nfft
amp = abs(fft(fieldd,[],1)/N);      %discrete 1D Fourier transform, points = N
zero = amp(1,:,:);                  %mean wavenumber or frequency                  
asym = 2*amp(2:half+1,:,:);         %asymmetric wavenumbers or frequencies
spec = cat(1,zero,asym);            %total spectrum


l = [0:1:half];
k = (2*pi*l)/(dx*N);                %wavenumber, km^-1
