% this is a test script to show you how to access information in Muller's
% compilation of ocean age.  The data base is in the form of a large
% matrix. Each line has a longitude, lattitude, age (Ma) and age
% uncertainty (Ma).  Each line increments lattitudes and longitudes in .1
% degree steps.  Here are the first four lines of the matrix:         
%         -180           90       55.569           15
%       -179.9           90       55.569           15
%       -179.8           90       55.568           15
%       -179.7           90       55.568           15      
% There are a total of 6,485,401 lines.  So this is a big matrix that will
% fill up a lot of computer memory if you load it and any search of it is
% going to take some time.   I give below an example of how you can search
% the data base to find ages and age uncertainties for a collection of
% lattudes and longitudes.

% Here are some locations where I want to know the ages:
 latdata=[ 40 40 40 40 40 40 ];
 longdata=[-130 -129 -128 -127 -126 -125 ];
 npts=length(latdata);  %need to know the number of points I have
 
 % now set up a couple variables to put ages and uncertainties into
 agedata=zeros(1,npts);
 uncertdata=agedata;
 
% Now I need to load the big file.  Expect this to be slow:

load ages

% the search assumes that the lats and longs are given only to the nearest
% 0.1 degree so that there is an exact match.  The loop searches the data
% base for each of your points:

 for i=1:npts,
     id=find(age(:,1)==longdata(i) & age(:,2)==latdata(i));
     agedata(i)=age(id,3);
     uncertdata(i)=age(id,4);
 end
 
 %lets calculate the distance between our lattitudes and longitudes:
 
dist=m_lldist(longdata,latdata);  % this calculates the ditance between the points
dist=[0; cumsum(dist)];  % this adds up the distances to give the total distance from the 1st point

%  Now plot with age uncertainty errobars:

errorbar(dist,agedata,uncertdata,'o')

% lets calculate a fit through these points:

[a,b,da,db]=linreg(dist,agedata,uncertdata);

% Now calculate this line (two points only needed for line)

ycalc=a + b*dist([1 end]);

hold on
plot(dist([1 end]),ycalc)
hold off


% OK we have calculated a slope as age/distance but we really want
% dist/age.  Easy enough to take 1/b to get the apparent spreading rate in
% km/Ma (which is the same as mm/y).  But what about the uncertainty?
% You cannot just take 1/db.  The easiest way to do this is the following:

rate=abs(1/b)  % take absolute value to make rate positive
 uncertrate= 1/b - 1/(b+db)  
 
 

% You will notice that for this somewhat random collection of data, there
% is a lot of scatter and large error bars.  It is not surprizing that the
% the uncertainty is large for this example.  If we decided that there is
% just one "bad" point that should not be fit, we can redo the fit:

[a,b,da,db]=linreg(dist([1 3 4 5 6]),agedata([1 3 4 5 6]),uncertdata([1 3 4 5 6]));


ycalc=a + b*dist([1 end]);

hold on
plot(dist([1 end]),ycalc)
hold off


% OK we have calculated a slope as age/distance but we really want
% dist/age.  Easy enough to take 1/b to get the apparent spreading rate in
% km/Ma (which is the same as mm/y).  But what about the uncertainty?
% You cannot just take 1/db.  The easiest way to do this is the following:

rate=abs(1/b)  % take absolute value to make rate positive
uncertrate= 1/b - 1/(b+db)  

% this fits the remaing points better but because the age is so uncertain,
% the rate remains uncertain for this very short profile.


 

