So that we can easily compare two materials.

It turns out that we don't document the plot orientation in MS_plot. We should and I recently sent this to a user:

In a sense the orientations in MSAT are up to the user. But if you generate a phase velocity
plot the directions are fixed. The positive X1 direction is up (12 o'clock), the positive X2
direction is to the left (9 o'clock) and the positive X3 direction is towards you, out of the
middle of the plot.

add this to the documentation

MS_rotEuler (and presumably the other rotation functions) allows row vectors for the three sets of angles but not column vectors. When using my scripts aimed at the VPSC code this means I need to do:

MS_rotEuler(C, eulers(1,:)', eulers(2,:)', eulers(3,:)')

Should handle this case cleanly in the code.

Can we add a small circle showing the SWS window (with an optional angular parameter) to MS_plot?

It would be nice to be able to customise the MS_plot plots to allow e.g. plots of the S1 and S2 phase velocities as well as or instead of the S wave anisotropy. We would probably want to do this alongside a better way of returning bits of the figure.

It is not clear how to input the layer thicknesses (vector) in the theory by Backus (1962) in MS_effective_medium. The error happens, for example, if I use "[1 2]" for 1m layers separated by 2m layers.

The documentation about the unit of cd (crack density) in the theory by "Hudson (1980, 1981) (hudson or crack) Cracks" in "MS_effective_medium" is not clear for me. Is it the density of cracks in the medium (the number of cracks per m^3)? I also cannot find information of the shape, size and direction of cracks and the characteristics of the medium.

We should add qz and other common crustal minerals.

I added support to MS_TI for Anderson's 1961 parameterisation in 68ec9df but we don't have test cases, or the inverse function in MS_TI_parameters. We should add these. It is probably time to rework MS_TI_parameters so that the output gets put in varargout rather than dumping a massive list of outputs (we can do this without breaking backwards compatibility).

The documentation in MS_effective_splitting_N does not specify the units of the splitting parameters, frequency, or initial polarization angle. This should be fixed and the rest of the documentation checked for similar issues.

Report from Alan: "I just tried your MS_effective_medium insaff option, and I may have found a bug. I think your znh input parameter, which you say is equal to zn/zt is actually zt/zn (the reciprocal)." - presumably this is the unmatched insaff branch.

Behaviour relating to overlaying data ('sdata' and 'pdata' options) in MS_plot is inconsistent. Current behaviour is that data is plotted in the same location regardless of the sign of the data incidence angle. However, for a default (upper hemisphere) plot, if the incidence angle is negative (i.e., it is on the lower hemisphere), this is not the correct place to compare data to the tensor (it should be placed at +180 degrees azimuth, and for splitting data the polarisation transformed).
Options to resolve:
(a) suppress plotting (with a warning) of data in the wrong hemisphere
(b) automatically transform data to be on the correct hemisphere
(c) default to (a) while providing a optional flag for (b)

Just like MS_plot, we should avoid evaling varargin in MS_sphere.

There are examples in the paper that can be made into test - and I've had a report that one of these does not work.

We don't run the C (elasticity matrix) input argument to MS_Phasevels through MS_checkC. If invalid (e.g. non-positive definite) C is provided we can then return complex numbers for the anisotropy - which breaks e.g. MS_plot. We may want to make this check optional.

One of the tests in test_MS_splitting fails on Matlab version R2013a but passes on R2009b. Specifically, the test around line 29:

 [time, T00, T90] = MS_make_trace(30.0, 0.2, 4.0);
 [~, tlag] = MS_measure_trace_splitting(time, T00, T90, 5.0, ...
      0.5, 4.0);
 assertElementsAlmostEqual([0.0 tlag], [0.0 0.0]) ;

gives tlag of 1.5 for R2013a - tlag is 0.0 on R2009b. The only obvious difference between the two version is that COVM(1,1) is -0.00 on R2013a and 0.00 on R2009b but I'm not sure why this should vary the behaviour (or, why the other three initial polarisation directions do not show the same difference.

Thank you for providing this really useful program!

I have one question about the example code you've offered (split_model.m).
According to the code of generating elasticities for lower layer (L2), it contains 2-L2_faln term in MS_VRH function, which is different with that of upper layer (L1) 1-L1_faln.

However, since L1_faln and L2_faln are both defined as volume fractions (up to my understanding) given as 0.3 in examples, I was confused why 2-L2_faln was used for lower layer instead of 1-L2_faln.

I attach below the part of the code related to my inquiry:

`%  ** generate layer elasticities:
%        This is a Voigt-Reuss-Hill average of the appropriately rotated olivine
%        tensor and its isotropic equivalent.
      [L1_C,~] = MS_VRH([L1_faln 1-L1_faln],...
         MS_rot3(Cani,0,-L1_dip,L1_aaz,'order',[3 2 1]),rh, Ciso, rh) ;

      [L2_C,~] = MS_VRH([L2_faln 2-L2_faln],...
         MS_rot3(Cani,0,-L2_dip,L2_aaz,'order',[3 2 1]),rh, Ciso, rh) ;`

It would be useful to be able to calculate group velocities (as opposed to just phase velocities) using MSAT. I've already written some python code to do this, and would be happy to rewrite it in matlab. Would it be better to write it as a new stand alone function (e.g. MS_groupvels) or perhaps as an optional output for MS_phasevels?

This was queried by a user, my answer was:

I think the vector is in the coordinate system you are after and you don't need to transform it. The
reference frame for the unit vectors describing the polarisation direction (displacement vector) of the
fast quasi-shear wave, SF (and the slow one, SS) is the same as that used to describe the elastic
constants tensor. So if SF = [ 1 0 0 ] and SS = [ 0 1 0 ] the sense of particle motion is to strain the
material in the sense resisted by C(1,1) and C(2,2), respectively. The wave propagation direction will
be along (or, if the symmetry is low enough, close to) the X3 direction (inc = 90.0 and the result should > be independent of azimuth).

but some clarification in the documentation would be useful.

Also, why don't we output the qP polarisation direction?

I am looking at the two example matrices in Browaeys and Chevrot (2004) [BC2004], for olivine (Section 4.1.1) and for enstatite (Section 4.1.2). I'll call these 6x6 Voigt matrices Colivine and Censtatite. The suggested approach in MSAT to the BC2004 decomposition is (see MSAT Walker paper, Listing 2): MS_axes, MS_decomp, MS_norms. I was expecting both of the BC2004 matrices to be in the SCCS (symmetry cartesian coordinate system) form (Section 3.2), which is supposedly achieved by using MS_axes.m. But when I run MS_axes(Colivine) and MS_axes(Censtatite), I get different matrices from the input ones. I'll refer to these as Colivine_MSaxes and Censtatite_MSaxes.

Running MS_decomp/MS_norms on Colivine and Colivine_MSaxes gives the same percentages as those listed in BC2004. Running MS_decomp/MS_norms on Censtatite gives the same percentages as those listed in the paper. Running MS_decomp/MS_norms on on Censtatite_MSaxes gives different percentages.

This may be more of a question for Browaeys and Chevrot, but I wanted to post the issue here, given that the BC decomposition is implemented in MSAT.

With Colivine and Censtatite devined, the Matlab code is

Colivine
Colivine_MSaxes = MS_axes(Colivine)
[Ciso,Chex,Ctet,Cort,Cmon,Ctri] = MS_decomp(Colivine);
Polivine = MS_norms(Colivine,Ciso,Chex,Ctet,Cort,Cmon,Ctri)
[Ciso,Chex,Ctet,Cort,Cmon,Ctri] = MS_decomp(Colivine_MSaxes);
Polivine_MSaxes = MS_norms(Colivine_MSaxes,Ciso,Chex,Ctet,Cort,Cmon,Ctri)

Censtatite
Censtatite_MSaxes = MS_axes(Censtatite)
[Ciso,Chex,Ctet,Cort,Cmon,Ctri] = MS_decomp(Censtatite);
Penstatite = MS_norms(Censtatite,Ciso,Chex,Ctet,Cort,Cmon,Ctri)
[Ciso,Chex,Ctet,Cort,Cmon,Ctri] = MS_decomp(Censtatite_MSaxes);
Penstatite_MSaxes = MS_norms(Censtatite_MSaxes,Ciso,Chex,Ctet,Cort,Cmon,Ctri)

The output is

Colivine =
   192    66    60     0     0     0
    66   160    56     0     0     0
    60    56   272     0     0     0
     0     0     0    60     0     0
     0     0     0     0    62     0
     0     0     0     0     0    49
Colivine_MSaxes =
   160    66    56     0     0     0
    66   192    60     0     0     0
    56    60   272     0     0     0
     0     0     0    62     0     0
     0     0     0     0    60     0
     0     0     0     0     0    49
Polivine =
    0.7930    0.1516    0.0034    0.0520         0         0
Polivine_MSaxes =
    0.7930    0.1516    0.0034    0.0520         0         0
Censtatite =
   225    54    72     0     0     0
    54   214    53     0     0     0
    72    53   178     0     0     0
     0     0     0    78     0     0
     0     0     0     0    82     0
     0     0     0     0     0    76
Censtatite_MSaxes =
   178    53    72     0     0     0
    53   214    54     0     0     0
    72    54   225     0     0     0
     0     0     0    76     0     0
     0     0     0     0    82     0
     0     0     0     0     0    78
Penstatite =
    0.9081    0.0426    0.0043    0.0450         0         0
Penstatite_MSaxes =
    0.9081    0.0210    0.0027    0.0682         0         0
>> 

The following has been reported:

% Make VTI tensor
A=MS_VTI(2000,1000,1000,0.4,0.3,0.3)
% Do MS_axes
[B,r1] = MS_axes(A)
% compare A and B 
MS_plot(A,1000)
MS_plot(B,1000) % notice the girdle in the plane of view, symmetry axis along x1.
% now do Browaeys decomposition.
[Aiso,Ahex] = MS_decomp(A) % do it on the starting tensor which should be in the right orientation for comparison.
[Biso,Bhex] = MS_decomp(B)
% look at results 
MS_plot(Aiso+Ahex,1000) % looks as expected
MS_plot(Biso+Bhex,1000) % looks odd
% Notice the Biso+Bhex has a completely different anisotropy to what we started with (A).

After a little debugging it turns out that numd from:

[nud,id] = ndistinct(valv,norm(valv)./1000) ;

returns '3' not '2'. Reducing the threshold (to norm(valv), which seems huge) fixes this. We either have a bug in ndistinct, or need to estimate the accuracy of the eigenvalues.

This seems to be working on my interp branch. However we need an example. And I need to think about tests for ortho.

Example

[Cjd, den] = MS_elasticDB('jd')
[Cdi, ~] = MS_elasticDB('di')

CjdR = MS_rot3(Cjd, 30, 30, 0);

% Element wise interpolation
[~, ~, CV, ~] = MS_VRH([0.5 0.5], CjdR, 1, Cdi, 1)

% Proper interp
Cint = MS_interpolate(CjdR, Cdi, 0.5)

MS_plot(Cdi, den)
MS_plot(CjdR, den)
MS_plot(Cint, den)
MS_plot(CV, den)

Tests

The first (but not the next two) of these work,
and I don't know if any should... Look OK using the
same as above.

function test_MS_interpolate_values_ortho
% Test the interpolator works for different
% matrices without rotation
[Col, ] = MS_elasticDB('olivine');
[Cens, ] = MS_elasticDB('enstatite');
[,,Cmix25,] = MS_VRH([0.25, 0.75], Col, 1.0, Cens, 1.0);
[,,Cmix50,] = MS_VRH([0.5, 0.5], Col, 1.0, Cens, 1.0);
[,,Cmix67,~] = MS_VRH([0.67, 0.33], Col, 1.0, Cens, 1.0);
assertElementsAlmostEqual(MS_interpolate(Col, Cens, 0.5), Cmix50);
assertElementsAlmostEqual(MS_interpolate(Cens, Col, 0.5), Cmix50);
assertElementsAlmostEqual(MS_interpolate(Col, Cens, 0.25), Cmix25);
assertElementsAlmostEqual(MS_interpolate(Col, Cens, 0.67), Cmix67);
end

function test_MS_interpolate_values_ang_ortho_x
% Test the interpolator works for different
% matrices with rotation
[Col, ] = MS_elasticDB('olivine');
[Cens, ] = MS_elasticDB('enstatite');
CensR = MS_rot3(Cens, 45, 0, 0);
[,,Cmix50,~] = MS_VRH([0.5, 0.5], Col, 1.0, Cens, 1.0);
Cmix50R = MS_rot3(Cmix50, 22.5, 0, 0);
assertElementsAlmostEqual(MS_interpolate(Col, CensR, 0.5), Cmix50R);
assertElementsAlmostEqual(MS_interpolate(CensR, Col, 0.5), Cmix50R);
end

function test_MS_interpolate_values_ang_ortho_y
% Test the interpolator works for different
% matrices with rotation
[Col, ] = MS_elasticDB('olivine');
[Cens, ] = MS_elasticDB('enstatite');
CensR = MS_rot3(Cens, 0, 45, 0);
[,,Cmix50,~] = MS_VRH([0.5, 0.5], Col, 1.0, Cens, 1.0);
Cmix50R = MS_rot3(Cmix50, 0, 22.5, 0);
assertElementsAlmostEqual(MS_interpolate(Col, CensR, 0.5), Cmix50R);
%assertElementsAlmostEqual(MS_interpolate(CensR, Col, 0.5), Cmix50R);
end

function test_MS_interpolate_values_ang_ortho_z
% Test the interpolator works for different
% matrices with rotation
[Col, ] = MS_elasticDB('olivine');
[Cens, ] = MS_elasticDB('enstatite');
CensR = MS_rot3(Cens, 0, 0, 45);
[,,Cmix50,~] = MS_VRH([0.5, 0.5], Col, 1.0, Cens, 1.0);
Cmix50R = MS_rot3(Cmix50, 0, 0, -22.5);
assertElementsAlmostEqual(MS_interpolate(Col, CensR, 0.5), Cmix50R);
%assertElementsAlmostEqual(MS_interpolate(CensR, Col, 0.5), Cmix50R);
end

Thanks for a great program! However, one always wants more...

I think there is a problem with the documentation for MS_plot when plotting one's own data, in terms of the forms of the matrices/vectors called. The document for how to plot your own observations on the stereonets for MS_plot says this:

MS_plot(..., 'sdata', azimuth, inclination, polarisation, avs)
Add data points to the P-wave velocity ('VP') or S-wave
polarisation plot ('AVSPOL'), respectively. This can be used,
for example, to compare an elastic model with shear-wave
splitting measurements. In each case data points relate
to a ray propagating in a direction described by the azimuth
and inclination. The point is coloured according to either
the P-wave velocity, vp, or S-wave anisotropy, avs, in km/s
or %. In the case of S-wave anisotropy, the fast polarisation
direction is also shown. Each argument (azimuth, inclination,
vp, polarisation and avs) is an array and each must be the
same length.

However, the program bombs if they are all the same 1xn or nx1 dimensions. I found I needed to make the variables azimuth and inclination have the same dimensions, but polarisation and avs need to be the transpose, i.e. if azimuth and inclination are each 1xn, then polarisation and avs need to be nx1.

I have moved to Leeds and have a new email address and web page. This needs updating on the website and possibly here too.