function [ ] = Func_Seismic_multiplot_MTEX4(C,rho,plot_name,plot_x,Vp_polar,Vs1_polar,Vs2_polar,ref_axis_1,ref_axis_2,ref_axis_3)
% Seismic_multiplot
% plot standard seismic anisotropy figures for publications
% 
% David Mainprice (20/11/2014)
%
% Tested with MTEX 4.0.6
%
% MTEX4 optimization thanks to Ralf Heilscher
%
% Input
% C - MTEX object elastic stiffness tensor
%
% rho - density (gm/cm3)
%
% plot_name - plot name (can include path '/MatLab_Programs/Test_plot.eps')
%
% plot_x : 1= X to north, 2= X to East, 3= X to South, 4= X to West
%
% Vp_polar : 1= plots Vp polarizations in black if 0 does not plot them on Vp polefigure
%
% Vs1_polar : 1= plots Vs1 polarizations in black if 0 does not plot them on AVs polefigure
%
% Vs2_polar : 1= plots Vs2 polarizations in white if 0 does not plot them on AVs polefigure
%
% ref_axis_1,ref_axis_2,ref_axis_3 = text reference axes labels e.g.;
% 'X','Y','Z' or '[100]','[010]','[001]'
%
%%
%**************************************************************************
% Plotting section
%**************************************************************************
% plotting convention
% set the default plot direction of the X-axis
%
% default north
plotx2north;
%
switch plot_x
    case 1
    plotx2north;
    case 2
    plotx2east;
    case 3
    plotx2south;
    case 4
    plotx2west;
    otherwise
    error('Default - use plotx2north')
end
%% Plotting section
% close all open graphics
close all
%
%**************************************************************************
% Setup multiplot
%**************************************************************************
% set colour map to seismic color map
% convection red = slow velocity and blue = fast velocity: red2blueColorMap
setMTEXpref('defaultColorMap',red2blueColorMap);
% max and min marker size
Marker_size_max_min = 20;
blackMarker = {'Marker','s','MarkerSize',Marker_size_max_min,...
  'MarkerEdgeColor','white','MarkerFaceColor','black'};
whiteMarker = {'Marker','o','MarkerSize',Marker_size_max_min,...
  'MarkerEdgeColor','black','MarkerFaceColor','white'};
% linewidth for polarization
linewidth_pol = 3;
% some global options for the titles
% set font size
global_fontSize = 40;
titleOpt = {'FontSize',global_fontSize,'visible','on'};
% reset default fontsize for larger numbers on colorbar
setMTEXpref('FontSize',30)
% Setup multiplot
% define plot size [origin X,Y,Width,Height]
mtexFig = mtexFigure('position',[0 0 1000 1000]);
% set up spacing between subplots default is 10 pixel
%mtexFig.innerPlotSpacing = 20;
%
% Standard Seismic plot with 8 subplots in 3 by 3 matrix
%
% Plot matrix layout
%        1 Vp        2 AVs      3 S1 or/& S2 polarizations
%        4 Vs1       5 Vs2      6 dVs
%        7 Vp/Vs1    8 Vp/Vs2
%
%% Computation of velocities and polarisations

% Generate velocities and polarizations on a fine pole figure grid
v_grid = plotS2Grid(C,'complete','upper');

[vp,vs1,vs2] = velocity(C,v_grid,rho);

% The polarizations we need only at a finer grid
p_grid = quiverGrid(C,'complete','upper');
% polarization = partical movement
[~,~,~,pp,ps1,ps2] = velocity(C,p_grid,rho);

%%
%**************************************************************************
% Vp : Plot P-wave velocity (km/s)
%**************************************************************************

% Plot P-wave velocity (km/s)
contourf(v_grid,vp,'parent',mtexFig.gca)
mtexTitle('Vp (km/s)',titleOpt{:})
% plot Vp polarizations
if(Vp_polar==1)
  hold on
  quiver(p_grid,pp,'linewidth',linewidth_pol,'color','black','parent',mtexFig.gca)
  hold off
end
% mark maximum with black square and minimum with white circle
[Vp_max_value,index_max] = max(vp(:));
Vp_max_vector = v_grid(index_max);

[Vp_min_value,index_min] = min(vp(:));
Vp_min_vector = v_grid(index_min);

% Anisotropy percent
AVp=200*(Vp_max_value-Vp_min_value)/(Vp_max_value+Vp_min_value);

AVpS = ['Vp Anisotropy (%) = ',num2str(AVp,'%6.1f')];

% N.B. x and y reversed in subplot
if (plot_x==1||plot_x==3)
   ylabel(AVpS,titleOpt{:})
end
if (plot_x==2||plot_x==4) 
   xlabel(AVpS,titleOpt{:})
