% GLG 410 Fall 2011
% Yellowstone Seismicity Script

% Christy Baker
% September 28, 2011

% This script will plot seismicity data for Yellowstone.

% Load text file
%----------------------------------------------------------------
myData = load ('YellowstoneSeismicitydata.txt', '-ascii');


EqYear = myData (:,1);
EqDate = myData (:,2);
EqLat = myData (:,3);
EqLong = myData (:,4);
EqDepth = myData (:,5);
EqMag = myData (:,6);

%------------------------------------------------------------------
% Creat plot of earthquakes and histogram displaying frequency

plot (EqLong, EqLat, 'r.');
title ('Yellowstone Seismicity 1973 - 2009');
figure (1)
ylabel ('Lattitude (degrees)');
xlabel ('Longitude (degrees)');
% %-------------------------------------------------------------------
% % %Histogram of Earthquakes per year
% % figure (2)
% % title ('Yellowstone Seismicity 1973 - 2009');
% % edges = 1972:1:2010;
% % n = histc(EqYear,edges);
% % bar(edges,n,'histc')
% % bar(edges,n,'histc');
% % xlabel('Year')
% % ylabel('Frequency')
% % title('Yellowstone Earthquake Histogram');
% % %--------------------------------------------------------------------
% % % Seismic activity in Yellowstone illustrates two events where earthquakes were more active:  1973 and 2008.
% % % Between '73 and '08, earthquake activity lulled and the seismic data reveals that earthquakes
% % % were common around two general clusters.

hold on
stepsize = 0.25;
long_range=min(EqLong):stepsize:max(EqLong);
lat_range=min(EqLat):stepsize:max(EqLat);
for long = long_range
    for lat = lat_range
        plot(long,lat,'bo')
         tf = EqLong <= long+stepsize./2 & EqLong >= long-stepsize./2 & EqLat <= lat+stepsize./2 & EqLat >= lat-stepsize./2;
       locations = find(tf);
       mean_depth = mean(EqDepth(locations));
       s = sprintf('%2.1f',mean_depth);
       text(long, lat, s);
    end
end




