andreww
commented
3 years ago Yes - this needs to be fixed in the documentation. But, more importantly, I think we need some error checking in the code too. You are supposed to be able to use this in two ways. In the first way you should be able to pass in three 1D arrays with each element corresponding to the thickness, P-wave velocity and S-wave velocity of a layer. For example for repeated layers of three materials where the first material is 1.1 m thick has Vp of 6.1 km/s and Vs of 3.1 km/s, and density 2.1 kg/m^3, a second layer of thickness 10.2 m with Vp = 6.2 km/s, Vs = 3.2 km/s, density = 2.2 kg/m^3 and a third layer with thickness 2.3 m with Vp of 6.3 km/s, Vs of 3.3 km/s and density = 2.3 kg/m^3 the effective elastic constants would be generated by doing:
[Ceff, rh] = MS_effective_medium('backus', [1.1 10.2 2.3], [6.1 6.2 6.3], [3.1 3.2 3.3], [2.1 2.2 2.3])The symmetry of Ceff will be hexagonal with the with the unique axis pointing in the z (3) direction. This means that the layers are stacked in this direction (so if you think of z as depth the layers are horizontal). The units of the thickness do not actually matter (as far as I can see) but the other units do if you want Ceff in GPa and rh in kg/m^3 (you do want this).
Now, the problem is that the code does not check that the lengths of those four arrays are the same. It ought to, probably in this if clause. At present passing in a scalar to one of those arguments will work (I think) and (if I'm following the code right) will treat that value as the same for each layer. However, I don't think anybody will really want to do that and different length arrays should lead to a runtime crash...
We probably need some more tests two. I'll put this on my list of things that need fixing.
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.