clear;
load del13C.dat; % these are from Veizer's compilation
a=del13C(:,1); %
h=del13C(:,2); % these are the del 13C values
h=h+100; % a temporary offset so that all values are positive

% this routine finds the mean del18O values for each 1 myr year period
last=floor(max(a));

t=(0:last);
mh=zeros(size(t));
w=0;

for i=0:last-1;
    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 mh 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=t(q); % time vector corresponding to teh above

h=h-100; % restores it to original values
mh1=mh1-100; % ditto above
figure(1); 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);