end
% mark maximum with black square and minimum with white circle
hold on
plot(Vp_max_vector,blackMarker{:},'parent',mtexFig.gca)
plot(Vp_min_vector,whiteMarker{:},'parent',mtexFig.gca)
hold off
%
%**************************************************************************
% AVS : Plot S-wave anisotropy percentage for each propagation direction
% defined as AVs = 200*(Vs1-Vs2)/(Vs1+Vs2)
%**************************************************************************
%
% calculate S-wave anisotropy AVs for all propagation directions
AVs = 200*(vs1-vs2)./(vs1+vs2);
% create a new axis
mtexFig.nextAxis
% Plot S-wave anisotropy (percent)
contourf(v_grid,AVs,'parent',mtexFig.gca);
mtexTitle('S-wave anisotropy',titleOpt{:})

% find maximum and minimum values
[AVs_max_value,index_max] = max(AVs(:));
AVs_max_vector = v_grid(index_max);

[~,index_min] = min(AVs(:));
AVs_min_vector = v_grid(index_min); 

% Max percentage anisotropy
AVsS = ['Max Vs Anisotropy (%) = ',num2str(AVs_max_value,'%6.1f')];

if (plot_x==1||plot_x==3)
   ylabel(AVsS,titleOpt{:})
end
if (plot_x==2||plot_x==4) 
   xlabel(AVsS,titleOpt{:})
end
% mark reference axes
hold on
text([xvector,yvector,zvector],{ref_axis_1,ref_axis_2,ref_axis_3},...
  'backgroundcolor','w','FontSize',25,'parent',mtexFig.gca);

% mark maximum with black square and minimum with white circle
hold on
plot(AVs_max_vector,blackMarker{:},'parent',mtexFig.gca)
plot(AVs_min_vector,whiteMarker{:},'parent',mtexFig.gca)
hold off
%
%**************************************************************************
% S1 Polarization: Plot S-wave (Vs1 or Vs2) polarization directions
%**************************************************************************
%
% create a new axis
mtexFig.nextAxis
% plot
contourf(v_grid,AVs,'parent',mtexFig.gca);
mtexTitle('Vs1 polarization',titleOpt{:})
if(Vs2_polar==1 && Vs2_polar==0 )
  mtexTitle('Vs1 polarization',titleOpt{:})
end
if(Vs2_polar==0 && Vs2_polar==1 )
  mtexTitle('Vs2 polarization',titleOpt{:})
end
if(Vs2_polar==1 && Vs2_polar==1 )
  mtexTitle('Vs1 and Vs2 polarizations',titleOpt{:})
end
% plot Vs2 polarizations
if(Vs2_polar==1)
  hold on
  quiver(p_grid,ps2,'linewidth',linewidth_pol,'color','white','parent',mtexFig.gca)
  hold off
end
% plot Vs1 polarizations
if(Vs1_polar==1)
  hold on
  quiver(p_grid,ps1,'linewidth',linewidth_pol,'color','black','parent',mtexFig.gca)
  hold off
end

%
%**************************************************************************
% Vs1 : Plot Vs1 velocities (km/s)
%**************************************************************************
%
% create a new axis
mtexFig.nextAxis
% plot
contourf(v_grid,vs1,'parent',mtexFig.gca);
mtexTitle('Vs1 (km/s)',titleOpt{:})
% plot Vs1 polarizations
hold on
quiver(p_grid,ps1,'linewidth',linewidth_pol,'color','black','parent',mtexFig.gca)
hold off

% S1-wave velocity (km/s)
% Maximum and Minimum values and their directions
[Vs1_max_value,index_max] = max(vs1(:));
Vs1_max_vector = v_grid(index_max);

[Vs1_min_value,index_min] = min(vs1(:));
Vs1_min_vector = v_grid(index_min);

% Anisotropy percent
AVs1=200*(Vs1_max_value-Vs1_min_value)/(Vs1_max_value+Vs1_min_value);
% Percentage anisotropy
AVs1S = ['Vs1 Anisotropy (%) = ',num2str(AVs1,'%6.1f')];

if (plot_x==1||plot_x==3)
   ylabel(AVs1S,titleOpt{:})
end
if (plot_x==2||plot_x==4) 
   xlabel(AVs1S,titleOpt{:})
end

% mark maximum with black square and minimum with white circle
hold on
plot(Vs1_max_vector,blackMarker{:},'parent',mtexFig.gca)
plot(Vs1_min_vector,whiteMarker{:},'parent',mtexFig.gca)
hold off
%
%**************************************************************************
% Vs2 : Plot Vs2 velocities (km/s)
%**************************************************************************
%
% create a new axis
mtexFig.nextAxis
% plot
contourf(v_grid,vs2,'parent',mtexFig.gca);
mtexTitle('Vs2 (km/s)',titleOpt{:})
% plot Vs2 polarizations
hold on
quiver(p_grid,ps2,'linewidth',linewidth_pol,'color','black','parent',mtexFig.gca)
hold off

% Maximum and Minimum values and their directions
[Vs2_max_value,index_max] = max(vs2(:));
Vs2_max_vector = v_grid(index_max);

[Vs2_min_value,index_min] = min(vs2(:));
Vs2_min_vector = v_grid(index_min);

% Anisotropy percent
AVs2=200*(Vs2_max_value-Vs2_min_value)/(Vs2_max_value+Vs2_min_value);

