% this file has a demonstration of several steps needed for this lab.
%  (1) plotting  data on a map 
%  (2) drawing small circles (based on rotations about a pole)
%  (3) calculate distances between points on a map
%  (4) fit a line through a set of points


% things you need to do
%  (1) add data for other hot spot tracks
%  (2) make maps showing small circles through other hot spot tracks
%  (3) make age distance plots for other hot spot tracks (in addition to
%       Hawaii

% for "fun" and totally optional
%  (4) save your best tracks to a kml file and import into GOOGLE EARTh to
%       plot.  Type "help pwr_kml" to find out how to do this
%   (5) make an "oblique mercator" projection using your pole of rotation


%The file containing all seamount data from the paper is in the folder
% I'll put some notes at the bottom about this


% Here is the age location data for the lab:


% from: Ages of seamounts, islands and plateaus on the Pacific plate
%         version 2.1 - 15 May 2004
%      Valérie Clouard(1) and Alain Bonneville(2)

% Long.Lat.   Age(Ma) Uncertainty   island/seamount dating method original source
data=[
174.30 31.80  38.70 0.90   %Abbott Hawaii-Emperor K/Ar [20]
175.90 30.90  38.60 0.30   %Colahan Hawaii-Emperor K/Ar [20]
172.30 32.10  42.40 2.30   %Daikakuji Hawaii-Emperor K/Ar [21]
203.80 20.60  0.75 0.04     %Haleakala Hawaii-Emperor K/Ar [22]
171.20 38.40  55.40 0.90   %Jingu Hawaii-Emperor K/Ar [23]
203.40 20.50  1.02 0.18     %Kahoolawe Hawaii-Emperor K/Ar [22]
200.50 22.00  5.10 0.20     %Kauai Hawaii-Emperor K/Ar [24]
204.75 19.40  0.00 0.01     %Kilauea Hawaii-Emperor observation
171.60 33.70  39.90 1.20   %Kimmei Hawaii-Emperor K/Ar [21]
204.25 20.10  0.43 0.02     %Kohala Hawaii-Emperor K/Ar [25]
171.67 35.10  48.10 0.80    %Koko Hawaii-Emperor K/Ar [21]
193.70 23.60  12.00 0.40    %La Pérouse Pinnacle Hawaii-Emperor K/Ar [26]
203.00 20.80  1.28 0.04      %Lanai Hawaii-Emperor K/Ar [27]
188.00 25.70  19.90 0.30    %Laysan Hawaii-Emperor K/Ar [28]
204.74 18.93  0.00 0.01                %Loihi Hawaii-Emperor observation
204.50 19.75  0.38 0.05      %Mauna Kea Hawaii-Emperor K/Ar [29]
182.70 28.30  27.70 0.60    %Midway Hawaii-Emperor K/Ar [30]
203.30 21.20  1.83 0.07      %Molokai Hawaii-Emperor K/Ar [22]
195.50 23.50  10.30 0.40    %Necker Hawaii-Emperor K/Ar [26]
198.00 23.00  7.20 0.30      %Nihoa Hawaii-Emperor K/Ar [26]
170.20 41.20  56.20 0.60     %Nintoku Hawaii-Emperor K/Ar [23]
188.00 25.30  26.60 2.70     %Northampton Bank Hawaii-Emperor K/Ar [28]
202.20 21.40  2.60 0.10       %Oahu Hawaii-Emperor K/Ar [31]
170.30 37.50  55.20 0.70     %Ojin Hawaii-Emperor K/Ar [23]
184.10 27.90  20.60 0.50    %Pearl and Hermes reef Hawaii-Emperor K/Ar [32]
170.00 44.00  64.70 1.10     %Suiko Hawaii-Emperor K/Ar [23]
204.00 20.80  1.32 0.04       %West Maui Hawaii-Emperor K/Ar [31]
172.00 32.70  43.40 1.60     %Yuryaku Hawaii-Emperor K/Ar [32]
];

% this gets the "size" of data (how many rows(np) and columns (nc)).
[np,nc]=size(data);


% lets make a map of the north Pacific and plot the locations in the table
% on it.  This is the first item on the list


% start your map:
% remember you can change any of the arguments below to "improve" your map
figure(1)
clf
m_proj('miller','lat',[10 60],'long',[160 210]);
m_coast('patch',[.7 1 .7],'edgecolor','none');
m_grid('box','fancy','linestyle','none');
m_elev('contour',[-10000:200:4000])


% turns out that to plot individual points using m_map we need to use a loop
% there are np points to plot so do a loop from 1 to np
for i=1:np
   m_line(data(i,1),data(i,2),'Marker','square','markersize',6,'color','r'); 
end
%try other marker types or sizes to "improve" your map


% the double percent sets off a cell.  You can execute cells independently
%%

% now lets draw a small circle starting from the current big Island of
% Hawaii (it is the 3rd row in the data)

% first calculate a small circle based on a guessed axis of rotation:
% our guess for the lattitude and longitude of the pole:
% you are free to change these numbers and the ones here are not 'optimal'.
%  you can do a guess and check or do some reading to find what poles
%  others have used.

pole_lat=59;
pole_long=-140;


% you need a vector of angles to rotate through: (minus to rotate in right
% direction

% you ought to form the angle vector in a physical way - a rotational
% velocity (degrees per million years) time an age
% since you know that the bend in the chain is at about 43 Ma.
 ang_vel=1;  % this is 1 degree per million years
 rot_age=0:43;  % this vector would be the ages along the small circle 
                % ending at an age that is about where the bend occurs
 ang_vec=-ang_vel*rot_age;

% starting point (long and lat) of rotation at Kilauea Hawaii
long_start=data(8,1);
lat_start=data(8,2);


% here is the command that gives us a set of points along a small circle
% generated by an "Euler rotation"
[long_out, lat_out]=Euler_Rotation(pole_long,pole_lat,ang_vec,long_start,lat_start);
 
% we are working in an awkward area (around the break between E and W
% longitudes.  This causes some difficulty for m_map.  Here is a fix.
% Insure that the output of Euler_Rotation is consistent with m_map:

long_out= 360-(long_out+180);  % can you see what this does?  -)

%now plot a line through points on the small circle
m_line(long_out,lat_out,'color','r','linewi',2); 


% The default values above give a small circle that neither follows the
% chain adequately nor extends for the full distance of the chain.  You can
% change the pole location and the angular velocity to get a better fit.
% It is up to you to decide how to judge the quality of the fit.  Should
% the small circle go exactly through all points? How far misfit is
% acceptable?  Change numbers and re-run just this cell to see new results
% while keeping the old.

 %%    
     

%  On to #3.How to calculate distances between points on a map?
% m_map has a function.  What we will do is calculate the distance
% between Kilauea and other older volcanoes.  we will need a loop
% lets just deal with data out to the "bend"

% this command find indeces of the data that are younger than 43 Ma
id =find(data(:,3)<43);


% range is going to be a vector of distances between  Kilauea and other
% dated islands and seamounts
range=zeros(size(id));
    for i=1:length(id)      
      range(i)=m_lldist([data(8,1) data(id(i),1)],[data(8,2) data(id(i),2)]); 
    end


% for #4, fitting a line: Let x be the ages and y the locations.  
%   Prblems to think about: The table has a single location but we know 
%that the volcanoes are typcially active at the same time over a region 
%of 10s of km.  Then at one location the volcano may be active over a long
%period of time. 
 
xdata=range;
ydata=data(id,3);
uncert_ydata=data(id,4);  % here are dating uncertainties for our error analysis

% The following does a linear fit to data having uncertain values for y.
% Output are intercepts and slopes and their uncertainties
[intercept,slope,uncert_intercept,uncert_slope]=linreg(xdata,ydata,uncert_ydata)
% note that the slope here is in Ma per km (slowness or inverse velocity).  A  linear velocity
% Note that the fitted uncertainty is also in slowness.  To determine the
% plate velocity uncertainty, determine the percentage uncertainty in the
% slowness and then apply that percentage uncertainty to the inverse of the
% slowness. 


%plot the results:
figure(6)
xc=[1 3400];
yc=intercept+slope*xc;
plot(xc,yc)
hold on
errorbar(xdata,ydata,uncert_ydata,'o')
hold off
%needs labels and title.
xlabel('Distance (km)')
ylabel('Age (Ma')
title('My first x-y plot')  % change to something meaningfull

% note that for Hawaii, there is a good age-distance progression  - farther
% away as age increases.

% Note also that the data beyond 15 Ma begin to scatter more. What might
% that mean?  If the error bars are smaller than their misfit from the line, what does
% that mean?

%You need to make age distance plots for other hot spot tracks - If they are
%the same age as Hawaii and if plates are rigid and hot spots are fixed,
%then the same euler pole should work for the same rotational velocity and
%the age progression should be consistent.


%%
% to load data for all dated seamounts and islands type

 load PacSeamountAges
 
 %this gives you five vectors of data
 %  (1) Averageage
 %  (2) uncertainty in age
 %  (3) latitude
 %  (4) longitude
 %  (5) name of feature
 
 % you will need to consult a map and choose which feature should be
 % associated with which hot spot tracks.
 



    