From 803946b4f55e3ec58270c8a15006d27ae9ea205f Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Fri, 5 Feb 2016 13:43:50 +0000 Subject: [PATCH 1/7] added .DS_Stote to .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 6d725d9..43e41df 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ *.m~ +.DS_Store \ No newline at end of file From d57a98d9cfc616dfce8a6f0a165dd6ffa3ece319 Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Fri, 5 Feb 2016 16:10:31 +0000 Subject: [PATCH 2/7] added MS_groupvels (copied from MS_Phasevels with added function rayvel --- msat/MS_groupvels.m | 322 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 322 insertions(+) create mode 100644 msat/MS_groupvels.m diff --git a/msat/MS_groupvels.m b/msat/MS_groupvels.m new file mode 100644 index 0000000..62c7d44 --- /dev/null +++ b/msat/MS_groupvels.m @@ -0,0 +1,322 @@ +function [pol,avs,vs1,vs2,vp, S1P, S2P] = MS_groupvels(C,rh,inc,azi) + + if (length(inc)~=length(azi)) + error('MS:ListsMustMatch', ... + 'AZI and INC must be scalars or vectors of the same dimension'); + end + +% ** convert inc, azi to column vectors if necessary + inc = reshape(inc,length(inc),1) ; + azi = reshape(azi,length(azi),1) ; + + isotol = sqrt(eps); % Mbars + + % Check that C is valid (if check not suppressed) + MS_checkC(C); + + % ** convert GPa to MB file units (Mbars), density to g/cc + + C(:,:) = C(:,:) * 0.01 ; + rh = rh ./ 1e3 ; + + avs = zeros(size(azi)) ; + vp = zeros(size(azi)) ; + vs1 = zeros(size(azi)) ; + vs2 = zeros(size(azi)) ; + pol = zeros(size(azi)) ; + S1 = zeros(length(azi),3) ; + S1P = zeros(length(azi),3) ; + S2P = zeros(length(azi),3) ; + +% ** Handle isotropic case quickly + if isIsotropic(C, isotol) + vp(:) = sqrt(( ((1.0/3.0)*(C(1,1)+2*C(1,2)))+ ... + ((4.0/3.0)*C(4,4)) )/rh)*10.0; + vs1(:) = sqrt(C(4,4)/rh)*10.0; % Factor of 10 converts from + vs2 = vs1; % Mbar to Pa. + avs(:) = 0.0; + pol(:) = NaN; % Both waves have same velocity... meaningless. + S1P(:) = NaN; + return + end + +% ** start looping + for ipair = 1:length(inc) + cazi = azi(ipair) ; + cinc = inc(ipair) ; + +% ** create the cartesian vector + XI = cart2(cinc,cazi) ; + +% ** compute phase velocities + [V,EIGVEC]=velo(XI,rh,C) ; + +% ** pull out the eigenvectors + P = EIGVEC(:,1) ; + S1 = EIGVEC(:,2) ; + + if ~isreal(S1) + error_str = ['The S1 polarisation vector is not real!\n\n'... + sprintf('inc = %f, azi = %f\n\n',cinc,cazi) ... + sprintf('C = %f %f %f %f %f %f\n',C(1:6,1)) ... + sprintf(' %f %f %f %f %f %f\n',C(1:6,2)) ... + sprintf(' %f %f %f %f %f %f\n',C(1:6,3)) ... + sprintf(' %f %f %f %f %f %f\n',C(1:6,4)) ... + sprintf(' %f %f %f %f %f %f\n',C(1:6,5)) ... + sprintf(' %f %f %f %f %f %f\n\n',C(1:6,6)) ... + sprintf('S1 = %f %f %f\n',S1)]; + error('MS:PHASEVELS:vectornotreal', error_str) ; + end + S2 = EIGVEC(:,3) ; + +% ** calculate projection onto propagation plane + S1N = V_cross(XI,S1) ; + S1P(ipair,:) = V_cross(XI,S1N); + S2N = V_cross(XI,S2) ; + S2P(ipair,:) = V_cross(XI,S2N); + +% ** rotate into y-z plane to calculate angles +% (use functions optimised for the two needed +% rotations, see below). + [S1PR] = V_rot_gam(S1P(ipair,:),cazi) ; + [S1PRR] = V_rot_bet(S1PR,cinc) ; + + + + ph = atan2(S1PRR(2),S1PRR(3)) .* 180/pi ; + +% ** transform angle to between -90 and 90 + if (ph < -90.), ph = ph + 180.;end + if (ph > 90.), ph = ph - 180.;end + +% ** calculate some useful values + dVS = (V(2)-V(3)) ; + VSmean = (V(2)+V(3))/2.0 ; + + avs(ipair) = 100.0*(dVS/VSmean) ; + vp(ipair) = V(1) ; + vs1(ipair) = V(2) ; + vs2(ipair) = V(3) ; + + pol(ipair) = ph ; + end % ipair = 1:length(inc_in) + + % If any directions have zero avs (within machine accuracy) + % set pol to NaN - array wise: + isiso = real(avs > sqrt(eps)) ; % list of 1.0 and 0.0. + pol = pol .* (isiso./isiso) ; % times by 1.0 or NaN. + + S1P(:,1) = S1P(:,1) .* (isiso./isiso); + S1P(:,2) = S1P(:,2) .* (isiso./isiso); + S1P(:,3) = S1P(:,3) .* (isiso./isiso); + +return +%======================================================================================= + +function [c] = V_cross(a, b) + % This is ~10 times quicker than the Matlab cross() function + % for me (AMW). We assume the arguments are both 3-vectors and + % avoid the checks and reshaping needed for the more general case. + % (According to the prfiler, this moves cross from the most expensive + % child function costing ~50% of the time to the third most expensive + % child costing ~10% of the time). + + c = [a(2).*b(3)-a(3).*b(2) + a(3).*b(1)-a(1).*b(3) + a(1).*b(2)-a(2).*b(1)]; + +return + + +%======================================================================================= +% Rather useing a general case rotation about three angles we use two +% special case rotations around the c and b axes. This saves time as +% we don't need to convert 0 degrees to radians six times per call to +% phasevels, or do two double matrix multiplications. This saves ~60% +% of the time spent in vector rotation (and, with the cross thing above) +% reduces the 10000 evaulations needed for MS_anisotropy lma from ~19 +% secs to ~10 secs. + +function [VR] = V_rot_gam(V,gam) + +% Make rotation matrix +g = gam * pi/180. ; + +RR = [ cos(g) sin(g) 0 ; -sin(g) cos(g) 0 ; 0 0 1 ] ; +VR = V * RR ; + +return + +function [VR] = V_rot_bet(V,bet) + +% Make rotation matrix +b = bet * pi/180. ; + +RR = [ cos(b) 0 -sin(b) ; 0 1 0 ; sin(b) 0 cos(b) ] ; + +VR = V * RR ; + +return + +%======================================================================================= + +%======================================================================================= + function [X] = cart2(inc,azm) +%======================================================================================= +%c convert from spherical to cartesian co-ordinates +%c north x=100 west y=010 up z=001 +%c irev=+1 positive vector x +%c irev=-1 negative vector x +% NB: pre-converting azm and inc to radians and using +% cos and sin directly (instead to cosd and sind) +% is ~10x faster making the function ~4x faster. + azmr = azm.*(pi/180.0); + incr = inc.*(pi/180.0); + caz=cos(azmr) ; + saz=sin(azmr) ; + cinc=cos(incr) ; + sinc=sin(incr) ; + X=[caz*cinc -saz*cinc sinc] ; +%c normalise to direction cosines + r=sqrt(X(1)*X(1)+X(2)*X(2)+X(3)*X(3)) ; + + X = X./r ; + return +%======================================================================================= + +%======================================================================================= + function [V,EIGVEC]=velo(X,rh,C) +%======================================================================================= +% PHASE-VELOCITY SURFACES IN AN ANISOTROPIC MEDIUM +% revised April 1991 +% X(3) - DIRECTION OF INTEREST +% RHO - DENSITY +% V - PHASE VELOCITIES (1,2,3= P,S,SS) +% EIGVEC(3,3) - eigenvectors stored by columns +% +% Translated to MATLAB by James Wookey +% Revised, Andrew Walker 2012. + +% Form the 3x3 Christoffel tensor without converting the elasticity +% from 6x6 to 3x3x3x3 form using the formula from page 1076 of +% Winterstein 1990. This is > twice as fast as the quickest way I +% have found going via the full tensor form. + gamma = [X(1) 0.0 0.0 0.0 X(3) X(2) ; ... + 0.0 X(2) 0.0 X(3) 0.0 X(1) ; ... + 0.0 0.0 X(3) X(2) X(1) 0.0 ]; + T = gamma * C * gamma'; + +% determine the eigenvalues of symmetric tij + [EIVEC EIVAL] = eig(T) ; + +% calculate velocities and sort +% note that we could get a significant speedup if +% we could avoid the sort - LAPACK usually does sort +% in order of incresing eigenvectors but I've found +% cases where this does not happen (and we swap P and +% S wave velocities). Note that MATLAB does not "guarantee +% that the eignevalues are not returned in sorted order" +% http://www.mathworks.com/matlabcentral/newsreader/view_original/394371 +% so we have to sort here... + V_RAW = (sqrt([EIVAL(1,1) EIVAL(2,2) EIVAL(3,3)]./rh))*10. ; + [V IND] = sort(V_RAW,2,'descend') ; + EIGVEC = EIVEC ; % for dimensioning + for i=1:3 + EIGVEC(:,i) = EIVEC(:,IND(i)) ; + end + + return +%======================================================================================= + +function [ l ] = isIsotropic( C, tol ) + + % Are we isotropic - assume matrix is symmetrical at this point. + l = (abs(C(1,1)-C(2,2)) < tol) & (abs(C(1,1)-C(3,3)) < tol) & ... + (abs(C(1,2)-C(1,3)) < tol) & (abs(C(1,2)-C(2,3)) < tol) & ... + (abs(C(4,4)-C(5,5)) < tol) & (abs(C(4,4)-C(6,6)) < tol) & ... + (abs(C(1,4)) < tol) & (abs(C(1,5)) < tol) & (abs(C(1,6)) < tol) & ... + (abs(C(2,4)) < tol) & (abs(C(2,5)) < tol) & (abs(C(2,6)) < tol) & ... + (abs(C(3,4)) < tol) & (abs(C(3,5)) < tol) & (abs(C(3,6)) < tol) & ... + (abs(C(4,5)) < tol) & (abs(C(4,6)) < tol) & (abs(C(5,6)) < tol) & ... + (((C(1,1)-C(1,2))/2.0)-C(4,4) < tol); + +return + + + + + + +function VG =rayvel(C,SN,rho) +% TO CALCULATE THE RAY-VELOCITY VECTOR CORRESPONDING TO A SLOWNESS VECTOR. +% Original fortran by David Mainprice as part of the EMATRIX code. +% Converted to Python by Alan Baird +% +% C: Stiffness tensor in Voigt Notation (6X6). +% SN: Slowness vector (3). +% rho: Density +% +% returns VG: Group velocity vector (3) + +ijkl = [1,6,5; ... + 6,2,4; ... + 5,4,3] ; + + +gamma = [X(1) 0.0 0.0 0.0 X(3) X(2) ; ... + 0.0 X(2) 0.0 X(3) 0.0 X(1) ; ... + 0.0 0.0 X(3) X(2) X(1) 0.0 ]; + +F = gamma * C * gamma'-eye(3); + + +% Signed cofactors of F[i,k] +CF = zeros(3,3) + +CF(1,1)=F(2,2)*F(3,3)-F(2,3)**2 +CF(2,2)=F(1,1)*F(3,3)-F(1,3)**2 +CF(3,3)=F(1,1)*F(2,2)-F(1,2)**2 +CF(1,2)=F(2,3)*F(3,1)-F(2,1)*F(3,3) +CF(2,1)=CF(1,2) +CF(2,3)=F(3,1)*F(1,2)-F(3,2)*F(1,1) +CF(3,2)=CF(2,3) +CF(3,1)=F(1,2)*F(2,3)-F(1,3)*F(2,2) +CF(1,3)=CF(3,1) + + +% Derivatives of determinant elements +DF = zeros(3,3,3) +for i=1:3 + for j=1:3 + for k=1:3 + DF(i,j,k)=0.0 + for l=1:3 + DF(i,j,k) = DF(i,j,k) + (C(ijkl(i,j),ijkl(k,l))+ C(ijkl(k,j),ijkl(i,l)) ) * SN(l) + end + end + end +end + +% Components of Gradient +DFD = zeros(3) +for k=1:3 + DFD(k) = 0.0 + for i=1:3 + for j=1:3 + DFD(k)=DFD(k)+DF(i,j,k)+CF(i,j) + end + end +end + +% Normalize to obtain group velocity +denom = 0.0 +VG = zeros(3) +for i=1:3 + denom = denom+SN(i)*DFD(i) +end +for i=1:3 + VG(i) = DFD(i)/denom +end + +end % function + From cae5a846653c001d5ad4855d521ec5c2f016eeb8 Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Fri, 5 Feb 2016 17:39:55 +0000 Subject: [PATCH 3/7] MS_groupvels should now output group velocities --- msat/MS_groupvels.m | 107 +++++++++++++------------------------------- 1 file changed, 31 insertions(+), 76 deletions(-) diff --git a/msat/MS_groupvels.m b/msat/MS_groupvels.m index 62c7d44..05796dd 100644 --- a/msat/MS_groupvels.m +++ b/msat/MS_groupvels.m @@ -1,4 +1,4 @@ -function [pol,avs,vs1,vs2,vp, S1P, S2P] = MS_groupvels(C,rh,inc,azi) +function [PP ,S1P ,S2P ,SNP ,SNS1,SNS2,VGP ,VGS1,VGS2] = MS_groupvels(C,rh,inc,azi) if (length(inc)~=length(azi)) error('MS:ListsMustMatch', ... @@ -9,36 +9,26 @@ inc = reshape(inc,length(inc),1) ; azi = reshape(azi,length(azi),1) ; - isotol = sqrt(eps); % Mbars - - % Check that C is valid (if check not suppressed) - MS_checkC(C); + % ** convert GPa to MB file units (Mbars), density to g/cc C(:,:) = C(:,:) * 0.01 ; rh = rh ./ 1e3 ; - avs = zeros(size(azi)) ; - vp = zeros(size(azi)) ; - vs1 = zeros(size(azi)) ; - vs2 = zeros(size(azi)) ; - pol = zeros(size(azi)) ; - S1 = zeros(length(azi),3) ; - S1P = zeros(length(azi),3) ; - S2P = zeros(length(azi),3) ; + + + + PP = zeros(length(azi),3) ; + S1P = zeros(length(azi),3) ; + S2P = zeros(length(azi),3) ; + SNP = zeros(length(azi),3) ; + SNS1 = zeros(length(azi),3) ; + SNS2 = zeros(length(azi),3) ; + VGP = zeros(length(azi),3) ; + VGS1 = zeros(length(azi),3) ; + VGS2 = zeros(length(azi),3) ; -% ** Handle isotropic case quickly - if isIsotropic(C, isotol) - vp(:) = sqrt(( ((1.0/3.0)*(C(1,1)+2*C(1,2)))+ ... - ((4.0/3.0)*C(4,4)) )/rh)*10.0; - vs1(:) = sqrt(C(4,4)/rh)*10.0; % Factor of 10 converts from - vs2 = vs1; % Mbar to Pa. - avs(:) = 0.0; - pol(:) = NaN; % Both waves have same velocity... meaningless. - S1P(:) = NaN; - return - end % ** start looping for ipair = 1:length(inc) @@ -52,63 +42,25 @@ [V,EIGVEC]=velo(XI,rh,C) ; % ** pull out the eigenvectors - P = EIGVEC(:,1) ; - S1 = EIGVEC(:,2) ; - - if ~isreal(S1) - error_str = ['The S1 polarisation vector is not real!\n\n'... - sprintf('inc = %f, azi = %f\n\n',cinc,cazi) ... - sprintf('C = %f %f %f %f %f %f\n',C(1:6,1)) ... - sprintf(' %f %f %f %f %f %f\n',C(1:6,2)) ... - sprintf(' %f %f %f %f %f %f\n',C(1:6,3)) ... - sprintf(' %f %f %f %f %f %f\n',C(1:6,4)) ... - sprintf(' %f %f %f %f %f %f\n',C(1:6,5)) ... - sprintf(' %f %f %f %f %f %f\n\n',C(1:6,6)) ... - sprintf('S1 = %f %f %f\n',S1)]; - error('MS:PHASEVELS:vectornotreal', error_str) ; - end - S2 = EIGVEC(:,3) ; - -% ** calculate projection onto propagation plane - S1N = V_cross(XI,S1) ; - S1P(ipair,:) = V_cross(XI,S1N); - S2N = V_cross(XI,S2) ; - S2P(ipair,:) = V_cross(XI,S2N); - -% ** rotate into y-z plane to calculate angles -% (use functions optimised for the two needed -% rotations, see below). - [S1PR] = V_rot_gam(S1P(ipair,:),cazi) ; - [S1PRR] = V_rot_bet(S1PR,cinc) ; - - - - ph = atan2(S1PRR(2),S1PRR(3)) .* 180/pi ; - -% ** transform angle to between -90 and 90 - if (ph < -90.), ph = ph + 180.;end - if (ph > 90.), ph = ph - 180.;end + PP(ipair,:) = EIGVEC(:,1) ; + S1P(ipair,:) = EIGVEC(:,2) ; + S2P(ipair,:) = EIGVEC(:,3) ; + +% ** slowness vectors + SNP(ipair,:) = XI/V(1) + SNS1(ipair,:) = XI/V(2) + SNS2(ipair,:) = XI/V(3) + +% ** Group velocity vectors (need to convert density back first) + VGP(ipair,:) = rayvel(C,SNP(ipair,:),rh*1e3) + VGS1(ipair,:) = rayvel(C,SNS1(ipair,:),rh*1e3) + VGS2(ipair,:) = rayvel(C,SNS2(ipair,:),rh*1e3) -% ** calculate some useful values - dVS = (V(2)-V(3)) ; - VSmean = (V(2)+V(3))/2.0 ; + - avs(ipair) = 100.0*(dVS/VSmean) ; - vp(ipair) = V(1) ; - vs1(ipair) = V(2) ; - vs2(ipair) = V(3) ; - - pol(ipair) = ph ; end % ipair = 1:length(inc_in) - % If any directions have zero avs (within machine accuracy) - % set pol to NaN - array wise: - isiso = real(avs > sqrt(eps)) ; % list of 1.0 and 0.0. - pol = pol .* (isiso./isiso) ; % times by 1.0 or NaN. - S1P(:,1) = S1P(:,1) .* (isiso./isiso); - S1P(:,2) = S1P(:,2) .* (isiso./isiso); - S1P(:,3) = S1P(:,3) .* (isiso./isiso); return %======================================================================================= @@ -320,3 +272,6 @@ end % function + + + From 349301c5775f639e4adc8a8c70f2f3f769aa4b85 Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Fri, 4 Mar 2016 17:37:04 +0000 Subject: [PATCH 4/7] updated phasevels to output eigenvectors and prop vector, used it in groupvels --- msat/MS_groupvels.m | 315 +++++++++++++++++--------------------------- msat/MS_phasevels.m | 62 +++++++-- 2 files changed, 174 insertions(+), 203 deletions(-) diff --git a/msat/MS_groupvels.m b/msat/MS_groupvels.m index 05796dd..13de66e 100644 --- a/msat/MS_groupvels.m +++ b/msat/MS_groupvels.m @@ -1,27 +1,64 @@ -function [PP ,S1P ,S2P ,SNP ,SNS1,SNS2,VGP ,VGS1,VGS2] = MS_groupvels(C,rh,inc,azi) - - if (length(inc)~=length(azi)) - error('MS:ListsMustMatch', ... - 'AZI and INC must be scalars or vectors of the same dimension'); - end - -% ** convert inc, azi to column vectors if necessary - inc = reshape(inc,length(inc),1) ; - azi = reshape(azi,length(azi),1) ; - - - - % ** convert GPa to MB file units (Mbars), density to g/cc +% MS_PHASEVELS - Wave velocities in anisotropic media. +% +% // Part of MSAT - The Matlab Seismic Anisotropy Toolkit // +% +% Calculate the group velocity details for an elsticity matrix. +% +% [ VGP, VGS1, VGS2, ...] = MS_groupvels( C, rh, inc, azi ) +% +% Usage: +% [ VGP, VGS1, VGS2, ...] = MS_groupvels( C, rh, inc, azi ) +% Calculate group velocity vectors from elasticity matrix C (in GPa) and +% density rh (in kg/m^3) corresponding to a phase angle defined by +% an inclination and azimuth (both in degrees). Output +% details are given below. +% +% [ VGP, VGS1, VGS2, PE, S1E, S2E ] = MS_groupvels( C, rh, inc, azi ) +% Additionally output P, S1 and S2-wave polarisations in vector +% form. +% +% [ VGP, VGS1, VGS2, PE, S1E, S2E, SNP, SNS1, SNS2 ] = MS_groupvels( C, rh, inc, azi ) +% Additionally output P, S1 and S2 Slowness vectors +% - C(:,:) = C(:,:) * 0.01 ; - rh = rh ./ 1e3 ; - - - - - PP = zeros(length(azi),3) ; - S1P = zeros(length(azi),3) ; - S2P = zeros(length(azi),3) ; +% Copyright (c) 2016, Alan Baird +% 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 [ varargout ] = MS_groupvels(C,rh,inc,azi) + + [~,~,vs1,vs2,vp,~,~,PE,S1E,S2E,XIS] = MS_phasevels(C,rh,inc,azi); + SNP = zeros(length(azi),3) ; SNS1 = zeros(length(azi),3) ; SNS2 = zeros(length(azi),3) ; @@ -32,167 +69,55 @@ % ** start looping for ipair = 1:length(inc) - cazi = azi(ipair) ; - cinc = inc(ipair) ; - -% ** create the cartesian vector - XI = cart2(cinc,cazi) ; -% ** compute phase velocities - [V,EIGVEC]=velo(XI,rh,C) ; - -% ** pull out the eigenvectors - PP(ipair,:) = EIGVEC(:,1) ; - S1P(ipair,:) = EIGVEC(:,2) ; - S2P(ipair,:) = EIGVEC(:,3) ; - % ** slowness vectors - SNP(ipair,:) = XI/V(1) - SNS1(ipair,:) = XI/V(2) - SNS2(ipair,:) = XI/V(3) + SNP(ipair,:) = XIS(ipair,:)./vp(ipair); + SNS1(ipair,:) = XIS(ipair,:)./vs1(ipair); + SNS2(ipair,:) = XIS(ipair,:)./vs2(ipair); % ** Group velocity vectors (need to convert density back first) - VGP(ipair,:) = rayvel(C,SNP(ipair,:),rh*1e3) - VGS1(ipair,:) = rayvel(C,SNS1(ipair,:),rh*1e3) - VGS2(ipair,:) = rayvel(C,SNS2(ipair,:),rh*1e3) - - - - end % ipair = 1:length(inc_in) - - - -return -%======================================================================================= -function [c] = V_cross(a, b) - % This is ~10 times quicker than the Matlab cross() function - % for me (AMW). We assume the arguments are both 3-vectors and - % avoid the checks and reshaping needed for the more general case. - % (According to the prfiler, this moves cross from the most expensive - % child function costing ~50% of the time to the third most expensive - % child costing ~10% of the time). + VGP(ipair,:) = rayvel(C,SNP(ipair,:),rh); + VGS1(ipair,:) = rayvel(C,SNS1(ipair,:),rh); + VGS2(ipair,:) = rayvel(C,SNS2(ipair,:),rh); - c = [a(2).*b(3)-a(3).*b(2) - a(3).*b(1)-a(1).*b(3) - a(1).*b(2)-a(2).*b(1)]; + end -return - - -%======================================================================================= -% Rather useing a general case rotation about three angles we use two -% special case rotations around the c and b axes. This saves time as -% we don't need to convert 0 degrees to radians six times per call to -% phasevels, or do two double matrix multiplications. This saves ~60% -% of the time spent in vector rotation (and, with the cross thing above) -% reduces the 10000 evaulations needed for MS_anisotropy lma from ~19 -% secs to ~10 secs. + switch nargout + case 3 + varargout{1} = VGP ; + varargout{2} = VGS1 ; + varargout{3} = VGS2 ; -function [VR] = V_rot_gam(V,gam) -% Make rotation matrix -g = gam * pi/180. ; + case 6 + varargout{1} = VGP ; + varargout{2} = VGS1 ; + varargout{3} = VGS2 ; + varargout{4} = PE ; + varargout{5} = S1E ; + varargout{6} = S2E ; -RR = [ cos(g) sin(g) 0 ; -sin(g) cos(g) 0 ; 0 0 1 ] ; -VR = V * RR ; - -return - -function [VR] = V_rot_bet(V,bet) - -% Make rotation matrix -b = bet * pi/180. ; - -RR = [ cos(b) 0 -sin(b) ; 0 1 0 ; sin(b) 0 cos(b) ] ; - -VR = V * RR ; - -return - -%======================================================================================= - -%======================================================================================= - function [X] = cart2(inc,azm) -%======================================================================================= -%c convert from spherical to cartesian co-ordinates -%c north x=100 west y=010 up z=001 -%c irev=+1 positive vector x -%c irev=-1 negative vector x -% NB: pre-converting azm and inc to radians and using -% cos and sin directly (instead to cosd and sind) -% is ~10x faster making the function ~4x faster. - azmr = azm.*(pi/180.0); - incr = inc.*(pi/180.0); - caz=cos(azmr) ; - saz=sin(azmr) ; - cinc=cos(incr) ; - sinc=sin(incr) ; - X=[caz*cinc -saz*cinc sinc] ; -%c normalise to direction cosines - r=sqrt(X(1)*X(1)+X(2)*X(2)+X(3)*X(3)) ; - - X = X./r ; - return -%======================================================================================= - -%======================================================================================= - function [V,EIGVEC]=velo(X,rh,C) -%======================================================================================= -% PHASE-VELOCITY SURFACES IN AN ANISOTROPIC MEDIUM -% revised April 1991 -% X(3) - DIRECTION OF INTEREST -% RHO - DENSITY -% V - PHASE VELOCITIES (1,2,3= P,S,SS) -% EIGVEC(3,3) - eigenvectors stored by columns -% -% Translated to MATLAB by James Wookey -% Revised, Andrew Walker 2012. - -% Form the 3x3 Christoffel tensor without converting the elasticity -% from 6x6 to 3x3x3x3 form using the formula from page 1076 of -% Winterstein 1990. This is > twice as fast as the quickest way I -% have found going via the full tensor form. - gamma = [X(1) 0.0 0.0 0.0 X(3) X(2) ; ... - 0.0 X(2) 0.0 X(3) 0.0 X(1) ; ... - 0.0 0.0 X(3) X(2) X(1) 0.0 ]; - T = gamma * C * gamma'; - -% determine the eigenvalues of symmetric tij - [EIVEC EIVAL] = eig(T) ; - -% calculate velocities and sort -% note that we could get a significant speedup if -% we could avoid the sort - LAPACK usually does sort -% in order of incresing eigenvectors but I've found -% cases where this does not happen (and we swap P and -% S wave velocities). Note that MATLAB does not "guarantee -% that the eignevalues are not returned in sorted order" -% http://www.mathworks.com/matlabcentral/newsreader/view_original/394371 -% so we have to sort here... - V_RAW = (sqrt([EIVAL(1,1) EIVAL(2,2) EIVAL(3,3)]./rh))*10. ; - [V IND] = sort(V_RAW,2,'descend') ; - EIGVEC = EIVEC ; % for dimensioning - for i=1:3 - EIGVEC(:,i) = EIVEC(:,IND(i)) ; - end + + case 9 + varargout{1} = VGP ; + varargout{2} = VGS1 ; + varargout{3} = VGS2 ; + varargout{4} = PE ; + varargout{5} = S1E ; + varargout{6} = S2E ; + varargout{7} = SNP ; + varargout{8} = SNS1 ; + varargout{9} = SNS2 ; + + otherwise + error('MS:GROUPVELS:BadOutputArgs','Requires 3, 6, or 9 output arguments.') + end - return -%======================================================================================= -function [ l ] = isIsotropic( C, tol ) - % Are we isotropic - assume matrix is symmetrical at this point. - l = (abs(C(1,1)-C(2,2)) < tol) & (abs(C(1,1)-C(3,3)) < tol) & ... - (abs(C(1,2)-C(1,3)) < tol) & (abs(C(1,2)-C(2,3)) < tol) & ... - (abs(C(4,4)-C(5,5)) < tol) & (abs(C(4,4)-C(6,6)) < tol) & ... - (abs(C(1,4)) < tol) & (abs(C(1,5)) < tol) & (abs(C(1,6)) < tol) & ... - (abs(C(2,4)) < tol) & (abs(C(2,5)) < tol) & (abs(C(2,6)) < tol) & ... - (abs(C(3,4)) < tol) & (abs(C(3,5)) < tol) & (abs(C(3,6)) < tol) & ... - (abs(C(4,5)) < tol) & (abs(C(4,6)) < tol) & (abs(C(5,6)) < tol) & ... - (((C(1,1)-C(1,2))/2.0)-C(4,4) < tol); - return +%======================================================================================= @@ -215,62 +140,62 @@ 5,4,3] ; -gamma = [X(1) 0.0 0.0 0.0 X(3) X(2) ; ... - 0.0 X(2) 0.0 X(3) 0.0 X(1) ; ... - 0.0 0.0 X(3) X(2) X(1) 0.0 ]; +gamma = [SN(1) 0.0 0.0 0.0 SN(3) SN(2) ; ... + 0.0 SN(2) 0.0 SN(3) 0.0 SN(1) ; ... + 0.0 0.0 SN(3) SN(2) SN(1) 0.0 ]; -F = gamma * C * gamma'-eye(3); +F = gamma * C * gamma'-eye(3).*rho; % Signed cofactors of F[i,k] -CF = zeros(3,3) +CF = zeros(3,3); -CF(1,1)=F(2,2)*F(3,3)-F(2,3)**2 -CF(2,2)=F(1,1)*F(3,3)-F(1,3)**2 -CF(3,3)=F(1,1)*F(2,2)-F(1,2)**2 -CF(1,2)=F(2,3)*F(3,1)-F(2,1)*F(3,3) -CF(2,1)=CF(1,2) -CF(2,3)=F(3,1)*F(1,2)-F(3,2)*F(1,1) -CF(3,2)=CF(2,3) -CF(3,1)=F(1,2)*F(2,3)-F(1,3)*F(2,2) -CF(1,3)=CF(3,1) +CF(1,1)=F(2,2)*F(3,3)-F(2,3)^2; +CF(2,2)=F(1,1)*F(3,3)-F(1,3)^2; +CF(3,3)=F(1,1)*F(2,2)-F(1,2)^2; +CF(1,2)=F(2,3)*F(3,1)-F(2,1)*F(3,3); +CF(2,1)=CF(1,2); +CF(2,3)=F(3,1)*F(1,2)-F(3,2)*F(1,1); +CF(3,2)=CF(2,3); +CF(3,1)=F(1,2)*F(2,3)-F(1,3)*F(2,2); +CF(1,3)=CF(3,1); % Derivatives of determinant elements -DF = zeros(3,3,3) +DF = zeros(3,3,3); for i=1:3 for j=1:3 for k=1:3 - DF(i,j,k)=0.0 + DF(i,j,k)=0.0; for l=1:3 - DF(i,j,k) = DF(i,j,k) + (C(ijkl(i,j),ijkl(k,l))+ C(ijkl(k,j),ijkl(i,l)) ) * SN(l) + DF(i,j,k) = DF(i,j,k) + (C(ijkl(i,j),ijkl(k,l))+ C(ijkl(k,j),ijkl(i,l)) ) * SN(l); end end end end % Components of Gradient -DFD = zeros(3) +DFD = zeros(3,1); for k=1:3 - DFD(k) = 0.0 + DFD(k) = 0.0; for i=1:3 for j=1:3 - DFD(k)=DFD(k)+DF(i,j,k)+CF(i,j) + DFD(k)=DFD(k)+DF(i,j,k)+CF(i,j); end end end % Normalize to obtain group velocity -denom = 0.0 -VG = zeros(3) +denom = 0.0; +VG = zeros(3,1); for i=1:3 - denom = denom+SN(i)*DFD(i) + denom = denom+SN(i)*DFD(i); end for i=1:3 - VG(i) = DFD(i)/denom + VG(i) = DFD(i)./denom; end -end % function +return % function diff --git a/msat/MS_phasevels.m b/msat/MS_phasevels.m index 9d1d678..f133bb9 100644 --- a/msat/MS_phasevels.m +++ b/msat/MS_phasevels.m @@ -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,PE,S1E,S2E,XIS] = MS_phasevels(C,rh,inc,azi) if (length(inc)~=length(azi)) error('MS:ListsMustMatch', ... @@ -104,14 +104,20 @@ C(:,:) = C(:,:) * 0.01 ; rh = rh ./ 1e3 ; - avs = zeros(size(azi)) ; - vp = zeros(size(azi)) ; - vs1 = zeros(size(azi)) ; - vs2 = zeros(size(azi)) ; - pol = zeros(size(azi)) ; - S1 = zeros(length(azi),3) ; - S1P = zeros(length(azi),3) ; + avs = zeros(size(azi)) ; + vp = zeros(size(azi)) ; + vs1 = zeros(size(azi)) ; + vs2 = zeros(size(azi)) ; + pol = zeros(size(azi)) ; + S1P = zeros(length(azi),3) ; S2P = zeros(length(azi),3) ; + + % Eigenvectors + PE = zeros(length(azi),3) ; + S1E = zeros(length(azi),3) ; + S2E = zeros(length(azi),3) ; + % Cartesion propagation vectors + XIS = zeros(length(azi),3) ; % ** Handle isotropic case quickly if isIsotropic(C, isotol) @@ -132,6 +138,7 @@ % ** create the cartesian vector XI = cart2(cinc,cazi) ; + XIS(ipair,:)=XI ; % ** compute phase velocities [V,EIGVEC]=velo(XI,rh,C) ; @@ -153,6 +160,10 @@ error('MS:PHASEVELS:vectornotreal', error_str) ; end S2 = EIGVEC(:,3) ; + + PE(ipair,:) = P ; + S1E(ipair,:) = S1 ; + S2E(ipair,:) = S2 ; % ** calculate projection onto propagation plane S1N = V_cross(XI,S1) ; @@ -194,6 +205,41 @@ S1P(:,1) = S1P(:,1) .* (isiso./isiso); S1P(:,2) = S1P(:,2) .* (isiso./isiso); S1P(:,3) = S1P(:,3) .* (isiso./isiso); + + % switch nargout + % case 5 + % varargout{1} = pol ; + % varargout{2} = avs ; + % varargout{3} = vs1 ; + % varargout{4} = vs2 ; + % varargout{5} = vp ; + % + % case 7 + % varargout{1} = pol ; + % varargout{2} = avs ; + % varargout{3} = vs1 ; + % varargout{4} = vs2 ; + % varargout{5} = vp ; + % varargout{6} = S1P ; + % varargout{7} = S2P ; + % + % case 11 + % varargout{1} = pol ; + % varargout{2} = avs ; + % varargout{3} = vs1 ; + % varargout{4} = vs2 ; + % varargout{5} = vp ; + % varargout{6} = S1P ; + % varargout{7} = S2P ; + % varargout{8} = PE ; + % varargout{9} = S1E ; + % varargout{10} = S2E ; + % varargout{11} = XIS ; + % + % otherwise + % error('MS:PHASEVELS:BadOutputArgs','Requires 5, 7, or 11 output arguments.') + % end + return %======================================================================================= From 6270c9883d57e2477d7d3d229147a5789d85909b Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Thu, 10 Mar 2016 15:16:46 +0000 Subject: [PATCH 5/7] fixed some bugs in MS_groupvels --- msat/MS_groupvels.m | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/msat/MS_groupvels.m b/msat/MS_groupvels.m index 13de66e..292b73b 100644 --- a/msat/MS_groupvels.m +++ b/msat/MS_groupvels.m @@ -65,6 +65,8 @@ VGP = zeros(length(azi),3) ; VGS1 = zeros(length(azi),3) ; VGS2 = zeros(length(azi),3) ; + + % ** start looping @@ -75,11 +77,11 @@ SNS1(ipair,:) = XIS(ipair,:)./vs1(ipair); SNS2(ipair,:) = XIS(ipair,:)./vs2(ipair); -% ** Group velocity vectors (need to convert density back first) +% ** Group velocity vectors - VGP(ipair,:) = rayvel(C,SNP(ipair,:),rh); - VGS1(ipair,:) = rayvel(C,SNS1(ipair,:),rh); - VGS2(ipair,:) = rayvel(C,SNS2(ipair,:),rh); + VGP(ipair,:) = rayvel(C,SNP(ipair,:),rh./1e3); + VGS1(ipair,:) = rayvel(C,SNS1(ipair,:),rh./1e3); + VGS2(ipair,:) = rayvel(C,SNS2(ipair,:),rh./1e3); end @@ -180,7 +182,7 @@ DFD(k) = 0.0; for i=1:3 for j=1:3 - DFD(k)=DFD(k)+DF(i,j,k)+CF(i,j); + DFD(k)=DFD(k)+DF(i,j,k)*CF(i,j); end end end From e289626e48ac9675a2426be3a3f94b327eceb987 Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Fri, 11 Mar 2016 11:56:50 +0000 Subject: [PATCH 6/7] tidied up indentation --- msat/MS_groupvels.m | 98 +++++++++--------- msat/MS_phasevels.m | 239 +++++++++++++++++++------------------------- 2 files changed, 151 insertions(+), 186 deletions(-) diff --git a/msat/MS_groupvels.m b/msat/MS_groupvels.m index 292b73b..4e61d5e 100644 --- a/msat/MS_groupvels.m +++ b/msat/MS_groupvels.m @@ -137,65 +137,65 @@ % % returns VG: Group velocity vector (3) -ijkl = [1,6,5; ... - 6,2,4; ... - 5,4,3] ; - - -gamma = [SN(1) 0.0 0.0 0.0 SN(3) SN(2) ; ... - 0.0 SN(2) 0.0 SN(3) 0.0 SN(1) ; ... - 0.0 0.0 SN(3) SN(2) SN(1) 0.0 ]; - -F = gamma * C * gamma'-eye(3).*rho; - - + ijkl = [1,6,5; ... + 6,2,4; ... + 5,4,3] ; + + + gamma = [SN(1) 0.0 0.0 0.0 SN(3) SN(2) ; ... + 0.0 SN(2) 0.0 SN(3) 0.0 SN(1) ; ... + 0.0 0.0 SN(3) SN(2) SN(1) 0.0 ]; + + F = gamma * C * gamma'-eye(3).*rho; + + % Signed cofactors of F[i,k] -CF = zeros(3,3); - -CF(1,1)=F(2,2)*F(3,3)-F(2,3)^2; -CF(2,2)=F(1,1)*F(3,3)-F(1,3)^2; -CF(3,3)=F(1,1)*F(2,2)-F(1,2)^2; -CF(1,2)=F(2,3)*F(3,1)-F(2,1)*F(3,3); -CF(2,1)=CF(1,2); -CF(2,3)=F(3,1)*F(1,2)-F(3,2)*F(1,1); -CF(3,2)=CF(2,3); -CF(3,1)=F(1,2)*F(2,3)-F(1,3)*F(2,2); -CF(1,3)=CF(3,1); - - + CF = zeros(3,3); + + CF(1,1)=F(2,2)*F(3,3)-F(2,3)^2; + CF(2,2)=F(1,1)*F(3,3)-F(1,3)^2; + CF(3,3)=F(1,1)*F(2,2)-F(1,2)^2; + CF(1,2)=F(2,3)*F(3,1)-F(2,1)*F(3,3); + CF(2,1)=CF(1,2); + CF(2,3)=F(3,1)*F(1,2)-F(3,2)*F(1,1); + CF(3,2)=CF(2,3); + CF(3,1)=F(1,2)*F(2,3)-F(1,3)*F(2,2); + CF(1,3)=CF(3,1); + + % Derivatives of determinant elements -DF = zeros(3,3,3); -for i=1:3 - for j=1:3 - for k=1:3 - DF(i,j,k)=0.0; - for l=1:3 - DF(i,j,k) = DF(i,j,k) + (C(ijkl(i,j),ijkl(k,l))+ C(ijkl(k,j),ijkl(i,l)) ) * SN(l); + DF = zeros(3,3,3); + for i=1:3 + for j=1:3 + for k=1:3 + DF(i,j,k)=0.0; + for l=1:3 + DF(i,j,k) = DF(i,j,k) + (C(ijkl(i,j),ijkl(k,l))+ C(ijkl(k,j),ijkl(i,l)) ) * SN(l); + end end end end -end % Components of Gradient -DFD = zeros(3,1); -for k=1:3 - DFD(k) = 0.0; - for i=1:3 - for j=1:3 - DFD(k)=DFD(k)+DF(i,j,k)*CF(i,j); - end + DFD = zeros(3,1); + for k=1:3 + DFD(k) = 0.0; + for i=1:3 + for j=1:3 + DFD(k)=DFD(k)+DF(i,j,k)*CF(i,j); + end + end end -end % Normalize to obtain group velocity -denom = 0.0; -VG = zeros(3,1); -for i=1:3 - denom = denom+SN(i)*DFD(i); -end -for i=1:3 - VG(i) = DFD(i)./denom; -end + denom = 0.0; + VG = zeros(3,1); + for i=1:3 + denom = denom+SN(i)*DFD(i); + end + for i=1:3 + VG(i) = DFD(i)./denom; + end return % function diff --git a/msat/MS_phasevels.m b/msat/MS_phasevels.m index f133bb9..2e00b81 100644 --- a/msat/MS_phasevels.m +++ b/msat/MS_phasevels.m @@ -85,51 +85,51 @@ function [pol,avs,vs1,vs2,vp,S1P,S2P,PE,S1E,S2E,XIS] = MS_phasevels(C,rh,inc,azi) - if (length(inc)~=length(azi)) - error('MS:ListsMustMatch', ... - 'AZI and INC must be scalars or vectors of the same dimension'); - end + if (length(inc)~=length(azi)) + error('MS:ListsMustMatch', ... + 'AZI and INC must be scalars or vectors of the same dimension'); + end % ** convert inc, azi to column vectors if necessary - inc = reshape(inc,length(inc),1) ; - azi = reshape(azi,length(azi),1) ; - - isotol = sqrt(eps); % Mbars - - % Check that C is valid (if check not suppressed) - MS_checkC(C); - - % ** convert GPa to MB file units (Mbars), density to g/cc - - C(:,:) = C(:,:) * 0.01 ; - rh = rh ./ 1e3 ; - - avs = zeros(size(azi)) ; - vp = zeros(size(azi)) ; - vs1 = zeros(size(azi)) ; - vs2 = zeros(size(azi)) ; - pol = zeros(size(azi)) ; - S1P = zeros(length(azi),3) ; - S2P = zeros(length(azi),3) ; - - % Eigenvectors - PE = zeros(length(azi),3) ; - S1E = zeros(length(azi),3) ; - S2E = zeros(length(azi),3) ; - % Cartesion propagation vectors - XIS = zeros(length(azi),3) ; + inc = reshape(inc,length(inc),1) ; + azi = reshape(azi,length(azi),1) ; + + isotol = sqrt(eps); % Mbars + + % Check that C is valid (if check not suppressed) + MS_checkC(C); + + % ** convert GPa to MB file units (Mbars), density to g/cc + + C(:,:) = C(:,:) * 0.01 ; + rh = rh ./ 1e3 ; + + avs = zeros(size(azi)) ; + vp = zeros(size(azi)) ; + vs1 = zeros(size(azi)) ; + vs2 = zeros(size(azi)) ; + pol = zeros(size(azi)) ; + S1P = zeros(length(azi),3) ; + S2P = zeros(length(azi),3) ; + + % Eigenvectors + PE = zeros(length(azi),3) ; + S1E = zeros(length(azi),3) ; + S2E = zeros(length(azi),3) ; + % Cartesion propagation vectors + XIS = zeros(length(azi),3) ; % ** Handle isotropic case quickly - if isIsotropic(C, isotol) - vp(:) = sqrt(( ((1.0/3.0)*(C(1,1)+2*C(1,2)))+ ... - ((4.0/3.0)*C(4,4)) )/rh)*10.0; - vs1(:) = sqrt(C(4,4)/rh)*10.0; % Factor of 10 converts from - vs2 = vs1; % Mbar to Pa. - avs(:) = 0.0; - pol(:) = NaN; % Both waves have same velocity... meaningless. - S1P(:) = NaN; - return - end + if isIsotropic(C, isotol) + vp(:) = sqrt(( ((1.0/3.0)*(C(1,1)+2*C(1,2)))+ ... + ((4.0/3.0)*C(4,4)) )/rh)*10.0; + vs1(:) = sqrt(C(4,4)/rh)*10.0; % Factor of 10 converts from + vs2 = vs1; % Mbar to Pa. + avs(:) = 0.0; + pol(:) = NaN; % Both waves have same velocity... meaningless. + S1P(:) = NaN; + return + end % ** start looping for ipair = 1:length(inc) @@ -144,8 +144,8 @@ [V,EIGVEC]=velo(XI,rh,C) ; % ** pull out the eigenvectors - P = EIGVEC(:,1) ; - S1 = EIGVEC(:,2) ; + P = EIGVEC(:,1) ; + S1 = EIGVEC(:,2) ; if ~isreal(S1) error_str = ['The S1 polarisation vector is not real!\n\n'... @@ -159,88 +159,53 @@ sprintf('S1 = %f %f %f\n',S1)]; error('MS:PHASEVELS:vectornotreal', error_str) ; end - S2 = EIGVEC(:,3) ; - - PE(ipair,:) = P ; - S1E(ipair,:) = S1 ; - S2E(ipair,:) = S2 ; - + S2 = EIGVEC(:,3) ; + + PE(ipair,:) = P ; + S1E(ipair,:) = S1 ; + S2E(ipair,:) = S2 ; + % ** calculate projection onto propagation plane - S1N = V_cross(XI,S1) ; - S1P(ipair,:) = V_cross(XI,S1N); - S2N = V_cross(XI,S2) ; - S2P(ipair,:) = V_cross(XI,S2N); + S1N = V_cross(XI,S1) ; + S1P(ipair,:) = V_cross(XI,S1N); + S2N = V_cross(XI,S2) ; + S2P(ipair,:) = V_cross(XI,S2N); % ** rotate into y-z plane to calculate angles % (use functions optimised for the two needed % rotations, see below). - [S1PR] = V_rot_gam(S1P(ipair,:),cazi) ; - [S1PRR] = V_rot_bet(S1PR,cinc) ; + [S1PR] = V_rot_gam(S1P(ipair,:),cazi) ; + [S1PRR] = V_rot_bet(S1PR,cinc) ; - ph = atan2(S1PRR(2),S1PRR(3)) .* 180/pi ; + ph = atan2(S1PRR(2),S1PRR(3)) .* 180/pi ; % ** transform angle to between -90 and 90 - if (ph < -90.), ph = ph + 180.;end - if (ph > 90.), ph = ph - 180.;end + if (ph < -90.), ph = ph + 180.;end + if (ph > 90.), ph = ph - 180.;end % ** calculate some useful values - dVS = (V(2)-V(3)) ; - VSmean = (V(2)+V(3))/2.0 ; - - avs(ipair) = 100.0*(dVS/VSmean) ; - vp(ipair) = V(1) ; - vs1(ipair) = V(2) ; - vs2(ipair) = V(3) ; - - pol(ipair) = ph ; - end % ipair = 1:length(inc_in) - - % If any directions have zero avs (within machine accuracy) - % set pol to NaN - array wise: - isiso = real(avs > sqrt(eps)) ; % list of 1.0 and 0.0. - pol = pol .* (isiso./isiso) ; % times by 1.0 or NaN. - - S1P(:,1) = S1P(:,1) .* (isiso./isiso); - S1P(:,2) = S1P(:,2) .* (isiso./isiso); - S1P(:,3) = S1P(:,3) .* (isiso./isiso); - - % switch nargout - % case 5 - % varargout{1} = pol ; - % varargout{2} = avs ; - % varargout{3} = vs1 ; - % varargout{4} = vs2 ; - % varargout{5} = vp ; - % - % case 7 - % varargout{1} = pol ; - % varargout{2} = avs ; - % varargout{3} = vs1 ; - % varargout{4} = vs2 ; - % varargout{5} = vp ; - % varargout{6} = S1P ; - % varargout{7} = S2P ; - % - % case 11 - % varargout{1} = pol ; - % varargout{2} = avs ; - % varargout{3} = vs1 ; - % varargout{4} = vs2 ; - % varargout{5} = vp ; - % varargout{6} = S1P ; - % varargout{7} = S2P ; - % varargout{8} = PE ; - % varargout{9} = S1E ; - % varargout{10} = S2E ; - % varargout{11} = XIS ; - % - % otherwise - % error('MS:PHASEVELS:BadOutputArgs','Requires 5, 7, or 11 output arguments.') - % end + dVS = (V(2)-V(3)) ; + VSmean = (V(2)+V(3))/2.0 ; + + avs(ipair) = 100.0*(dVS/VSmean) ; + vp(ipair) = V(1) ; + vs1(ipair) = V(2) ; + vs2(ipair) = V(3) ; + + pol(ipair) = ph ; + end % ipair = 1:length(inc_in) + +% If any directions have zero avs (within machine accuracy) +% set pol to NaN - array wise: + isiso = real(avs > sqrt(eps)) ; % list of 1.0 and 0.0. + pol = pol .* (isiso./isiso) ; % times by 1.0 or NaN. + + S1P(:,1) = S1P(:,1) .* (isiso./isiso); + S1P(:,2) = S1P(:,2) .* (isiso./isiso); + S1P(:,3) = S1P(:,3) .* (isiso./isiso); - return %======================================================================================= @@ -271,28 +236,28 @@ function [VR] = V_rot_gam(V,gam) % Make rotation matrix -g = gam * pi/180. ; - -RR = [ cos(g) sin(g) 0 ; -sin(g) cos(g) 0 ; 0 0 1 ] ; -VR = V * RR ; + g = gam * pi/180. ; + + RR = [ cos(g) sin(g) 0 ; -sin(g) cos(g) 0 ; 0 0 1 ] ; + VR = V * RR ; return function [VR] = V_rot_bet(V,bet) % Make rotation matrix -b = bet * pi/180. ; - -RR = [ cos(b) 0 -sin(b) ; 0 1 0 ; sin(b) 0 cos(b) ] ; - -VR = V * RR ; + b = bet * pi/180. ; + + RR = [ cos(b) 0 -sin(b) ; 0 1 0 ; sin(b) 0 cos(b) ] ; + + VR = V * RR ; return %======================================================================================= %======================================================================================= - function [X] = cart2(inc,azm) +function [X] = cart2(inc,azm) %======================================================================================= %c convert from spherical to cartesian co-ordinates %c north x=100 west y=010 up z=001 @@ -312,11 +277,11 @@ r=sqrt(X(1)*X(1)+X(2)*X(2)+X(3)*X(3)) ; X = X./r ; - return +return %======================================================================================= %======================================================================================= - function [V,EIGVEC]=velo(X,rh,C) +function [V,EIGVEC]=velo(X,rh,C) %======================================================================================= % PHASE-VELOCITY SURFACES IN AN ANISOTROPIC MEDIUM % revised April 1991 @@ -332,13 +297,13 @@ % from 6x6 to 3x3x3x3 form using the formula from page 1076 of % Winterstein 1990. This is > twice as fast as the quickest way I % have found going via the full tensor form. - gamma = [X(1) 0.0 0.0 0.0 X(3) X(2) ; ... - 0.0 X(2) 0.0 X(3) 0.0 X(1) ; ... - 0.0 0.0 X(3) X(2) X(1) 0.0 ]; - T = gamma * C * gamma'; - + gamma = [X(1) 0.0 0.0 0.0 X(3) X(2) ; ... + 0.0 X(2) 0.0 X(3) 0.0 X(1) ; ... + 0.0 0.0 X(3) X(2) X(1) 0.0 ]; + T = gamma * C * gamma'; + % determine the eigenvalues of symmetric tij - [EIVEC EIVAL] = eig(T) ; + [EIVEC EIVAL] = eig(T) ; % calculate velocities and sort % note that we could get a significant speedup if @@ -349,14 +314,14 @@ % that the eignevalues are not returned in sorted order" % http://www.mathworks.com/matlabcentral/newsreader/view_original/394371 % so we have to sort here... - V_RAW = (sqrt([EIVAL(1,1) EIVAL(2,2) EIVAL(3,3)]./rh))*10. ; - [V IND] = sort(V_RAW,2,'descend') ; - EIGVEC = EIVEC ; % for dimensioning - for i=1:3 - EIGVEC(:,i) = EIVEC(:,IND(i)) ; - end + V_RAW = (sqrt([EIVAL(1,1) EIVAL(2,2) EIVAL(3,3)]./rh))*10. ; + [V IND] = sort(V_RAW,2,'descend') ; + EIGVEC = EIVEC ; % for dimensioning + for i=1:3 + EIGVEC(:,i) = EIVEC(:,IND(i)) ; + end - return +return %======================================================================================= function [ l ] = isIsotropic( C, tol ) From a1eada44dd4064ab396354e59cff096bafa16d9a Mon Sep 17 00:00:00 2001 From: Alan Baird Date: Fri, 11 Mar 2016 13:31:09 +0000 Subject: [PATCH 7/7] added ref for groupvels and updated phasevel documentation --- msat/MS_groupvels.m | 13 ++++++++++++- msat/MS_phasevels.m | 8 ++++++-- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/msat/MS_groupvels.m b/msat/MS_groupvels.m index 4e61d5e..4c029f4 100644 --- a/msat/MS_groupvels.m +++ b/msat/MS_groupvels.m @@ -1,4 +1,4 @@ -% MS_PHASEVELS - Wave velocities in anisotropic media. +% MS_GROUPVELS - Wave velocities in anisotropic media. % % // Part of MSAT - The Matlab Seismic Anisotropy Toolkit // % @@ -20,6 +20,17 @@ % [ VGP, VGS1, VGS2, PE, S1E, S2E, SNP, SNS1, SNS2 ] = MS_groupvels( C, rh, inc, azi ) % Additionally output P, S1 and S2 Slowness vectors % +% +% Notes: +% This code is based on EMATRIX6 by D. Mainprice, recoded in +% MATLAB by Alan Baird. +% +% Reference: Mainprice D. (1990). An efficient +% FORTRAN program to calculate seismic anisotropy from +% the lattice preferred orientation of minerals. +% Computers & Gesosciences, vol16, pp385-393. +% +% See also: MS_PHASEVELS % Copyright (c) 2016, Alan Baird % All rights reserved. diff --git a/msat/MS_phasevels.m b/msat/MS_phasevels.m index 2e00b81..7c516b9 100644 --- a/msat/MS_phasevels.m +++ b/msat/MS_phasevels.m @@ -14,8 +14,12 @@ % 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. +% Additionally output fast and slow S-wave polarisation projected onto +% the plane normal to the propagation direction in vector form. +% +% [ pol, avs, vs1, vs2, vp, SF, SS, PE, S1E, S2E, XIS ] = MS_phasevels( C, rh, inc, azi ) +% Additionally output P, S1 and S2-wave polarisations in vector +% form, and the plane wave propagation vector. % % Notes: % Azi is defined as the angle in degrees from the +ve 1-axis in x1-x2