ElasticMatrix: A MATLAB Toolbox for Anisotropic Elastic Wave Propagation in Layered Media
Danny R. Ramasawmy Ben T. Cox Bradley E. Treeby

Radio Resource Dimensioning with Cox Process Based User Location Distribution
    Ridha Nasri, 
	Jalal Rachad and
	Laurent Decreusefond
	
	Ridha Nasri and Jalal Rachad are with Orange Labs, 40-48 Avenue de la République, 92320 Châtillon, France. Emails: {ridha.nasri,jalal.rachad}@orange.com
	Jalal Rachad and Laurent Decreusefond are with LTCI Telecom Paris-Tech, université Paris Saclay, 23 Avenue d'Italie,75013, Paris, France. Emails: jalal.rachad@telecom-paristech.fr, laurent.decreusefond@mines-telecom.fr
	

    December 30, 2023
====================================================================================================================================================================================================================================================================================================================================================================================================================================================







Simulating the propagation of elastic waves in multi-layered media has many applications. A common approach is to use matrix methods where the elastic wave-field within each material layer is represented by a sum of partial-waves along with boundary conditions imposed at each interface. While these methods are well-known, coding the required matrix formation, inversion, and analysis for general multi-layered systems is non-trivial and time-consuming. Here, a new open-source toolbox called  is described which solves the problem of acoustic and elastic wave propagation in multi-layered media for isotropic and transverse-isotropic materials where the wave propagation occurs in a material plane of symmetry. The toolbox is implemented in MATLAB using an object oriented programming framework and is designed to be easy to use and extend. Methods are provided for calculating and plotting dispersion curves, displacement and stress fields, reflection and transmission coefficients, and slowness profiles.


  










