%%

% This analyzes the exxon sea level curve 
load exxon.dat;   % input the time series
t=exxon(:,2);% Myr
sl=exxon(:,3);% sea level

figure(1);
plot(t,sl);
title('Vail/Haq Sea Level Curve')
ylabel('m.s.l.m')
xlabel('Age (Ma)');

%% FFT

[P,f] = dbfft(t,sl,.05,2^16);

figure(2);
plot(f,P);
xlabel('frequency cycles per Myr')
ylabel('spectral power')
mark1x
%% Band-pass

newsl=sqbandpass(t,sl,30,1);

figure(3);
plot(t,sl,t,newsl);
%% FFT of band-pass
[P2,f2] = dbfft(t,newsl,.5,2^16);
figure(4);
plot(f2,P2);
xlabel('frequency cycles per Myr')
ylabel('spectral power')
mark1x