diff --git a/examples/lowermost/PoH.png b/examples/lowermost/PoH.png new file mode 100644 index 0000000..10d2b10 Binary files /dev/null and b/examples/lowermost/PoH.png differ diff --git a/examples/lowermost/lmm_anisotropy.m b/examples/lowermost/lmm_anisotropy.m index 28e686c..6531878 100644 --- a/examples/lowermost/lmm_anisotropy.m +++ b/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 diff --git a/msat/MS_ice.m b/msat/MS_ice.m new file mode 100644 index 0000000..afdad75 --- /dev/null +++ b/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 + + + diff --git a/msat/MS_phasevels.m b/msat/MS_phasevels.m index 9d1d678..ed107c5 100644 --- a/msat/MS_phasevels.m +++ b/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) diff --git a/msat/MS_sphere.m b/msat/MS_sphere.m index d6dc1bb..41f4de7 100644 --- a/msat/MS_sphere.m +++ b/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')