§ MOTIVATION AND SIGNIFICANCE



    Matrix models of wave propagation in multi-layered elastic solids have had a significant contribution to research areas such as acoustics, geophysics and electromagnetics.  
    A few examples include: structural health monitoring <cit.>, characterisation of interface bonding <cit.>, detection of debonding in joints <cit.>, measuring material properties <cit.>, designing composite layered structures <cit.>,  mode sorting of guided waves <cit.>, the physical interpretation of guided wave structures <cit.>, modelling the directional response of Fabry-Pérot ultrasound sensors <cit.>, reflection and transmission of plane waves <cit.>, elastography of layered soft tissues <cit.>, and ice detection on wind turbines <cit.>. 
    
    Matrix methods, in particular the partial-wave and global matrix method, represent the stress and displacement fields as a sum of partial-waves for each material of the layered-structure. Each partial-wave represents an upward or downward travelling (quasi-)compressional or (quasi-)shear wave. By invoking boundary conditions at the interfaces of adjacent layers, the partial-wave amplitudes and field properties of the first layer can be related to the last in the form of a `global' matrix. The resulting matrix equation can be used in two different ways. Firstly, the roots of the equation can be found which give the modal solutions or dispersion curves. Secondly, a subset of partial-wave amplitudes can be defined and the remaining amplitudes solved for. This can be used to calculate the displacement and stress fields within the multi-layered structure when a plane wave is incident. This method will be discussed further in Section <ref>. 
    
    Despite its usefulness, there are few available implementations of the partial-wave method. The current state-of-the-art implementation is Disperse <cit.>. This software has been in development since 1990 and is primarily focused on calculating the dispersion solutions for multi-layered structures. 
    The Disperse software was originally based the partial-wave method, however, is currently being updated to use the spectral collocation method <cit.>. The main limitation with Disperse is that it is closed-source. For this reason it is not easily adaptable for applications that are not dispersion analysis, for example, extracting reflection coefficients or slowness profiles. There currently is only a single open-source code modelling the partial-wave method, (LAMB <cit.>), however, this is limited to modelling only an isotropic plate.
    
    In this paper, a new open-source toolbox called  is introduced which uses the partial-wave method for multi-layered structures with an arbitrary number of isotropic and transverse-isotropic materials. Where possible, it is validated against existing literature and has been implemented so that it is both easy to use and extend. Some potential uses of this software are: 1) plotting the slowness profiles of materials, 2) determining the reflection and transmission coefficients of multi-layered structures, 3) finding the dispersion curves of multi-layered structures, 4) plotting the displacement and stress fields, 5) extending the toolbox for other applications, for example modelling the directional response of Fabry-Perot ultrasound sensors <cit.>. 
	
	A brief overview of the underlying mathematical model is described in Section <ref>. A selection of code snippets and examples are shown in Section <ref>. (More extensive examples are available with the toolbox documentation.) The impact and conclusions are described in Section <ref>.




  	




§ MODEL DESCRIPTION


    
    

 §.§ Overview

    
    
     uses the partial-wave method to model wave propagation in multi-layered elastic solids. The method describes elastic plane-wave propagation along a plane of symmetry for n-layers of rigidly bonded transverse-isotropic materials. 
    
    An example of an isotropic material is glass, where the material properties are the same when measured from every direction. An example of a transverse-isotropic material is a bundle of fibres, where the properties have translational symmetry axially along the fibre, and are isotropic in the plane perpendicular to this. 
    
     can model layered transverse-isotropic materials if they are aligned such that they have rotational symmetry about the axis perpendicular to the plane of each layer and the wave-vector of the propagating wave lies in the plane of symmetry. In this case, the multi-layered structure can be modelled in two-dimensions. This is illustrated in Fig. <ref>. 
    
    
    Each partial-wave represents the superposition of waves that have been multiply reflected or transmitted at the interfaces between each layer in a steady-state. The polarisation vector and wave-vector of each of these partial-waves can be found from the Christoffel Equation which is described in Section <ref>. The degree of reflection and transmission depends on the boundary conditions at the interfaces and material properties of each layer. The coupled equations that arise from the boundary conditions can be combined into a `global-matrix' which allows them to be solved simultaneously, which is discussed in Section <ref>. This global matrix approach can be used to tackle various problems in elastic wave propagation. For example, the singularities of the global matrix give the dispersion curves, and by specifying an incident wave, the resulting wave-field throughout the structure can be calculated. More detailed descriptions of the partial-wave and global-matrix method can be found in <cit.>.
    
    

 §.§ Wave-vectors and Polarisation

    
    
    Firstly, the solution for a plane wave propagating in an unbounded medium is derived. This is needed to calculate the polarisation and wave-vectors for each partial-wave component and the process is repeated independently for every layer.
	The wave-equation for an anisotropic unbounded medium is
  	
  	
    ρ∂^2 u_i/∂ t^2 = C_ijkl∂^2 u_l/∂ x_j ∂ x_k,

	
    where the indices i,j,k,l ∈{1,2,3}, x and t are the spatial and temporal variables and Einstein summation notation is used. The variable u_i is the displacement in direction i. The elastic properties of each material are described by the density ρ and the stiffness-tensor C_ijkl. The stiffness tensor has 81 components which can be reduced to 21 independent coefficients to describe a fully-anisotropic medium <cit.>. Here, the analysis is restricted to materials that are either isotropic or transverse-isotropic, which reduces the number of independent coefficients further. As described previously, the wave-vectors of the partial-waves lie in a plane of material symmetry, (x_1, x_3). Here ∂/∂ x_2 = 0 and the expanded form of the wave-equation Eq. (<ref>) is
	
    
    ρ∂^2 u_1/∂ t^2   =  
        	C_11∂^2 u_1/∂ x_1^2 + C_55∂^2 u_1/∂ x_3^2 + (C_13 + C_55)( ∂^2 u_3/∂ x_1 ∂ x_3) 
    ρ∂^2 u_2/∂ t^2   = 
        	C_66∂^2 u_2/∂ x_1^2 + C_44∂^2 u_2/∂ x_3^2
    ρ∂^2 u_3/∂ t^2   = 
        	C_55∂^2 u_3/∂ x_1^2 + C_33∂^2 u_3/∂ x_3^2 + (C_13 + C_55)( ∂^2 u_1/∂ x_1 ∂ x_2),

    
    where Voigt notation has been used to contract the indices of the stiffness-matrix (where 11 → 1, 22 → 2, 33 → 3, 23 → 4, 13 → 5, 12 → 6). A single-frequency plane wave can be written in the form
    
    
    u_i = A_i exp(𝗂(ζ x_1 + ζα x_3 - ω t)),

    
    where i ∈{ 1,2,3 }, ω is the circular frequency, α is the ratio of the vertical and horizontal (ζ) wavenumbers, and A_i is the polarisation unit vector which describes the direction of displacement relative to the direction of wave propagation. Substituting Eq. (<ref>) into Eq. (<ref>) gives the Christoffel equation
    
    
    Γ_ij(α) A_j = 0,

    where the components of the Christoffel matrix (Γ) are
    
    Γ_11   = (C_11  - ρν^2 + C_55α^2)     Γ_22   = (C_66 - ρν^2 + C_44α^2)  
    Γ_33   = (C_55 - ρν^2 + C_33α^2)    Γ_13   = Γ_31 =  (C_13 + C_55) α
    Γ_12   = Γ_21 = Γ_32 =Γ_23 = 0.

   	
    The phase velocity ν along the x_1 axis is calculated from the relation ν = ω/ζ. Solving Eq. (<ref>) admits three solutions for α^2 and therefore 6 solutions for α.  From here, the notation α_q, where q ∈{1,2,...,6 }, will be used to indicate each solution.
   
   It can be seen from Eq. (<ref>) that the plane wave component A_2 is only dependent on Γ_22, hence displacement occurring in the (x_1, x_3) plane is independent of displacement in x_2. 
   Four solutions (q = 1,2,3,4) of α_q Eq. (<ref>) can be found when
    
    
    |
            [ Γ_11(α) Γ_13(α); Γ_13(α) Γ_33(α) ]| = 0.

    
    These describe upwards and downward travelling quasi-shear-vertical (qSV) and quasi-longitudinal (qL) waves with the displacement restricted to the (x_1,x_3) plane.
    
    The remaining two solutions (q = 5,6) are found from Γ_22=0, and correspond to upward- and downward- travelling quasi-shear-horizontal (qSV) waves. 
    
    The notation A_iq will be used to indicate polarisation vector for each solution q. 
    The displacement field can now be written as
    
    u_i = ∑_q A_iq B_qexp(𝗂(ζ x_1 + ζα_q x_3 - ω t)),

    where B_q is the amplitude of each partial-wave. Additionally, the stress field within the unbounded medium can be found using Hooke's law
    
    σ_ij = C_ijkl( ∂ u_k/∂ x_l + ∂ u_l/∂ x_k),

   	where i,j,k,l ∈{1,2,3}. For a multi-layered medium, the Christoffel equation Eq. (<ref>) is solved independently for every layer to calculate the polarisation vector and wave-vector of each partial-wave. However, the amplitude B_q of each partial-wave is solved by invoking the boundary conditions at the interfaces of adjacent layers. This is discussed in the following section.
   	
   

 §.§ Boundary Conditions and Partial-wave Amplitudes

   
    As mentioned previously, the wave-vector of the plane waves are constrained to a plane of symmetry of the transverse-isotropic material reducing the analysis to two dimensions, (x_1, x_3). The normal and transverse displacement and stress describing (quasi-)longitudinal and (quasi-)shear-vertical waves for a single layer is written in the form 
   
    [ [  u_1;  u_3; σ_33; σ_13 ]] =
    	    [ [ A_11 A_12 A_13 A_14; A_31 A_32 A_33 A_34; D_11 D_12 D_13 D_14; D_21 D_22 D_23 D_24 ]]
    	     [ [ e_1            ;     e_2        ;         e_3    ;             e_4 ]]
    	    [  [ B_1; B_2; B_3; B_4 ]],

     where 
    
    D_1q   = (C_13A_1q + C_33α_q A_3q )(𝗂ζ)     D_2q   = C_55(α_qA_1
        q+  A_3q)(𝗂ζ), 
    
        e_q    = exp(𝗂 (ζ x_1 + ζα_q x_3 - ω t)).

    Only the first four solutions of α_q are needed as the motion is restricted to two dimensions (x_1, x_3).
    
    The left hand vector of Eqs. (<ref>) contains the components of the displacement and stress, and the right hand vector contains the amplitude of the partial-wave components. The product of the matrices in Eq. (<ref>) will be written as a field matrix 𝐅. At an interface x_3 = d between material layers in welded contact, the normal transverse stress and displacement must be continuous across the interface. 
    
    Therefore, the product of the field matrix and wave amplitudes at the interface of one layer is set equal to the field matrix and wave amplitudes of the adjacent layer. This process is repeated for every interface of the layered medium. For n-layers, there are 4(n-1) boundary conditions and 4n wave amplitudes which can be arranged into a global matrix. For example, for a medium consisting of 4 layers, the global matrix equation is written
    
    [ [   𝐅_1^1  -𝐅_2^1                ;           𝐅_2^2 - 𝐅_3^2        ;                   𝐅_3^3  -𝐅_4^3;         ]]
    	    [  [ 𝐁_1; 𝐁_2; 𝐁_3; 𝐁_4 ]] = 0.

  	Here, 𝐅_n^N is the 4 × 4 field matrix and 𝐁_n a 4×1 vector of partial-wave amplitudes of layer n at interface N. By assigning values to four of the partial-wave amplitudes, Eq. (<ref>), can be rearranged and solved for the remaining partial-wave amplitudes.
  	
  	For example, if a compressional wave in the first medium is incident on the layered-structure, the downward-travelling partial-wave amplitude relating to shear (B_2^1, Figure <ref>) in the first layer and upward-travelling partial-wave amplitudes relating to compressional (B_3^n) and shear (B_4^n) waves in the last (nth) layer are set to zero. Finally, the downward-travelling partial-wave amplitude relating to a compressional wave in the first layer is set to an arbitrary value (B^1_1 = 1). In this case, the solved amplitudes describe the solution for an incident single-frequency plane-wave at an angle θ or wavenumber ζ and frequency f. For a four-layered medium this would be 
  	    
    [ [  𝐅_1^1+  -𝐅_2^1                ;           𝐅_2^2 - 𝐅_3^2        ;                   𝐅_3^3 -𝐅_4^3-;         ]]
    	    [  [ 𝐁_1^+;   𝐁_2;   𝐁_3; 𝐁_4^- ]] =      [ [ -𝐅_1^1-                        ;                                ;                          𝐅_4^3+;         ]]
    	    [  [ 𝐁_1^-;      ;      ; 𝐁_4^+ ]]  .

    
    Here, + and - superscripts indicate the upwards and downwards travelling partial-wave amplitudes and their respective columns in the field matrices. For example, 𝐅_1^1+ would be the third and fourth columns of 𝐅_1^1, and 𝐁^+_1 would be the third and fourth elements of 𝐁_1. In the example described above, the first element of 𝐁_1^- is 1  and all the elements of 𝐁_4^+ and the second element of 𝐁_1^- are 0.
  	
  	Additionally, once the remaining wave-amplitudes for each layer are found, Eq. (<ref>) can be used to find the displacement and stress anywhere in the layered structure. Alternatively, the dispersion curves can be extracted from the model by setting the incident wave-amplitudes of the layered structure to zero and finding the frequency-wavenumber pairs in which the resulting left-hand-side matrix of Eq. (<ref>) becomes singular. The algorithm used is described further in Section <ref> 
  	
  	

 §.§ Shear-Horizontal Waves

    Shear-horizontal waves propagate independently of (quasi-)shear-vertical and (quasi-)compressional waves and their propagation is analogous to compressional waves in a liquid <cit.>. Solutions 5 and 6 of α_q correspond to shear-horizontal waves and the displacement and shear stress can be written in the form
    
    [ [  u_2; σ_23 ]] = 
            [ [ A_25 A_26; D_35 D_36 ]]
            [ [ e_5    ;     e_6 ]]
            [ [ B_5; B_6 ]],

    where
    
    D_3q = C_44α_q A_2q (iζ).

    By setting the properties of a new medium C' so that C'_11, C'_22, C'_33 = C_44 and setting the remaining coefficients to zero, Eq. (<ref>) reduces to Eq. (<ref>). Hence, the propogation of shear-horizontal waves can be modelled with Eq. (<ref>) by reassigning the relevant components of the stiffness matrix.
  	
  
 	

 §.§ Implementation Details

 	
 	To construct the global matrix, Eq. (<ref>), a field matrix 𝐅_𝐢 must be calculated for each layer. To improve the conditioning of the matrix, rows relating to displacement equations are scaled by ζ the horizontal wavenumber, and rows relating to stress are scaled by ρω^2. The system matrix is constructed by looping over each interface, calculating the 4×4 field matrices above and below each interface and arranging them into a single matrix. This leads to a rectangular matrix which has 4n columns and 4n-1 rows. Additionally there are 4n partial-wave amplitudes. Four partial-wave amplitudes are defined and the global-matrix in Eq. (<ref>) is rearranged to be square Eq. (<ref>). The resulting equation is solved using the  function in MATLAB. This function solves a system of linear equations using the fastest algorithm based on the matrix structure. However, the global-matrix becomes singular at values of ζ and ω on or close-to dispersion curve solutions.
  	
  	The computation of dispersion curves follows the algorithm described in <cit.>. Firstly, the wavenumber parameter is fixed and the determinant of the global-matrix is found over a range of frequencies. Close to dispersion solutions the determinant of the global-matrix tends to zero. Using these as starting points and taking a limit either side, the exact frequency and wavenumber of the dispersive solution is found using a bisection algorithm.  makes use of MATLAB's  function for this. These solutions are the starting points for each dispersion curve. To find the second point on each dispersion curve, the fixed value of wavenumber is increased and the search is performed again. The algorithm then uses linear interpolation to estimate the location of the third, fourth and fifth points on the dispersion curve, similarly using a bisection algorithm to find the exact frequency-wavenumber pairs. After five points have been found, a higher-order polynomial interpolation scheme is used to more accurately predict points on the dispersion curve. The algorithm implemented in  only searches in the real domain of ζ which is a good-estimate for simple plate structures in a vacuum, however, it may be inaccurate for leaky solutions, for example an plate embedded in soil.

    Slowness profiles are calculated by defining a range of phase-speeds and calculating the horizontal and vertical component of wavenumber by solving the Christoffel equation Eq. (<ref>). The values from the calculation may be complex, however, only the real values are plotted.








