function map_model(filename,map_depth,showgrid,showevents)

% Extract index corresponding to chosen slice depth
load(filename,'zmod')
iz = find(zmod==map_depth);
if isempty(iz); error(['Invalid depth chosen for map slice. Available depths in km are: ', sprintf('%d, ', zmod)]) ; end

% Initialize figure
figure;
set(gcf,'Color','w');clf;
pos = get(gcf,'Position');
set(gcf,'Position',[pos(1),pos(2),1456,448])

maskdws=true; % Apply mask based on dws threshold

% Plot synthetic model
subplot(1,3,1)
vpvs_map(filename,iz,'Vp',maskdws,showgrid,showevents)
title(['Vp at ',num2str(map_depth),' km depth'],'fontsize',16)

% Plot recovered model
subplot(1,3,2);
vpvs_map(filename,iz,'Vs',maskdws,showgrid,showevents)
title(['Vs at ',num2str(map_depth),' km depth'],'fontsize',16)

subplot(1,3,3);
vpvs_map(filename,iz,'Vp/Vs',maskdws,showgrid,showevents)
title(['Vp/Vs at ',num2str(map_depth),' km depth'],'fontsize',16)

end

function vpvs_map(filename,iz,type,maskdws,showgrid,showevents)
global shorelines

% Load model data
switch type
    case 'Vp'
        load(filename,'Vp')
        V=Vp;
    case 'Vs'
        load(filename,'Vs')
        V=Vs;
    case 'Vp/Vs'
        load(filename,'Vp','Vs')
        V=Vp./Vs;
        
end
load(filename,'zmod','gridlat','gridlon')
if maskdws; load(filename,'dws');end
if showevents; load(filename,'evlat','evlon','evz','event_ids');end

% velocity slice
velslice=squeeze(V(:,:,iz));

% DWS slice
if maskdws
    % Set DWS threshold
    dwsThresh=0.02.*max(vec(dws(dws>3)));
    % Extract slice
    dwsslice=dws(:,:,iz);
end

% Initialize map axes
hm=worldmap([48.75,51], [-128.5,-125]);
setm(hm,'PLabelLocation',1,'MLabelLocation',1,'LabelUnits','degree','Plinelocation',3,'Grid','on','Plabelround',0,'Mlabelround',0,'FontSize',18)
setm(hm,'FFaceColor','w') % set the frame fill to white

% Interpolate on finer grid for better display
[glon,glat]=meshgrid(linspace(-130,-123.5,500),linspace(48,52,500));
velsliceq = griddata(gridlon,gridlat,velslice,glon,glat);
if maskdws
    dwssliceq = griddata(gridlon,gridlat,dwsslice,glon,glat);
    % Apply DWS threshold to get mask
    dwsmask = zeros(size(dwssliceq));
    dwsmask(dwssliceq>=dwsThresh)=1;
    % Smooth
    hfilt = fspecial('average', [15 15]);
    dwsmask = filter2(hfilt, dwsmask);
end

% Plot velocity map
h=pcolorm(glat,glon,velsliceq);
if maskdws % Apply transparency mask for dws values above threshold
    set(h,'AlphaData',dwsmask,'facealpha','flat','edgecolor','none');
end

% Set colormap and colorbar
colormap(jet(124))
cb=colorbar('Location','EastOutside','tag','colorbar');
switch type
    case 'Vp'
        colormap(gca,jet(124))
        cb=colorbar('Location','EastOutside','tag','colorbar');
        clo=5; cup=8;
        set(gca, 'CLim', [clo,cup])
        set(cb,'XTick', clo:cup,'fontsize',14);
        ylabel(cb,'Vp (km/s)','fontsize',14)
    case 'Vs'
        colormap(gca,jet(124))
        cb=colorbar('Location','EastOutside','tag','colorbar');
        clo=3; cup=5;
        set(gca, 'CLim', [clo,cup])
        set(cb,'XTick', clo:0.5:cup,'fontsize',14);
        ylabel(cb,'Vs (km/s)','fontsize',14)
    case 'Vp/Vs'
        load('RdBu_colormap.mat','cmap')
        colormap(gca,cmap)
        cb=colorbar('Location','EastOutside','tag','colorbar');
        clo=1.6; cup=1.85;
        set(gca, 'CLim', [clo,cup])
        set(cb,'XTick', 1.5:0.1:2,'fontsize',14);
        ylabel(cb,'Vp/Vs','fontsize',14)
end

% Plot shorelines
geoshow(shorelines,'FaceColor','none','EdgeColor','k')
lake = ([shorelines.Level]==2);
geoshow(shorelines(lake),'FaceColor','none','EdgeColor','k')

% Plot inversion grid nodes
if showgrid
    plotm(gridlat(:),gridlon(:),'go','linewidth',2);
end

% Plot relocated events
if showevents
    % Earthquakes
    ind=find(evz>=zmod(iz) & event_ids> 100000001);
    plotm(evlat(ind),evlon(ind),'ko','MarkerFaceColor','k','MarkerSize',3)
    % LFEs
    indlfe=find(evz>zmod(iz) & event_ids< 100000001);
    if ~isempty(indlfe)
        plotm(evlat(indlfe),evlon(indlfe),'kd','MarkerSize',6,'MarkerFaceColor','c')
    end
    % Blasts
    indblast=find(evz==0.01);
    if ~isempty(indblast)
        plotm(evlat(indblast),evlon(indblast),'rp','MarkerSize',6,'MarkerFaceColor','r')
    end
end

% Plot stations
load('stations.mat','slat','slon')
plotm(slat,slon,'ks','MarkerFaceColor','y','MarkerSize',8,'linewidth',1.5)

end