clear all;

load DivExtOr.dat;
t=DivExtOr(:,1); % age in myr
ag=DivExtOr(:,2); % Number of All Genera
wr=DivExtOr(:,3); % Number of Well-Resolved Genera
e=DivExtOr(:,4); % Origination extinction %
o=DivExtOr(:,5); % Origination intensity %
d=DivExtOr(:,6); % WR generea detrended by removal of a cubic %

figure(1);
subplot(4,1,1); plot(t,ag); title('all genera');
subplot(4,1,2); plot(t,wr,t,d); title('well-resolved genera; detrended');
subplot(4,1,3); plot(t,e); title('extinction %');
subplot(4,1,4); plot(t,o); title('origination %');

% detrend
ag=detrend(ag);
wr=detrend(wr);
e=detrend(e);
o=detrend(o);
d=detrend(d);

%%

%calculate frequency array
Nyquist=0.5/abs(t(2)-t(1)); 
npt=2^16;
f=linspace(0,Nyquist,npt/2);
fmax=.05;
num=round(npt/2*fmax/Nyquist);
freq=f(1:num);

%calculate FFTs and Power after padding

H1=fft(ag,npt);
P1=H1.*conj(H1);
P1=P1/mean(P1(1:npt/2));

H2=fft(wr,npt);
P2=H2.*conj(H2);
P2=P2/mean(P2(1:npt/2));

H3=fft(e,npt);
P3=H3.*conj(H3);
P3=P3/mean(P3(1:npt/2));

H4=fft(o,npt);
P4=H4.*conj(H4);
P4=P4/mean(P4(1:npt/2));

H5=fft(d,npt);
P5=H5.*conj(H5);
P5=P5/mean(P5(1:npt/2));

%plot

figure(2);
plot(freq,P1(1:num),freq,P2(1:num),freq,P3(1:num),freq,P4(1:num),freq,P5(1:num))
title('Muller-Rohde Diversity Data','FontSize',14)
xlabel('frequency in cycles per myr')
ylabel('spectral power')
legend('all genera','well-resolved genera','% extinction','% origination','WR gen - cubic')
mark1x

%%

% Diversity Band-Passes


% bandpass filter to remove noise because there is a long-period trend to
% the data
%first define the range of frequencies to keep, between f1 and f2
f1=1/75; % 1/longer period
f2=1/40;% 1/shorter period
n=length(t);
fc=fft(wr);
fre=linspace(0,2*Nyquist,n)';
[m,k1]=min(abs(fre-f1));
[m,k2]=min(abs(fre-f2));
k3=n-k2+2;
k4=n-k1+2;
% zero all but the desired band
fc(1:k1)=zeros(k1,1);
fc(k2:k3)=zeros(k3-k2+1,1);
fc(k4:n)=zeros(k1-1,1);

%take inverse FFT; remove the imaginary part
newwr=real(ifft(fc));
figure(3); plot(t,wr,t,newwr); title('band-pass of well-resolved data');
%%


load Sr_ratio.dat;
a=Sr_ratio(:,1);
h=Sr_ratio(:,2);

% this routine finds the mean Sr ratios for each 0.1 myr year period

t1=(0:.1:520);
mh=zeros(size(t1));
w=0;

for i=0:.1:519;
    ind=find(i<a & a<=(i+.1));
    w=w+1;
    k=sum(h(ind));
    n=sum(h(ind)>0);
    mh(w+1)=k/n;
end

% the resulting vector nh which is the mean isotope ratio for each myr
% contains some NaN results if there are no data for that particular myr,
% so these have to be removed

q=find(mh>0); % indices of mh that are real numbers greater than 0
mh1=mh(q); % vecor of means of real numbers
t1=t1(q); % time vector corresponding to the above

figure(4); plot(t1,mh1,a,h);


% next comes a linear interpolation to evenly spaced data points

%interpolation to evenly spaced meters
t2=linspace(min(t1),max(t1),4*length(t1));
mh2=interp1(t1,mh1,t2);


% subtract mean
mh2=detrend(mh2);

%calculate frequency array
Nyquist=0.5/abs(t2(2)-t2(1)); 
npt=2^16;
f=linspace(0,Nyquist,npt/2);
fmax=.05;
num=round(npt/2*fmax/Nyquist);
freq=f(1:num);
%calculate FFT and Power after padding

H=fft(mh2,npt);
P=H.*conj(H);

%normalize to unit mean power
P=P/mean(P(1:npt/2));
%plot

figure(5);
plot(freq,P(1:num))
title('Sr isotopes','FontSize',14)
xlabel('frequency in cycles per myr')
ylabel('spectral power')
mark1x
%%

% bandpass filter to remove noise because there is a long-period trend to
% the data
%first define the range of frequencies to keep, between f1 and f2
f1=1/75; % 1/longer period
f2=1/40;% 1/shorter period
n=length(t2);
fc=fft(mh2);
fre=linspace(0,2*Nyquist,n)';
[m,k1]=min(abs(fre-f1));
[m,k2]=min(abs(fre-f2));
k3=n-k2+2;
k4=n-k1+2;
% zero all but the desired band
fc(1:k1)=zeros(k1,1);
fc(k2:k3)=zeros(k3-k2+1,1);
fc(k4:n)=zeros(k1-1,1);

%take inverse FFT; remove the imaginary part
newmh=real(ifft(fc));
figure(5); plot(t2,mh2,t2,newmh);

figure(6); plotyy(t2,newmh,t,newwr); title('Sr BP vs well-Resolved genera BP')

%%

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

figure(7);
plot(t3,c3);

c3=detrend(c3);

%calculate frequency array
npt=2^16;
Nyquist=0.5/abs((t3(2)-t3(1)));
f=linspace(0,Nyquist,npt/2);
fmax=.05;
num3=round(npt/2*fmax/Nyquist);
freq3=f(1:num3);

% bandpass filter to remove noise
%first define the range of frequencies to keep, between f1 and f2
f1=1/75; % 1/longer period
f2=1/40;% 1/shorter period
n=length(t3);
fc=fft(c3);
fre=linspace(0,2*Nyquist,n)';
[m,k1]=min(abs(fre-f1));
[m,k2]=min(abs(fre-f2));
k3=n-k2+2;
k4=n-k1+2;
% zero all but the desired band
fc(1:k1)=zeros(k1,1);
fc(k2:k3)=zeros(k3-k2+1,1);
fc(k4:n)=zeros(k1-1,1);

%take inverse FFT; remove the imaginary part
newc3=real(ifft(fc)); % invert the bp to go back to time domain
figure(8);
plot(t3,c3,t3,newc3);

%calculate FFT and Power after padding

Hn=fft(newc3,npt);
Pn=Hn.*conj(Hn); 

%normalize to unit mean power
Pn=Pn/mean(Pn(1:npt/2));


%figure(3);
%plot(freq3,Pn(1:num3));title('BP filtered Exxon sea level','FontSize',14);
%xlabel('frequency cycles per Myr')
%ylabel('spectral power')
%mark1x

%calculate FFT and Power after padding

H=fft(c3,npt);
P=H.*conj(H); 

%normalize to unit mean power
P=P/mean(P(1:npt/2));


figure(9);
plot(freq3,P(1:num3));title(' Exxon sea level','FontSize',14);
xlabel('frequency cycles per Myr')
ylabel('spectral power')
mark1x
