% 2007 version.

% read the height data, and put in vector "z"
fid = fopen('500_hght.dat');
z = fscanf(fid,'%g',inf);
fclose(fid);

Nx = length(z); Nxm = Nx/2;

% take FFT 
fz = fft(z);

% notes on fz:
% fz(1) is the mean.
% fz(2:181) are the Fourier coeficients for zonal wavenumbers 
% 1--181.
% fz(182:360) are complex conjugates of fz(2:181), in reverse 
% order.
%

plotfz = fz(2:181);

% plot with a logarithmic y-axis because the range of values 
% is large.

semilogy(abs(plotfz));
title('Fourier Amplitude Spectrum')
xlabel('Planetary Wavenumber')
ylabel('Fourier Amplitude')

%---------------------------------------------------------
% filter by planetary wavenumber
s = [0 1:Nxm -(Nxm-1):1:-1]';

% long waves (s <= 4)
sf = abs(s)<=4;

%
% end of hints for this part!
%

% fourier derivative operator
lat = 40;
g = 9.8;
f = 2*7.292e-5*sind(lat);
km = 1000.;
Lx = 6380.*km*2*pi*cosd(lat);

% spectral derivative operator
%d = 

%
% end of hints for this part!
%

% do it again, but with grid point derivative operator
dx = Lx/Nx;
%D = 

%
% end of hints for this part!
%