§ SOFTWARE DESCRIPTION AND EXAMPLES



    

 §.§ Overview

    The  toolbox implements the partial-wave method using an object-oriented framework in MATLAB. This allows the toolbox to be used with either a simple scripting or command line interface, and makes it easy to use and expand. The software is divided into three classes. The first class, , defines the multi-layered geometry and material properties of each layer. The second class  is initialised by a  object. This class contains the partial-wave method implementation and methods for extracting additional details such as dispersion curves and reflection coefficients. By default, all the calculations use 64 bit precision. The final class, , is an example of how  numerical models can be built from the  and  objects. This class inherits  and can be used to model the directional response of a Fabry-Pérot ultrasound sensor. More details can be found in <cit.>. Each class in the toolbox inherits the MATLAB handle class. Consequently, the object does not need to be reassigned when a method is called. The classes and their respective attributes and methods can be seen in Fig. <ref>. The toolbox is self contained and has been tested with MATLAB 2016a and above.
    
    This section presents a small selection of code snippets and examples. More detailed examples can be found in the   folder and  documentation can be accessed through the MATLAB help and clicking . There are three steps to using the toolbox. Firstly, the geometry of the layered medium must be defined. Secondly, the input parameters to the model should be defined, which are generally a range of angles, frequencies or wavenumbers. Finally, the model can be solved and details such as the reflection coefficients and dispersion curves can be extracted. Note, for clarity in the code implementation, the x_1 and x_3 coordinates are referred to as x and z, respectively.
    
   
 
    

 §.§ 

    The  class is used to define the material properties and thickness of each layer. The class is initialised by calling the class constructor with input arguments of the material name followed by its thickness. However, the thickness of the first and last layers are semi-infinite and their values should be set with the  keyword. The  class will automatically set the thickness of the first and last layer to  if another value is used. An example is given below.
        
    Here,  is an object array and every index in the object array corresponds to a different layer in the medium. In the current example,  will return a object with the material properties and thickness associated with PVDF. The  keyword can used for a material which is not predefined. The material properties and names can be set using their respective  functions. User defined materials can be added to the script . 

	

 §.§ Slowness profiles

	Slowness profiles are a plot of the inverse-phase velocity of each bulk wave component. They can be used to determine the angles of reflection and transmission between multi-layered media as well as the direction of energy propagation and skew angle <cit.>. Slowness profiles are found by solving the Christoffel equation, Eq. <ref>, and only depend on the material properties of each material. The method  is part of the  class, and calls the function 	
  	 
  	
  	which is an implementation of Eqs. (<ref>) and (<ref>). This takes input arguments of the material properties and phase-velocity and returns the polarisation and wave-vectors. The slowness profiles given by this function are plotted in terms of k_x/ω vs k_z/ω. For an isotropic material, the slowness profiles for each bulk wave are spherical, however, this is not true for an anisotropic material. An example of the slowness profiles for isotropic-glass and transverse-isotropic beryl is shown in Fig. <ref>. This figure has been reproduced from <cit.>. The slowness profiles of the (quasi-)longitudinal, (quasi-)shear-vertical and (quasi-)shear-horizontal bulk waves are shown. As glass is an isotropic material, the slowness profiles are spherical and the magnitudes of L,SV and SH when k_x/ω = 0 or k_z/ω = 0 are equal to the reciprocal of the compressional- and shear-speeds of glass. For the transverse-isotropic case, when k_x/ω = 0, the value of qL is equal to √(ρ/C_33) and qSV is equal to √(ρ/C_55). When k_z/ω = 0, the value of qL is equal to √(ρ/C_11) and the value of qSV is equal to √(ρ/C_55). These have been checked in the toolbox example script and all are within numerical precision.
	
	 
	 
      	
  	

 §.§ 

    The medium class is used to initialise the  class which runs the partial-wave method over a range of frequencies, wavenumbers, phasespeeds and angles. Two of these must be defined using the  functions. The  method is then used to run the partial-wave procedure. If the properties are not set before the  method is called, the model will run over a predefined range of frequencies and angles. 
    The  method constructs, rearranges and solves the global-matrix, Eq. (<ref>), using the function
  	
  	
  	
  	This function takes input arguments of the material properties and the parameters to calculate over (angles, frequencies, wavenumbers). It returns the determinant of the system matrix and the stresses and displacements at the layer interfaces. Each individual field-matrix is calculated using the function  
    
  	
  	
  	which is an implementation of Eq. (<ref>). This takes input arguments of the material properties, the wave-vector components, polarisation components and the phase velocity and returns the field-matrix. The default calculation is to find the partial-wave amplitudes and interface stresses and displacements when there is a single-frequency compressional wave incident on the structure from the first layer. An example is given below for a titanium plate.
    
    
    
	

 §.§ Reflection and Transmission Coefficients

	For a plane wave incident at an oblique angle on a multi-layered structure, the reflection and transmission coefficients determine the amplitude of the wave that is reflected and transmitted at each interface. Knowing these coefficients is useful for a number of applications. For example, selecting the appropriate launch angle when coupling energy into particular modes in a wave-guide, or determining the thickness and material properties of matching layers for ultrasonic transducers <cit.>. 
	
	The angles of refraction at the interfaces between multi-layered media can be found by studying the slowness profiles. However, slowness profiles do not take into account the boundary conditions at the interfaces. Consequently, the magnitude of each of the refracted waves cannot be calculated directly. For a plane wave incident on a multi-layered structure, the magnitude of the reflection and transmission coefficients are found by normalising the partial-wave amplitudes B_i^n by the incident plane wave amplitude B_1^1. This is automatically calculated when using the  method. 
	
	An example of the reflection and transmission coefficients at a PVDF-aluminium interface is given below and shown in Fig. <ref>. For a plane compressional wave incident on a PVDF-aluminium interface, there are four resulting refracted waves. These are a reflected R and transmitted T compressional L and shear S wave. The reflection and transmission coefficients have been compared to the analytic solutions for a two-layered elastic-medium from <cit.> and have an average error less than 1e^-15 which is within numerical precision for a 64 bit floating point number. For clarity the analytical solutions have not been plotted but can be seen in the toolbox examples folder.
	
	 
	 
    
	

 §.§ Dispersion Curves

	Dispersion curves describe the modal solutions of the multilayer structure and describe a wave-mode which propagates parallel to the layer-interfaces independently of a bulk wave. As one example, exciting these modes is essential in ultrasonic inspection. Knowledge of the dispersion curves is useful for determining the most appropriate modes to excite and for optimising the inspection process. 
	
	As mentioned in Section <ref>, the modal solutions are found when the global matrix becomes singular. The  software can calculate dispersion curves for simple layered structures (i.e., a plate in a vacuum or water). However, it is not robust for very-leaky cases, for example a plate embedded in soil. For these types of cases either Disperse, or other techniques based on the spectral-collocation method or semi-analytic finite element method are more appropriate <cit.>. An example of the dispersion curves for a 1 mm titanium plate in a vacuum is shown in Fig. <ref>(a). The dispersion curves are plotted on a graph of frequency vs wavenumber and show the first three symmetric S and anti-symmetric A Lamb modes. The results from Disperse are also plotted and have excellent agreement.
		
	 
	 
    


 §.§ Displacement and Stress Fields

    More information about the wave-physics and guided wave structures can be taken from dispersion curves by plotting the displacement and stress fields at different points. In the  software implementation, the x and z ranges over which to plot the displacement or stress fields must be specified. The  method returns a structure with the input ranges and field values at each point of the resulting grid. The values of the structure can be plotted independently or given as an argument to the  method. An example is given below for the displacement field within an titanium plate for a symmetric and anti-symmetric mode. The resulting plot can be seen in Fig. <ref>(b). 


    

 §.§ 

	One application of the toolbox is to model the directional response of Fabry-Pérot ultrasound sensors <cit.>.  The  is a child class of  and is an example of how the  toolbox may be expanded. The are additional inputs to this class which are described in more detail in the  folder, and a description of the modelling process can be found in <cit.>. The modelled directional response was found to have good agreement with the measured directional response. An example of the modelled and measured directional response for a glass etalon Fabry-Pérot sensor can be seen in Fig. <ref>. Features of the directional response correspond to symmetric and anti-symmetric Lamb modes propagating within the sensor. 
	
	 

    



§ IMPACT AND CONCLUSIONS


This paper introduces a new open-source toolbox called  which models elastic wave propagation in multi-layered media with anisotropic materials with isotropic or transverse-isotropic symmetry. The toolbox uses the partial-wave method which allows the calculation of slowness profiles, reflection and transmission coefficients, dispersion curves and stress and displacement fields. The software has been implemented using the object-oriented capabilities of MATLAB allowing for a simple command line or scripting interface. The implementation allows researchers to add functionality and integrate the software into other projects. For example, this toolbox has already been used to model the directional response of Fabry-Pérot ultrasound sensors <cit.>. It is anticipated the research user-base will actively contribute to  and add to the functionality. 



§ ACKNOWLEDGEMENTS


This work was supported in part by the Engineering and Physical Sciences Research Council (EPSRC), UK, Grant Nos. EP/P008860/1, EP/L020262/1 and EP/S026371/1, in part by the EPSRC-funded UCL Centre for Doctoral Training in Medical Imaging (EP/L016478/1) and the Department of Health’s NIHR-funded Biomedical Research Centre at NIHR Biomedical Research Centre at University College London Hospitals.

elsarticle-num 




§ CURRENT CODE VERSION





