Ice anisotropy and updated plot warning

examples/lowermost/lmm_anisotropy.m

@@ -47,6 +47,8 @@ function lmm_anisotropy()
 %  ** plot each of the derived properties on a globe
       plot_cprop(LON,LAT,Cani,[0 0.1],'Universal Elastic Anisotropy Index') ;
       plot_cprop(LON,LAT,PercH,[0 1.0],'Proportion of Hexagonality') ;
+      fprintf('Note: Refer to the accompanying file PoH.png rather than Walker\n') 
+      fprintf('      and Wookey (2012) Fig. 2b for the correct PoH plot.\n')
 
 end
 

msat/MS_ice.m

@@ -0,0 +1,168 @@
+% MS_ICE - Generate different ice fabrics
+%
+% // Part of MSAT - The Matlab Seismic Anisotropy Toolkit //
+%
+% Produce elastic tensor for different ice fabrics. 
+%
+% [C,rh] = MS_ice(...)
+%
+% Usage: 
+%     [ C,rh ] = MS_ice()                    
+%     [ C,rh ] = MS_ice('iso')
+%         Return elastic tensor for isotropic ice.
+%
+%     [ C,rh ] = MS_ice('cluster',phi)
+%         Return elastic tensor for an ice 'cluster' fabric, with a cluster
+%         opening angle of phi. The cluster is around the x3 axis.
+%
+%     [ C,rh ] = MS_ice('girdle',xi)
+%         Return elastic tensor for an ice 'girdle' fabric, with a cluster
+%         thickness angle of xi. The cluster is around the x1 axis.
+% 
+% Reference: Maurel et al. 2015, Proc. R. Soc. A 471:20140988. 
+%            doi:10.1098/rspa.2014.0988 
+
+% Copyright (c) 2011-2020, James Wookey and Andrew Walker
+% Copyright (c) 2007-2011, James Wookey
+% All rights reserved.
+% 
+% Redistribution and use in source and binary forms, 
+% with or without modification, are permitted provided 
+% that the following conditions are met:
+% 
+%    * Redistributions of source code must retain the 
+%      above copyright notice, this list of conditions 
+%      and the following disclaimer.
+%    * Redistributions in binary form must reproduce 
+%      the above copyright notice, this list of conditions 
+%      and the following disclaimer in the documentation 
+%      and/or other materials provided with the distribution.
+%    * Neither the name of the University of Bristol nor the names 
+%      of its contributors may be used to endorse or promote 
+%      products derived from this software without specific 
+%      prior written permission.
+% 
+% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 
+% AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED 
+% WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
+% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
+% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL 
+% THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY 
+% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
+% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF 
+% USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
+% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
+% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE 
+% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
+% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+function [C,rh] = MS_ice(varargin)
+
+   ice_mode = 0  ; % isotropic (default)
+          % = 1    % cluster
+          % = 2    % girdle
+
+   phi = NaN ;
+   xi  = NaN ;
+
+%  ** process the optional arguments
+   iarg = 1 ;
+   while iarg <= (length(varargin))
+      switch lower(varargin{iarg})
+         case 'iso'
+            ice_mode = 0 ;
+            iarg = iarg + 1 ;
+         case 'cluster'
+            ice_mode = 1 ;
+            phi = varargin{iarg+1} ;
+            iarg = iarg + 2 ;
+         case 'girdle'
+            ice_mode = 2 ;
+            xi = varargin{iarg+1} ;
+            iarg = iarg + 2 ;               
+         otherwise 
+            error(['Unknown option: ' varargin{iarg}]) ;   
+      end   
+   end
+
+   rh = 917.0 % same for all 
+
+   if ice_mode==0
+      C = MS_iso(3.871810,1.966181,rh) ;
+   elseif ice_mode==1
+      assert((phi<=90 & phi>=0), ...
+         'MS:ICE:BAD_PHI_VALUE','phi should be 0-90 degrees');
+      C = ice_cluster(phi) ;
+   elseif ice_mode==2
+      assert((phi<=90 & phi>=0), ...
+         'MS:ICE:BAD_XI_VALUE','xi should be 0-90 degrees');      
+      C = ice_girdle(xi) ;
+   else
+      error('MS:ICE:UNSUPPORTED_MODE','Unsupport fabric mode selected');
+   end
+
+   return
+end
+
+function [CC] = ice_girdle(xi)
+% Reference: Maurel et al. 2015, Proc. R. Soc. A 471:20140988. 
+%            doi:10.1098/rspa.2014.0988 
+% uses hardwired elastic coefficients
+
+% calculate the fabric using the equations
+A = 14.06 ;
+C = 15.24 ;
+L = 3.06  ;
+N = 3.455 ;
+F = 5.88  ;
+
+s=sind(xi) ;
+CC(1,1)=(1/15.) * (A*(15-10*s^2+3*s^4)+3*C*s^4+2*(2*L+F)*(5*s^2-3*s^4)) ;
+CC(2,2)=(1/120.) * (A*(45+10*s^2+9*s^4)+3*C*(15-10*s^2+3*s^4)+2*(2*L+F)*(15+10*s^2-9*s^4)) ;
+CC(4,4)=(1/120.) * ((A+C-2*F)*(15-10*s^2+3*s^4)+12*L*(5-s^4)+40*N*s^2) ;
+CC(5,5)=(1/30.) * ((A+C-2*F)*(5*s^2-3*s^4)+3*L*(5-5*s^2+4*s^4)+5*N*(3-s^2)) ;
+CC(1,2)=(1/30.) * (3*A*(5-s^4)+(C-4*L)*(5*s^2-3*s^4)-10*N*(3-s^2)+F*(15-5*s^2+6*s^4)) ;
+CC(1,3)=CC(1,2) ;
+CC(2,1)=CC(1,2) ;
+CC(2,3)=CC(2,2)-2*CC(4,4) ;
+CC(3,1)=CC(1,2) ;
+CC(3,2)=CC(2,2)-2*CC(4,4) ;
+CC(3,3)=CC(2,2) ;
+CC(6,6)=CC(5,5) ;
+
+end
+
+function [CC] = ice_cluster(phi)
+% Reference: Maurel et al. 2015, Proc. R. Soc. A 471:20140988. 
+%            doi:10.1098/rspa.2014.0988 
+% uses hardwired elastic coefficients
+
+% calculate the fabric using the equations
+A=14.06 ;
+C=15.24 ;
+L=3.06 ;
+N=3.455 ;
+F=5.88 ;
+
+X=1+cosd(phi)+cosd(phi)^2 ;
+Y=cosd(phi)^3+cosd(phi)^4 ;
+
+CC = zeros(6,6) ;
+CC(1,1) = (1/120.)*(A*(45+19*X+9*Y)+3*C*(15-7*X+3*Y)+2*(2*L+F)*(15+X-9*Y)) ;
+CC(3,3) = (1/15.)*(A*(15-7*X+3*Y)+3*C*(X+Y) + 2*(2*L+F)*(2*X-3*Y)) ;
+CC(4,4) = (1/30.)*((A+C-2*F)*(2*X-3*Y)+3*L*(5-X+4*Y)+5*N*(3-X)) ;
+CC(6,6) = (1/120.)*((A+C-2*F)*(15-7*X+3*Y)+12*L*(5-X-Y)+40*N*X) ;
+CC(1,3) = (1/30.)*(3*A*(5-X-Y)+(C-4*L)*(2*X-3*Y)-10*N*(3-X)+F*(15+X+6*Y)) ;
+CC(1,2) = CC(1,1) - 2.*CC(6,6) ;
+CC(2,1) = CC(1,1) - 2.*CC(6,6) ;
+CC(2,2) = CC(1,1) ;
+CC(2,3) = CC(1,3) ;
+CC(3,1) = CC(1,3) ;
+CC(3,2) = CC(1,3) ;
+CC(5,5) = CC(4,4) ;
+
+end 
+
+
+

