%This is a script for making various types of plots (mesh, surface, sclaed,
%gradded, etc.).  In this case using crustal thickness variations etc.  It
%calls on data 'ears_overview.txt'.  Call this data using the code written
%by Kevin Eagar called ImportCSV(filename)
%
%Wendy Bohon, March 12-14, 2009

for n = 1:length(Column03)
    if isempty(Column03{n})
        NewColumn03(n,1) = NaN;
    else
        NewColumn03(n,1) = str2num(Column03{n});
    end
end

for n = 1:length(Column04)
    if isempty(Column04{n})
        NewColumn04(n,1) = NaN;
    else
        NewColumn04(n,1) = str2num(Column04{n});
    end
end


for n = 1:length(Column05)
    if isempty(Column05{n})
        NewColumn05(n,1) = NaN;
    else
        NewColumn05(n,1) = str2num(Column05{n});
    end
end


for n = 1:length(Column06)
    if isempty(Column06{n})
        NewColumn06(n,1) = NaN;
    else
        NewColumn06(n,1) = str2num(Column06{n});
    end
end


for n = 1:length(Column07)
    if isempty(Column07{n})
        NewColumn07(n,1) = NaN;
    else
        NewColumn07(n,1) = str2num(Column07{n});
    end
end


for n = 1:length(Column08)
    if isempty(Column08{n})
        NewColumn08(n,1) = NaN;
    else
        NewColumn08(n,1) = str2num(Column08{n});
    end
end


for n = 1:length(Column09)
    if isempty(Column09{n})
        NewColumn09(n,1) = NaN;
    else
        NewColumn09(n,1) = str2num(Column09{n});
    end
end


for n = 1:length(Column10)
    if isempty(Column10{n})
        NewColumn10(n,1) = NaN;
    else
        NewColumn10(n,1) = str2num(Column10{n});
    end
end


for n = 1:length(Column11)
    if isempty(Column11{n})
        NewColumn11(n,1) = NaN;
    else
        NewColumn11(n,1) = str2num(Column11{n});
    end
end


for n = 1:length(Column12)
    if isempty(Column12{n})
        NewColumn12(n,1) = NaN;
    else
        NewColumn12(n,1) = str2num(Column12{n});
    end
end

A = [NewColumn03 NewColumn04 NewColumn05 NewColumn06 NewColumn07 NewColumn08 NewColumn08 NewColumn09 NewColumn10 NewColumn11];

s = size (A)
n = s(1)
m=s(2)
latmin = 31;
latmax = 42;
longmin = -120;
longmax = -105;
k = 0
B=[]

for i = 1:n
    if (A(i, 1) > latmin && A (i, 1) < latmax) && (A(i, 2) > longmin && A(i, 2)<longmax)
        k = k+1
        B(k,:) = A (i, :);
    end
end

NewColumn03 = B(:, 1);
NewColumn04 = B(:, 2);
NewColumn05 = B(:, 3);
NewColumn06 = B(:, 4);
NewColumn07 = B(:, 5);
NewColumn08 = B(:, 6);
NewColumn09 = B(:, 7);
NewColumn10 = B(:, 8 );
NewColumn11 = B(:, 9);
NewColumn12 = B(:, 10);

%Plots the original data
figure(1)
clf
plot (NewColumn04, NewColumn03, 'k+')
xlabel ('Longitude')
ylabel('Latitude')
title ('Seismic Stations')
print('seismicstations.jpg','-djpeg');

%Clears the figure and makes the new one
figure(2)
clf

%mesh plot for triangles
%delaunay function computes the triangles
tri = delaunay(NewColumn04, NewColumn03);

%trimesh function plots the triangles computed by delaunay
trimesh(tri,NewColumn04,NewColumn03)
colormap([0 0 0])
view(2)
hold on

%labels the axis and names the image
xlabel('Longitude')    
ylabel('Latitude')   
title('Seismic Station Mesh Plot')
axis equal
print('meshmap_seismicstation.jpg','-djpeg');


