search

mtex-toolbox / mtex

MTEX is a free Matlab toolbox for quantitative texture analysis. Homepage:
http://mtex-toolbox.github.io/
GNU General Public License v2.0
286 stars 184 forks source link

Problem with velocity using C.velocity('harmonic') #417

Closed DavidMainprice closed 4 years ago

DavidMainprice commented 5 years ago

MAC os Matlab 2015b MTEX 5.2.0.beta3

When using with this elastic constants the value for max(Vp) is bigger than MATLAB can handle, results MATLAB crashes. By running the script using XY_grid = equispacedS2Grid('upper','resolution',1*degree); [vp,vs1,vs2,~,~,~] = velocity(C,XY_grid); it does not crash and all results are good.

%** % Define elastic stiffness tensor (GPa) %** % % Lasted Checked : 21/05/2019 % % Reference Elastic constants % A.B. Belonoshko, N.V. Skorodumova, S.Davis, A.N.Osiptsov, A.Rosengren, B.Johansson,(2007) % Origin of the Low Rigidity of the Earths Inner Core % SCIENCE VOL 316 15 JUNE 2007 pp.1603-1605. % % Cylindrical average about axis [111] parallel to Z-axis or Earth rotation axis % Specimen symmetry reduced to triclinic (-1) % % Define density (g/cm3) rho=13.559; % % Define Cartesian Tensor crystal symmetry & frame cs_Tensor = crystalSymmetry('-1',[2.866 2.866 2.866],... [90.00 90.00 90.00]*degree,'x||a','z||c',... 'mineral','bcc-iron P=356.7 GPa T=6000K Cylindrical averaged [111] 2007'); % % elastic Cij stiffness tensor (GPa) as matrix M M =... [[ 1870.3 1345.2 1242.3 0.0 0.0 0.0];... [ 1345.2 1870.3 1345.2 0.0 0.0 0.0];... [ 1242.3 1345.2 1870.3 0.0 0.0 0.0];... [ 0.0 0.0 0.0 1973.2 0.0 0.0];... [ 0.0 0.0 0.0 0.0 365.1 0.0];... [ 0.0 0.0 0.0 0.0 0.0 365.1]]; % % M as stiffness tensor C with MTEX tensor command C = stiffnessTensor(M,cs_Tensor,'density',rho) % % Calculate seismic anisotropic 3 parameters: AVp,maxVs,minAVs % Tested with MTEX 5.2.beta3 and MATLAB 2016a % David Mainprice 23/05/2019 % INPUT % C = stiffnessTensor(M,cs_Tensor,'density',rho); %** % Equations of 3 anisotropic parameters based on phase velocities % AVp(%) = 200(maxVp-minVp)./(maxVp+minVp) % maxAVs(%) = 200(maxS1-minS1)./(maxS1+minS1) % minAVs(%) = min(AVs) %** % Calculate phase velocities using harmonic method with crystal symmetry %[vp,vs1,vs2,~,~,~] = C.velocity('harmonic') % ERROR using C.velocity('harmonic') max(Vp) is higher number than MATLAB % can handle! and crashes !!! XY_grid = equispacedS2Grid('upper','resolution',1degree); [vp,vs1,vs2,~,~,~] = velocity(C,XY_grid); % AVp - percentage anisotropy (%) maxVp = max(vp) minVp = min(vp) AVp = 200(maxVp-minVp)./(maxVp+minVp) % S-wave anisotropy (%) AVs = 200*(vs1-vs2)./(vs1+vs2); maxAVs = max(AVs) minAVs = min(AVs) % If minAVs is negative set to zero if(minAVs<0) minAVs=0; end

all the best David

kilir commented 5 years ago

Hi David, on macOS 10.13 Matlab 2018b, latest mtex development version, I don't have any issues with the following syntax. 

[vp,vs1,vs2] = C.velocity('harmonic')

vp = S2FunHarmonic (show methods, plot)
 mineral: bcc-iron P=356.7 GPa T=6000K Cylindrical averaged [111] 2007 (-1, X||a*, Y||b, Z||c*)
 bandwidth: 48
 antipodal: true

vs1 = S2FunHarmonic (show methods, plot)
 mineral: bcc-iron P=356.7 GPa T=6000K Cylindrical averaged [111] 2007 (-1, X||a*, Y||b, Z||c*)
 bandwidth: 48
 antipodal: true

vs2 = S2FunHarmonic (show methods, plot)
 mineral: bcc-iron P=356.7 GPa T=6000K Cylindrical averaged [111] 2007 (-1, X||a*, Y||b, Z||c*)
 bandwidth: 48
 antipodal: true

Anything more of information you could add? Any chance to try it with a more recent Matlab version?

Cheers, Rüdiger