% Percentage anisotropy
AVs2S = ['Vs2 Anisotropy (%) = ',num2str(AVs2,'%6.1f')];
if (plot_x==1||plot_x==3)
   ylabel(AVs2S,titleOpt{:})
end
if (plot_x==2||plot_x==4) 
   xlabel(AVs2S,titleOpt{:})
end

% mark maximum with black square and minimum with white circle
hold on
plot(Vs2_max_vector,blackMarker{:},'parent',mtexFig.gca)
plot(Vs2_min_vector,whiteMarker{:},'parent',mtexFig.gca)
hold off
%
%**************************************************************************
% dVs : Plot Velocity difference Vs1-Vs2 (km/s) plus Vs1 polarizations
%**************************************************************************
% calculate S-wave velocity difference
dVs = vs1-vs2;
% create a new axis
mtexFig.nextAxis
% plot
contourf(v_grid,dVs,'parent',mtexFig.gca);
mtexTitle('dVs=Vs1-Vs2 (km/s)',titleOpt{:})

% Maximum and Minimum values and thier directions
[dVs_max_value ,index_max] = max(dVs(:));
dVs_max_vector = v_grid(index_max);

[~,index_min] = min(dVs(:));
dVs_min_vector = v_grid(index_min);


% Max percentage anisotropy
AdVsS = ['Max dVs (km/s) = ',num2str(dVs_max_value,'%6.2f')];
if (plot_x==1||plot_x==3)
   ylabel(AdVsS,titleOpt{:})
end
if (plot_x==2||plot_x==4) 
   xlabel(AdVsS,titleOpt{:})
end

% mark maximum with black square and minimum with white circle
hold on
plot(dVs_max_vector,blackMarker{:},'parent',mtexFig.gca)
plot(dVs_min_vector,whiteMarker{:},'parent',mtexFig.gca)
hold off
%
%**************************************************************************
% Vp/Vs1 : Plot Vp/Vs1 ratio (no units)
%**************************************************************************
% calculate Vp/Vs1 velocity ratio
vpvs1 = vp./vs1;
% create a new axis
mtexFig.nextAxis
% plot
contourf(v_grid,vpvs1,'parent',mtexFig.gca);
mtexTitle('Vp/Vs1',titleOpt{:})

% Maximum and Minimum values and thier directions
[VpVs1_max_value,index_max] = max(vpvs1(:));
VpVs1_max_vector = v_grid(index_max);

[VpVs1_min_value,index_min] = min(vpvs1(:));
VpVs1_min_vector = v_grid(index_min);

% Anisotropy percent
AVpVs1=200*(VpVs1_max_value-VpVs1_min_value)/(VpVs1_max_value+VpVs1_min_value);

AVpVs1S = ['Vp/Vs1 Anisotropy (%) = ',num2str(AVpVs1,'%6.1f')];

if (plot_x==1||plot_x==3)
   ylabel(AVpVs1S,titleOpt{:})
end
if (plot_x==2||plot_x==4) 
   xlabel(AVpVs1S,titleOpt{:})
end
% mark maximum with black square and minimum with white circle
hold on
plot(VpVs1_max_vector,blackMarker{:},'parent',mtexFig.gca)
plot(VpVs1_min_vector,whiteMarker{:},'parent',mtexFig.gca)
hold off

%
%**************************************************************************
% Vp/Vs2 : Plot Vp/Vs2 ratio (no units)
%**************************************************************************
% calculate Vp/Vs2 velocity ratio
vpvs2 = vp./vs2;

% create a new axis
mtexFig.nextAxis
% plot
contourf(v_grid,vpvs2,'parent',mtexFig.gca);
mtexTitle('Vp/Vs2',titleOpt{:})

% Maximum and Minimum values and their directions
[VpVs2_max_value,index_max] = max(vpvs2(:));
VpVs2_max_vector = v_grid(index_max);

[VpVs2_min_value,index_min] = min(vpvs2(:));
VpVs2_min_vector = v_grid(index_min);

% Anisotropy percent
AVpVs2=200*(VpVs2_max_value-VpVs2_min_value)/(VpVs2_max_value+VpVs2_min_value);

AVpVs2S = ['Vp/Vs2 Anisotropy (%) = ',num2str(AVpVs2,'%6.1f')];
if (plot_x==1||plot_x==3)
   ylabel(AVpVs2S,titleOpt{:})
end
if (plot_x==2||plot_x==4) 
   xlabel(AVpVs2S,titleOpt{:})
end

% mark maximum with black square and minimum with white circle
hold on
plot(VpVs2_max_vector,blackMarker{:},'parent',mtexFig.gca)
plot(VpVs2_min_vector,whiteMarker{:},'parent',mtexFig.gca)
hold off
%
%
%**************************************************************************
% all colorbars and save file
%**************************************************************************
%
% add colorbars to all plots
colorbar
drawNow(gcm,'figSize','large')
%%
%
%
% save plot file
saveFigure(plot_name);
%
%
% reset colour map to default for MTEX 
setMTEXpref('defaultColorMap',WhiteJetColorMap);
%
disp(' ')
disp(' Plot has successfully terminated')
disp(' ')
%
% End of function
%
end
