% DON'T RUN THIS WHOLE SCRIPT.  
% It contains separate things for you to work with.  Look for the cells and
% run one cell at a time.

% m_map has access to a higher resolution database of elevations and
% bathymetry (at 2 minute resolution - an even higher resolution is
% available but not necessary for this lab).  The command:
%[elevations,longitudes,latitudes]=m_etopo2([lonMIN lonMAX latMIN latMAX])
%returns matrixes of points on the earth within the box of minimum and
%maximum longitudes and latitudes.  To get an elevation profile across a
%subduction zone along the coast of Sumatra at the equator, here is an example (this is
%essentially an east-west profile:


%%
% here you load two databases - I got these from the NEIC and the
% Smithsonian.  The lab writup has the web link.  I went ahead and helped
% you by saving some data in the MATLAB format.  When you load these files,
% you will find a couple "ReadMe's" in the list of variables - type who to
% see what variables have been loaded, type the readme names to see what
% data you have.

load IndoEQdata
load VolcanoData

%%

figure(1)
lonMIN=90;
lonMAX=105;
latMIN=0;
latMAX=.1;

[elev,lon,lat]=m_etopo2([lonMIN lonMAX latMIN latMAX]);

%this returns 4x451 matrixes (many more longitude points since the latitude
%range is small)
% to plot the first row (at the latitidue 0), first covert to distances
% from the starting point (same step as in the first lab:

dels=m_lldist(lon(1,:),lat(1,:));
dist=[0; cumsum(dels)];
plot(dist,elev(1,:))
xlabel('Distance east from 90E (km)')
ylabel('Elevation/Bathymetry (m)')

%Note that you can find the trench, the forearc ridge and basin, the
%volcanic arc, and the back arc in this profile.  I expect you to make some
%profiles of your own at other places.  You also would need to show where
%these profiles are made on your base map.



%%

% This cell assumes you have loaded the data base information above.

%here is an example of making a 3D plot to show earthquakes and volcanos
%along one segment of the subduction zone.  you can change the search
%criteria and view other sections.  Click on the rotate icon in the figure
%window to rotate the image to view from all angles.

%here is the search box:
LatMin=-12;
LatMax=0;
LongMin=110;
LongMax=114;

% LongMin=100;  % this shows a bigger chunk of the arc.
% LongMax=120;

% There are a bunch of earthquakes that have a "default" depth of 33 km -
% these are events that were not located in depth and were just assigned a
% "default" depth.  They are annoying in the plots so I give you a way to
% remove them - see below

% the database has almost a million events, whic is slow to work with so,
% by setting a threshold it is possible to work with only larger events.
% Play with this (up or down) to see what works for you
EQThresh=5.5;


%Here is the search (each row is a single event, so the search is over all
%elements of the various columns
% see how we "find" things that meet criteria that are connected by the
% logical "and" which in matlab is "&" 

id=find(EQdata(:,1)>LatMin & EQdata(:,1)<LatMax & EQdata(:,2)>LongMin & EQdata(:,2)<LongMax & EQdata(:,4)>EQThresh &  EQdata(:,3)~=33);
idv=find(volcanoLatLon(:,1)>LatMin & volcanoLatLon(:,1)<LatMax & volcanoLatLon(:,2)>LongMin & volcanoLatLon(:,2)<LongMax);

%Now do a 3-D plot:
figure(1)
plot(EQdata(id,1),-EQdata(id,3),'o')
plot3(EQdata(id,1),EQdata(id,2),-EQdata(id,3),'o')
hold on
plot(volcanoLatLon(idv,1),volcanoLatLon(idv,2),'ro')
hold off


%This does not have any coastlines, you might want to add them or plot
%these earthquakes separately on an m_map

% this is just one place in the Indonesian area - you proably ought to look
% at others.


%%
% This is a fun plot - the whole area - you can look at the big picture
%  here, You can use toolbar tools on the
% figure window to zoom in or rotate or move over different region.
% You can also change the limits of latitude and longitude to look at just
% one area at a time
lonMIN=80;
lonMAX=160;
latMIN=-20;
latMAX=40;

% chose the volcanos in the area of interest, the earthquakes in the
% database were chosen for the same area.

idv=find(volcanoLatLon(:,1)>latMIN & volcanoLatLon(:,1)<latMAX & volcanoLatLon(:,2)>lonMIN & volcanoLatLon(:,2)<lonMAX);

%Now do a 3-D plot:
figure(2)
clf

% I have another simple to use function that plots coastlines on any plot.
% You need to set some variable:
args.file='world';
args.map_limits=[ lonMIN lonMAX latMIN latMAX];
% one more thing that matlab allows are "structures" - a varibale name that
% has sub variables embedded in it.  Here the structure "args" is set to contain two
% variables - "file" which tells us which datafile (a mat file) to use for making maps -
% this is the whole world - latter you will use a higher resolution file
% that has only the Pacific Northwest.  and the variable "map_limits" contains
% those as a vector. 

% here is the function that makes the base map of the chosen area
mapp(args);

%now hold that map and overlay the other stuff:
hold on

%plot the volcanos
plot(volcanoLatLon(idv,2),volcanoLatLon(idv,1),'ro')

% Below, notice that I divide the depth(in km) by 111.  Since the horizontal axis is in
% degrees of latitude and longitude and I want the map to have no vertical
% exageration, I just converted the depth to degrees - in a final
% publication I'd edit the vertical scale tic marks to show the true depth (you can
% do that using the figure editor, but it is not required for your report (but
% the caption must say what the depth range is in km.)

% plot3 does a 3D plot of vectors X,Y,and Z:
plot3(EQdata(:,2),EQdata(:,1),-EQdata(:,3)/111,'o')


%the following force the plot to have no vertical exageration and to remain
% fixed as you rotate the plot to look at it from various directions. 
% Try commenting this out to see what it looks like with lots of vertical
% exageration - sometimes that is a useful thing to do.
axis equal
axis vis3d
hold off

%There is another useful command:
%        view(az,el) 
%  where az is the azimuth of your viewing direction and el is
%  the angle of the view (90 is the default straight down view, 0 for on
%  edge).  If you find a view you like by rotating by hand, you can force
%  it the next time you plot by using this command. When you hold the
%  button down over the map the azimuth and elevation are shown on the
%  bottom left of the plot.


%%

% here is a short script to put the earthquakes on an m_map.

%Careful about doing it to all data since it will really tie up your
%computer and use up memory.  You might try plottin a smaller number of
%earthquakes.

figure(3)
clf 
m_proj('miller','lat',[-20 40],'long',[80 160]);

set(gca,'color',[.9 .99 1]);     
m_elev('contour',[-10000:500:8000])
m_coast('patch',[.2 1 .7],'edgecolor','k');
m_grid('box','fancy','linestyle','none');

[npts,nc]=size(EQdata);

for i=1:1000,  %1:npts   just plotting the first 1000 instead of all npts points.
    m_line(EQdata(i,2),EQdata(i,1),'Marker','square','markersize',2,'color','r'); 
end



%%



% here is a short script to put the earthquake focal mechanisms on an m_map.


% NOTE IT WILL NOT RUN UNTIL YOU FILL IN SOME VALUES  You can fiddle with
% other parameters as well - ie looking at one area in more detail or
% investigating different depth ranges of earthquakes (the third column of
% the data.  or change the size and color of the beachballs


figure(3)
clf 
m_proj('miller','lat',[-25 45],'long',[75 170]);
m_coast('patch',[.2 1 .7],'edgecolor','k');
m_grid('box','fancy','linestyle','none');


load IndoFMdata   % I downloaded all focal mechanisms in the NEIC catalog for the 
%                    study area (as an excel file) and then loaded it into
%                    MATLAB and saved it for you as a .mat file (only saved
%                    the minimum stuff (location, depth, size, and
%                    mechanisms)

% OK you will need to uncomment the following line and decide what range of
% earthquake size to map.  If you plot all the earthquakes, the map will be
% too confusing with too many small earthquakes that have random fault
% directions. If you use only the very larges, you may have too few to see
% trends. 

%idsize=find(FMdata(:,4)> SOME Magnitude value for a lower cutoff);

fm=FMdata(idsize,5:7);  % these are the three numbers for the focal mechanism
centerX=FMdata(idsize,2); % the lats and longs
centerY=FMdata(idsize,1);
np=length(centerX);
diam=1.5;  % the diameter of the symbol
ta=1;  % a flag that tells the function to do the plot on the map
color='r'; % color of your beachball

for i=1:np
beachball(fm(i,:), centerX(i), centerY(i), diam, ta, color)
end
%%

