function [plon,plat,zm,zt] = read_moho_cmm(lon1,lat1,lon2,lat2,n)
% READ_MOHO_CMM plots 2D profiles of an estimated Moho surface
% 
% This function plots 2D profiles of the estimated Moho surface,
% California Moho Model, v.1.0 (CMM-1.0).
% REFERENCE: Tape-Plesch-Shaw-Gilbert SRL 2012, Electronic Seismologist
%
% Note that if you simply want to load the estimated Moho surface, you can
% use some of the commands below.
%
% The function calls functions that are within the Mapping Toolbox.
%
% INPUT:
%   lon1,lat1   starting point of profile
%   lon2,lat2   finishing point of profile
%   n           number of points along profile
% 
% OUTPUT:
%   plon,plat   profile points
%   zm          elevation of moho, km
%   zt          elevation of topography, km
% 
% EXAMPLES:
%   [plon,plat,zm,zt] = read_moho_cmm(-122,30.7,-112,37,150);
%   [plon,plat,zm,zt] = read_moho_cmm(-122,32,-116,39,150);
%   [plon,plat,zm,zt] = read_moho_cmm(-121,30.5,-112,35,150);
%

% directory containing the data files
ddir = './';

% estimated Moho surface CMM-1.0, saved on a regular lon-lat mesh
% (this means gridding is trivial with the 'reshape' command)
ifile = [ddir 'cal_moho01_q02_q08_ir02_id01_plot.dat'];
[xlon0,ylat0,zdep_km0] = textread(ifile,'%f%f%f');
melev_km = -zdep_km0;
nx = length(unique(xlon0));
ny = length(unique(ylat0));
Xm = reshape(xlon0,ny,nx);
Ym = reshape(ylat0,ny,nx);
Zm = reshape(melev_km,ny,nx);

disp('min/max for moho surface:');
dinfo(xlon0,ylat0,melev_km);

% ETOPO1 topography and bathymetry (Amante and Eakins, 2009)
% extracted onto a regular lon-lat mesh using GMT command
% grd2xyz ETOPO1_Ice_g.grd -R-128.05/-111.95/29.95/43.05 > ETOPO.xyz
ifile = [ddir 'ETOPO.xyz'];
[tlon,tlat,telev_m] = textread(ifile,'%f%f%f');
telev_km = telev_m/1000;
nx = length(unique(tlon));
ny = length(unique(tlat));
Xt = reshape(tlon,nx,ny)';
Yt = reshape(tlat,nx,ny)';
Zt = reshape(telev_km,nx,ny)';

disp('min/max for topo surface:');
dinfo(tlon,tlat,telev_km);

% compute profile of target points
%[plat,plon] = gcwaypts(lat1,lon1,lat2,lon2,n);
[plat,plon] = track2(lat1,lon1,lat2,lon2,[],'degrees',n);

% compute distance along profile, in km
d = deg2km(distance(repmat(plat(1),n,1),repmat(plon(1),n,1),plat,plon));

% find values of original model at CUVM gridpoints
zm = interp2(Xm,Ym,Zm,plon,plat);
zt = interp2(Xt,Yt,Zt,plon,plat);

ifig = 1;
if ifig==1 
    figure; hold on;
    plot(d,zt,'k','linewidth',3);
    plot(d,zm,'k','linewidth',3);
    %legend('Topo (ETOPO1)','Moho (CMM-1.0)','location','southwest');
    plot(d,zeros(n,1),'--','linewidth',2,'color',0.6*[1 1 1]);
    xlim([min(d) max(d)]);
    xlabel('Distance, km'); ylabel('Elevation, km');
    title(sprintf('%.1f km profile from (%.2f, %.2f) to (%.2f, %.2f)',d(end),lon1,lat1,lon2,lat2));

%     % list the vertical exaggeration (square region only)
%     ax1 = axis;
%     xran = ax1(2)-ax1(1);
%     yran = ax1(4)-ax1(3);
%     vex = xran/yran;
%     text(0.05*xran+ax1(1),0.05*yran+ax1(3),sprintf('H:V = %i',round(vex)),'fontsize',14);
%     axis square
    
    for kk=1:2
        figure; hold on;
        if kk==1
            pcolor(Xm,Ym,Zm); shading flat; title('Moho elevation, km'); colorbar;
        else
            pcolor(Xt,Yt,Zt); shading flat; title('Topography and Bathymetry, km'); 
        end
        colorbar; plot(plon,plat,'ko','markersize',6,'markerfacecolor','k');

        [socalfaults,ilabs,labels] = read_delimited_file('jennings_more.xy',2);
        plot(socalfaults(:,1),socalfaults(:,2),'k','linewidth',1);
        axis tight
    end
end

%--------------------------------------------------------------------------

function dinfo(x,y,z)
% list the extreme values

[~,iminz] = min(z);
[~,imaxz] = max(z);
disp(sprintf('minimum of %.3f km at (%.2f, %.2f)',z(iminz),x(iminz),y(iminz)));
disp(sprintf('maximum of %.3f km at (%.2f, %.2f)',z(imaxz),x(imaxz),y(imaxz)));

%--------------------------------------------------------------------------

function [data,ilabs,labels,inds] = read_delimited_file(filename,ncol)
% read in fault file -- perhaps dlmread provides a simpler option

% read all lines
lines = textread(filename,'%s','delimiter','\n','whitespace','');
nlines = length(lines);

% read in the LUR file and string delimiters
jj = 0; kk = 0;
data = NaN * ones(nlines,ncol);
for ii = 1:nlines
    tline = lines(ii);
    temp = str2num(char(tline));
    if ~isempty(temp)
        jj = jj + 1;
        data(ii,:) = temp;
    else
        kk = kk + 1;
        labels(kk) = tline;
        ilabs(kk) = ii;
    end
end
ilabs = ilabs(:);
labels = labels(:);

% indices for each segment
inds = [ilabs(1:end-1)+1 ilabs(2:end)-1];

%==========================================================================