Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Groupvel #30

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Groupvel #30

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
803946b
added .DS_Stote to .gitignore
alanfbaird Feb 5, 2016
d57a98d
added MS_groupvels (copied from MS_Phasevels with added function rayvel
alanfbaird Feb 5, 2016
cae5a84
MS_groupvels should now output group velocities
alanfbaird Feb 5, 2016
349301c
updated phasevels to output eigenvectors and prop vector, used it in …
alanfbaird Mar 4, 2016
6270c98
fixed some bugs in MS_groupvels
alanfbaird Mar 10, 2016
e289626
tidied up indentation
alanfbaird Mar 11, 2016
a1eada4
added ref for groupvels and updated phasevel documentation
alanfbaird Mar 11, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
*.m~
.DS_Store
204 changes: 204 additions & 0 deletions msat/MS_groupvels.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,204 @@
% 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
%
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

any reference to add?

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also add
% see also: MS_phasevels


% 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) ;
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you tidy up the indentation?

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) ;




% ** start looping
for ipair = 1:length(inc)

% ** slowness vectors
SNP(ipair,:) = XIS(ipair,:)./vp(ipair);
SNS1(ipair,:) = XIS(ipair,:)./vs1(ipair);
SNS2(ipair,:) = XIS(ipair,:)./vs2(ipair);

% ** Group velocity vectors

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

switch nargout
case 3
varargout{1} = VGP ;
varargout{2} = VGS1 ;
varargout{3} = VGS2 ;


case 6
varargout{1} = VGP ;
varargout{2} = VGS1 ;
varargout{3} = VGS2 ;
varargout{4} = PE ;
varargout{5} = S1E ;
varargout{6} = S2E ;


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 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 = [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);


% 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,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

% 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

return % function




62 changes: 54 additions & 8 deletions msat/MS_phasevels.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you add documentation for PE, S1E, S2E and XIS?


if (length(inc)~=length(azi))
error('MS:ListsMustMatch', ...
Expand All @@ -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)
Expand All @@ -132,6 +138,7 @@

% ** create the cartesian vector
XI = cart2(cinc,cazi) ;
XIS(ipair,:)=XI ;

% ** compute phase velocities
[V,EIGVEC]=velo(XI,rh,C) ;
Expand All @@ -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) ;
Expand Down Expand Up @@ -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
%=======================================================================================
Expand Down