msat/MS_phasevels.m

@@ -13,9 +13,9 @@
 %         an inclination and azimuth (both in degrees, see below). Output 
 %         details are given below.
 %
-%     [ pol, avs, vs1, vs2, vp, SF, SS ] = MS_phasevels( C, rh, inc, azi )                    
-%         Additionally output fast and slow S-wave polarisation in vector
-%         form.
+%     [ pol, avs, vs1, vs2, vp, SF, SS, PP ] = MS_phasevels( C, rh, inc, azi )                    
+%         Additionally output fast and slow S-wave, and P-wave, polarisation 
+%         in vector form.
 %
 % Notes:
 %     Azi is defined as the angle in degrees from the +ve 1-axis in x1-x2 
@@ -83,7 +83,7 @@
 % OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
 % SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-function [pol,avs,vs1,vs2,vp, S1P, S2P] = MS_phasevels(C,rh,inc,azi)
+function [pol,avs,vs1,vs2,vp, S1P, S2P, PP] = MS_phasevels(C,rh,inc,azi)
 
       if (length(inc)~=length(azi))
 			error('MS:ListsMustMatch', ...
@@ -112,6 +112,7 @@
 		S1 = zeros(length(azi),3) ;
 		S1P = zeros(length(azi),3) ;
       S2P = zeros(length(azi),3) ;
+      PP = zeros(length(azi),3) ;
 
 %   ** Handle isotropic case quickly
      if isIsotropic(C, isotol)
@@ -121,7 +122,8 @@
          vs2 = vs1;                     % Mbar to Pa.
          avs(:) = 0.0;
          pol(:) = NaN; % Both waves have same velocity... meaningless.
-         S1P(:) = NaN;
+         S1P(:) = NaN ;
+         PP(:) = NaN ;
          return
      end
 
@@ -138,6 +140,9 @@
 
 %  ** pull out the eigenvectors
       P  = EIGVEC(:,1) ;
+
+      PP(ipair,:) = P ;
+
       S1 = EIGVEC(:,2) ;
 
   		if ~isreal(S1)

msat/MS_sphere.m

@@ -12,7 +12,8 @@
 %     MS_sphere(CC,rh,mode)                    
 %         Mode can be 'p'(-wave), 's'(-wave), 's1' (fast s-wave) or 's2'
 %         (slow s-wave) for velocity plots or 'slowP','slowS1' or 'slowS2'
-%         for slowness plots (slowness is the inverse of velocity).
+%         for slowness plots (slowness is the inverse of velocity). 'pp'
+%         mode plots p-wave polarisation vectors.
 %
 %     MS_sphere(CC,rh,mode,...) 
 %          Further arguments are optional, can be combined in any order and
@@ -36,11 +37,11 @@
 %          coloured spots. 
 %
 %     MS_sphere(..., 'FSWTickLength', length)
-%          Set the length of the Fast shear wave direction markers, default
-%          is 0.08.
+%          Set the length of the shear wave direction / p-wave polarisations
+%          markers, default is 0.08.
 %
 %     MS_sphere(..., 'FSWMarkerSize', length)
-%          Set the size of the Fast shear wave direction markers, default
+%          Set the size of the direction markers, default
 %          is 4.
 %
 %     MS_sphere(..., 'cmap', CM)                    
@@ -231,7 +232,7 @@ function MS_sphere(CC,rh,mode,varargin)
 end
 
 if ischar(mode)
-   if sum(strcmpi({'p','s', 's1', 's2', 'slowp', 'slows1', 'slows2'},mode))~=1
+   if sum(strcmpi({'p','s', 's1', 's2', 'slowp', 'slows1', 'slows2','pp'},mode))~=1
       error('MS:SPHERE:badmode', ...
           'Mode must be ''S''(-wave), ''S1'', ''S2'', or ''P''(-wave) ')
    end
@@ -246,7 +247,7 @@ function MS_sphere(CC,rh,mode,varargin)
 end
 
 % check data overlay
-if iplot_pdata & ~strcmpi(mode,'p')
+if iplot_pdata & sum(strcmpi(mode,{'p','pp'}))==0
    error('MS:SPHERE:badpoverlay', 'P-wave data can only overlay a ''P'' mode plot.') ;
 end
 
@@ -275,7 +276,7 @@ function MS_sphere(CC,rh,mode,varargin)
 
 [~,avs,vs1,vs2,vp] = MS_phasevels(CC,rh,inc,az) ;
 
-if strcmpi(mode,'p')
+if sum(strcmpi(mode,{'p','pp'}))
    trisurf(faces,x,y,z,vp) ;
 elseif strcmpi(mode,'s1')
    trisurf(faces,x,y,z,vs1) ;
@@ -317,8 +318,8 @@ function MS_sphere(CC,rh,mode,varargin)
 
 if ~isnan(cax), caxis(cax), end
 
-% plot the shear-wave polarisation
-if sum(strcmpi({'s','s1','s2','slows1','slows2'},mode))==1 
+% plot the shear- or p-wave polarisation 
+if sum(strcmpi({'s','s1','s2','slows1','slows2','pp'},mode))==1 
    [x, y, z, ~, az, inc] =  get_mesh(polmesh);
 
 %  find the closest points to each specified direction   
@@ -335,7 +336,7 @@ function MS_sphere(CC,rh,mode,varargin)
       inc(ind) = din(idir)*180/pi ;
    end
 
-   [~,avs,vs1,vs2,vp, SF, SS] = MS_phasevels(CC,rh,inc,az) ;
+   [~,avs,vs1,vs2,vp, SF, SS, PP] = MS_phasevels(CC,rh,inc,az) ;
    % calculate PM vectors
    nv = length(x) ;
    if sum(strcmpi(mode,{'s','s1','slows1'}))==1 
@@ -355,6 +356,14 @@ function MS_sphere(CC,rh,mode,varargin)
           plot3(XI(1),XI(2),XI(3),'ko','MarkerSize',FSWMarkerSize,'MarkerFaceColor','k')
           plot3([X1(1) X2(1)],[X1(2) X2(2)],[X1(3) X2(3)],'k-','LineWidth',(FSWMarkerSize/2))
       end
+   elseif strcmpi(mode, 'pp')
+      for iv=1:nv
+          XI=[x(iv) y(iv) z(iv)] ;
+          X1 = XI - FSWTickLength.*PP(iv,:) ;
+          X2 = XI + FSWTickLength.*PP(iv,:) ;
+          plot3(XI(1),XI(2),XI(3),'ko','MarkerSize',FSWMarkerSize,'MarkerFaceColor','k')
+          plot3([X1(1) X2(1)],[X1(2) X2(2)],[X1(3) X2(3)],'k-','LineWidth',(FSWMarkerSize/2))
+      end
    else
       for iv=1:nv
 		  if strcmpi(mode, 'slows2')