%grids the data
%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command

x = NewColumn04 
y = NewColumn03
h = NewColumn06 

%This section will make a Contour Map w Data and Seismic Stations

%this brings up "help" and gives lots of instructions on griddata
help griddata

%grids the data
%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command
numsteps=100
xstep = abs((min(x)-max(x))/numsteps);
ystep = abs((min(y)-max(y))/numsteps);
xx = [min(x):xstep:max(x)];
yy = [min(y):ystep:max(y)];


%This grids the data.
%the n' means transpose N (switch it from a row vector to a column vector)
% This was orginial  [xi,yi,zi] = griddata(e,n,hh,ee,nn');
[xi,yi,zi] = griddata(x,y,h,xx,yy');

clf
figure(3)
%this was orginial [c,h] = contour(xi,yi,zi,v);
[c,h] = contour(xi,yi,zi);

%this labels the contours
clabel(c,h);
hold on
%show data:
plot(x,y,'o')
%show grid:
plot(xi,yi,'+')
%labels axis etc.
xlabel('Longitude')    
ylabel('Latitude')   
title('Contour Map w Data and Seismic Stations')
axis equal
print('CMap__Data_SS.jpg','-djpeg');

%This section will make a Contour Map of Crustal Thickness w Seismic Stations

%grids the data
%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command
numsteps=10
xstep = abs((min(x)-max(x))/numsteps);
ystep = abs((min(y)-max(y))/numsteps);
xx = [min(x):xstep:max(x)];
yy = [min(y):ystep:max(y)];

%THIS PART IS STILL MESSING UP.
%This grids the data.
%the 'n means transpose N (switch it from a row vector to a column vector)
% This was orginial  [xi,yi,zi] = griddata(e,n,hh,ee,nn');
[xi,yi,zi] = griddata(x,y,h,xx,yy');

%I have no idea what this does, or if I need it
%v = floor(min(hh)):ci:max(hh);
clf
figure(4)
%this was orginial [c,h] = contour(xi,yi,zi,v);
[c,h] = contour(xi,yi,zi);

%this labels the contours
clabel(c,h);
hold on
%show data:
plot(x,y,'o')
%labels axes
xlabel('Longitude')    
ylabel('Latitude')   
title('Contour Map of CT w Seismic Stations')
axis equal
print('CMap_CT_SS.jpg','-djpeg');

%grids the data
%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command




%This section will make a Contour Map w of V(p)/V(s) Data and Seismic
%Stations

%This defines the variables
x = NewColumn04 
y = NewColumn03
p = NewColumn08 

%this brings up "help" and gives lots of instructions on griddata
help griddata

%grids the data
%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command
numsteps=10
xstep = abs((min(x)-max(x))/numsteps);
ystep = abs((min(y)-max(y))/numsteps);
xx = [min(x):xstep:max(x)];
yy = [min(y):ystep:max(y)];


%This grids the data.
%the 'n means transpose N (switch it from a row vector to a column vector)
% This was orginial  [xi,yi,zi] = griddata(e,n,hh,ee,nn');
[xi,yi,zi] = griddata(x,y,p,xx,yy');

clf
figure(5)
%this was orginially [c,h] = contour(xi,yi,zi,v);
[c,p] = contour(xi,yi,zi);

%this labels the contours
clabel(c,p);
hold on
%show data:
plot(x,y,'o')
%show grid:
plot(xi,yi,'+')

xlabel('Longitude')
ylabel('Latitude')   
title('Contour Map of V(p)/V (s) w Data and Seismic Stations')
axis equal
print('CMap_V(p)/V(s)_Data_SS.jpg','-djpeg');

%This section will make a Contour Map w of V(p)/V(s) w/ seismic stations

%This defines the variables
x = NewColumn04 
y = NewColumn03
p = NewColumn08 

%this brings up "help" and gives lots of instructions on griddata
help griddata

%grids the data
%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command
numsteps=10
xstep = abs((min(x)-max(x))/numsteps);
ystep = abs((min(y)-max(y))/numsteps);
xx = [min(x):xstep:max(x)];
yy = [min(y):ystep:max(y)];

%This grids the data.
%the 'n means transpose N (switch it from a row vector to a column vector)
% This was orginial  [xi,yi,zi] = griddata(e,n,hh,ee,nn');
[xi,yi,zi] = griddata(x,y,p,xx,yy');

clf
figure(6)
%this was orginially [c,h] = contour(xi,yi,zi,v);
[c,p] = contour(xi,yi,zi);


%this labels the contours
clabel(c,p);
hold on
%show data:
plot(x,y,'o')
xlabel('Longitude')    
ylabel('Latitude')   
title('Contour Map of V(p)/V (s) w Seismic Stations')
axis equal
print('CM_Vps_ss.jpg','-djpeg');

%SCALED IMAGES
%Scaled image of Crustal thickness
figure(7)
clf
x = NewColumn04 
y = NewColumn03
h = NewColumn06 

%this brings up "help" and gives lots of instructions on griddata
help griddata

%grids the data
%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command
numsteps=100
xstep = abs((min(x)-max(x))/numsteps);
ystep = abs((min(y)-max(y))/numsteps);
xx = [min(x):xstep:max(x)];
yy = [min(y):ystep:max(y)];


%This grids the data.
%the n' means transpose N (switch it from a row vector to a column vector)
% This was orginial  [xi,yi,zi] = griddata(e,n,hh,ee,nn');
[xi,yi,zi] = griddata(x,y,h,xx,yy');


imagesc(xx',yy',zi)

hold on

%show data:
plot(x,y,'k.')

xlabel('Longitude')    
ylabel('Latitude')   
title('Scaled Image of Crustal Thickness')
axis equal 
colorbar

print('scalemap_CT.jpg','-djpeg');

%Scaled image of V(p)/V(s) with Seismic Stations
figure(8)
clf
x = NewColumn04 
y = NewColumn03
h = NewColumn08 

%this brings up "help" and gives lots of instructions on griddata
help griddata

%grids the data
%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command
numsteps=100
xstep = abs((min(x)-max(x))/numsteps);
ystep = abs((min(y)-max(y))/numsteps);
xx = [min(x):xstep:max(x)];
yy = [min(y):ystep:max(y)];


%This grids the data.
%the n' means transpose N (switch it from a row vector to a column vector)
% This was orginial  [xi,yi,zi] = griddata(e,n,hh,ee,nn');
[xi,yi,zi] = griddata(x,y,h,xx,yy');

imagesc(xx',yy',zi)

hold on

%show data:
plot(x,y,'k.')

xlabel('Longitude')    
ylabel('Latitude')   
title('Scaled Image of V(p)/V(s)')
axis equal 
colorbar

print('scalemap_Vps_ss.jpg','-djpeg');

%MAKES SURFACE MAPS
%surface map of V(p)/V(s)
surface
figure(10)
clf

x = NewColumn04 
y = NewColumn03
h = NewColumn08 

surfl(xx,yy,zi)
shading interp
colormap(gray)
hold on

plot3(x,y,h,'b.')
view(-30,60)
xlabel('Longitude')    
ylabel('Latitude')
zlabel('V(p)/V(s)')   
title('Surface Map of V(p)/V(s)')
axis equal

print('surfacemap_Vps.jpg','-djpeg');

%MAKES SURFACE MAP of crustal thickness
figure(9)
clf
surfl(xx,yy,zi)
shading interp
colormap(gray)
hold on

plot3(x,y,-h,'b.')
view(-30,60)
xlabel('Longitude')    
ylabel('Latitude')
zlabel('Crustal Thickness')   
title('Surface Map of Crustal Thickness')
axis equal

print('surfacemap_CT.jpg','-djpeg');