diff --git a/examples/decomp.ipynb b/examples/decomp.ipynb new file mode 100644 index 0000000..04d99bf --- /dev/null +++ b/examples/decomp.ipynb @@ -0,0 +1,228 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Decomposing elastic constants\n", + "\n", + "This example provides an illustration of the implementation of the decomposition of elastic constants following the approach of Browaeys and Chevrot (2004). Stishovite elastic constants are used as for the example. The implementation closely follows that in MSAT (Walker and Wookey, 2012) and decomposition follows three stages. First, the input elastic constants are rotated such that the axes lie along high symmetry directions. Second, the matrix is decomposed into contributions for each symmetry. Third, the magnitude of each symmetry component is calculated. \n", + "\n", + "### References\n", + "\n", + "J. T. Browaeys and S. Chevrot. Decomposition of the elastic tensor and geophysical applications. Geophysical Journal International, 159:667 – 678, 2004.\n", + "\n", + "A. M. Walker and J. Wookey. MSAT – a new toolkit for the analysis of elastic and seismic anisotropy. Computers and Geosciences, 49:81 – 90, 2012." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Module imports etc\n", + "import sys\n", + "sys.path.append('..')\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# We need to import decompose and rotate\n", + "import pytasa.decompose\n", + "import pytasa.rotate\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false + }, + "source": [ + "## Setup the input\n", + "\n", + "Here we use the elastic constants of stishovite, but we rotate them in an arbitary way so that the symmetry is not obvious." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Stishovite elastic constants:\n", + "[[ 453. 211. 203. 0. 0. 0.]\n", + " [ 211. 453. 203. 0. 0. 0.]\n", + " [ 203. 203. 776. 0. 0. 0.]\n", + " [ 0. 0. 0. 252. 0. 0.]\n", + " [ 0. 0. 0. 0. 252. 0.]\n", + " [ 0. 0. 0. 0. 0. 302.]]\n", + "Rotated stishovite elastic constants:\n", + "[[ 635. 87.3 150.9 77.6 1.6 -20.6]\n", + " [ 87.3 679.3 200.6 31.8 -13.4 -14.9]\n", + " [ 150.9 200.6 724.2 35.2 -24.2 10.6]\n", + " [ 77.6 31.8 35.2 250.4 13.9 -8.6]\n", + " [ 1.6 -13.4 -24.2 13.9 213.2 58.3]\n", + " [ -20.6 -14.9 10.6 -8.6 58.3 164.1]]\n" + ] + } + ], + "source": [ + "stishovite_cij = np.array([[453.0, 211.0, 203.0, 0.0, 0.0, 0.0],\n", + " [211.0, 453.0, 203.0, 0.0, 0.0, 0.0],\n", + " [203.0, 203.0, 776.0, 0.0, 0.0, 0.0],\n", + " [ 0.0, 0.0, 0.0, 252.0, 0.0, 0.0],\n", + " [ 0.0, 0.0, 0.0, 0.0, 252.0, 0.0],\n", + " [ 0.0, 0.0, 0.0, 0.0, 0.0, 302.0]])\n", + "\n", + "print(\"Stishovite elastic constants:\")\n", + "print(np.array_str(stishovite_cij, precision=1, suppress_small=True))\n", + "\n", + "# Rotate around three axes\n", + "stishovite_cij = pytasa.rotate.rot_3(stishovite_cij, 20, 30, 40)\n", + "stishovite_rho = 4290.0\n", + "\n", + "print(\"Rotated stishovite elastic constants:\")\n", + "print(np.array_str(stishovite_cij, precision=1, suppress_small=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Find the high symmetry oritentation\n", + "This works by computing the Voigt stiffness tensor and dilational stiffness tensor from the input elastic constants tensor (represented as a matrix using Voigt notation). For high symmetry materials the eigenvectors allign. The number of distinct eigenvalues provides information about the symmetry. " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "FIXME: need to issue a warning\n", + "Elastic constants matrix alligned with symmetry directions:\n", + "[[ 592.2 71.8 203. -0. -0. -76.3]\n", + " [ 71.8 592.2 203. 0. -0. 76.3]\n", + " [ 203. 203. 776. 0. -0. -0. ]\n", + " [ -0. 0. 0. 252. 0. -0. ]\n", + " [ -0. -0. 0. -0. 252. -0. ]\n", + " [ -76.3 76.3 0. -0. -0. 162.8]]\n" + ] + } + ], + "source": [ + "cij_on_axes, rotation_matrix = pytasa.decompose.axes(stishovite_cij)\n", + "print(\"Elastic constants matrix alligned with symmetry directions:\")\n", + "print(np.array_str(cij_on_axes, precision=1, suppress_small=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# FIXME! - this is clearly wrong!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Decompose and calculate norms\n", + "Project the rotated matrix into isotropic, hexagonal, tetragonal, orthorhombic, monoclinic and triclinic components, then calculate the norms." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "decomposed = pytasa.decompose.decompose(cij_on_axes)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# FIXME: this is a bit rubbish\n", + "norms = pytasa.decompose.norms(cij_on_axes, decomposed[0], decomposed[1], decomposed[2], \n", + " decomposed[3], decomposed[4], decomposed[5])" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ElasticNorms(isotropic=0.76694217270623366, hexagonal=0.058050715699155098, tetragonal=0.027521562819149636, orthorhombic=0.0, monoclinic=0.14748554877546138, triclinic=1.1102230246251565e-16)\n" + ] + } + ], + "source": [ + "print(norms)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.4.5" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/pytasa/decompose.py b/pytasa/decompose.py new file mode 100644 index 0000000..2a21219 --- /dev/null +++ b/pytasa/decompose.py @@ -0,0 +1,475 @@ +# -*- coding: utf-8 -*- +""" +pytasa.decompose - decompose tensors following Browaeys & Chevrot (2004) + +This module provides functions that implement the decomposition of elasticity +tensors according to their symmetry using the approach described by Browaeys +and Chevrot (2004). The code is more-or-less a direct translation from the +Matlab implementation in MSAT (Walker and Wookey, 2012) + +References: + + Browaeys, J. T. and S. Chevrot (2004) Decomposition of the elastic tensor + and geophysical applications. Geophysical Journal international, + 159:667-678 + Walker, A. M. and Wookey, J. (2012) MSAT -- a new toolkit for the analysis + of elastic and seismic anisotropy. Computers and Geosciences, + 49:81-90. + +""" +from __future__ import division +import collections + +import numpy as np + +from . import rotate +from . import fundamental + +ElasticNorms = collections.namedtuple('ElasticNorms', ['isotropic', 'hexagonal', + 'tetragonal', 'orthorhombic', 'monoclinic', 'triclinic']) + + +def axes(c, x3_stiff=False): + + det_thresh = 100.0 * np.sqrt(np.spacing(1)) + + # Dilational stiffness tensor + dst = np.array([[c[0,0] + c[0,1] + c[0,2], c[0,5] + c[1,5] + c[2,5], + c[0,4] + c[1,4] + c[2,4]], + [c[0,5] + c[1,5] + c[2,5], c[0,1] + c[1,1] + c[2,1], + c[0,3] + c[1,3] + c[2,3]], + [c[0,4] + c[1,4] + c[2,4], c[0,3] + c[1,3] + c[2,3], + c[0,2] + c[1,2] + c[2,2]]]) + + # Voigt stiffness tensor + vst = np.array([[c[0,0] + c[5,5] + c[4,4], c[0,5] + c[1,5] + c[3,4], + c[0,4] + c[2,4] + c[3,5]], + [c[0,5] + c[1,5] + c[3,4], c[5,5] + c[1,1] + c[3,3], + c[1,3] + c[2,3] + c[4,5]], + [c[0,4] + c[2,4] + c[3,5], c[1,3] + c[2,3] + c[4,5], + c[4,4] + c[3,3] + c[2,2]]]) + + # Eigenvectors of these symmetric arrays + # these are dst_evecs[:,0], dst_evecs[:,1] etc. + dst_vals, dst_evecs = np.linalg.eigh(dst) + vst_vals, vst_evecs = np.linalg.eigh(vst) + + assert _isortho(dst_evecs[:,0], dst_evecs[:,1], dst_evecs[:,2]),\ + 'Eigenvectors of dilational stiffness tensor not orthogonal' + assert _isortho(vst_evecs[:,0], vst_evecs[:,1], vst_evecs[:,2]),\ + 'Eigenvectors of Voigt stiffness tensor not orthogonal' + + # Number and order of distinct vectors + dst_ndist, dst_ind = _ndistinct(dst_vals, thresh=np.linalg.norm(dst_vals) + / 1000.0) + vst_ndist, vst_ind = _ndistinct(vst_vals, thresh=np.linalg.norm(vst_vals) + / 1000.0) + # FIXME: dd a unittest for _ndistinct + + do_rot = False + mono_or_tri = False + # The symmetry depends on the number of distinct eigenvalues. These + # should be the same for the two tensors for high symmetry cases. + if dst_ndist == 1: + # Isotropic case + x1 = np.array([1, 0, 0]) + x2 = np.array([0, 1, 0]) + x3 = np.array([0, 0, 1]) + + elif dst_ndist == 2: + # Hexagonal or tetragonal + if x3_stiff: + # x3 should be the stiff direction, treat like ortho + inds = np.argsort(dst_vals) + if not(inds[0] == 0 and inds[1] == 1 and inds[2] == 2): + do_rot = True + x1 = dst_evecs[:, inds[0]] + x2 = dst_evecs[:, inds[1]] + x3 = dst_evecs[:, inds[2]] + + else: + # x3 should be the distinct direction + x3 = dst_evecs[:, dst_ind] + if not((x3[0] == 0) and (x3[1] == 0)): + do_rot = True + # NB: this fails when x3 = [0 1 0] - check MSAT! + x2 = np.cross(x3, np.array([x3[0], x3[1], 0])) + x2 = x2/np.linalg.norm(x2) + x1 = np.cross(x3,x2) + x1 = x1/np.linalg.norm(x1) + else: + # we won't need to rotate, but we do need x1 and x2 + x1 = np.array([1, 0, 0]) + x2 = np.array([0, 1, 0]) + + elif dst_ndist == 3: + # orthorhombic or lower symmetry + # If the two tensors are alligned we are ortho, otherwise lower sym + neq = 0 + for i in range(3): + for j in range(3): + neq = neq + _veceq(dst_evecs[:,i], vst_evecs[:,j], 0.01) + # FIXME: add a unittest for _veceq + + if neq == 3: + # Ortho significant axes are the three eigenvectors + do_rot = True + x1 = dst_evecs[:, 0] + x2 = dst_evecs[:, 1] + x3 = dst_evecs[:, 2] + + else: + if x3_stiff: + # x3 should be the stiff direction, should use d + inds = np.argsort(dst_vals) + if not(inds[0] == 0 and inds[1] == 1 and inds[2] == 2): + do_rot = True + x1 = dst_evecs[:, inds[0]] + x2 = dst_evecs[:, inds[1]] + x3 = dst_evecs[:, inds[2]] + + else: + # monoclinic or triclinic. Here we have to make a 'best-guess'. + # Following Browraeys and Chevrot we use the bisectrix of each + # of the eigenvectors of d and its closest match in v. + x = _estimate_basis(dst_evecs, vst_evecs) + x1 = x[:,0] + x2 = x[:,1] + x3 = x[:,2] + do_rot = True + mono_or_tri = True + + else: + # This should be impossible + raise ValueError("Number of distinct eigenvalues is impossible") + + # Now apply the necessary rotation. The three new axes define the rotation + # matrix which turns the 3x3 unit matrix (original axes) into the best + # projection. So the rotation matrix we need to apply to the elastic + # tensor is the inverse of this. + r1 = np.stack((x1, x2, x3), axis=1) + rr = r1.T + + # We don't want to flip the axes (and make a LH system). If this + # is what we would do avoid it by reflecting X1. + if (np.linalg.det(rr) + 1 < det_thresh): + print("FIXME: need to issue a warning") + r1 = np.stack((-x1, x2, x3), axis=1) + rr = r1.T + + cr = rotate.rot_c(c, rr) + + if mono_or_tri: + # if crystal is monoclinic or triclinic, we have some more work to do. + # We need to make sure that differently rotated crystals return the + # same result. So they are consistent we try 180 degree rotations + # around all the principle axes. For each flip, we choose the result + # that gives the fastest P-wave velocity in the [1 1 1] direction + cr, rr = _tryrotations(cr, rr) + + return cr, rr + + +def norms(c_ref, c_iso, c_hex, c_tet, c_ort, c_mon, c_tri): + + x_ref = _c2x(c_ref) + n = np.sqrt(np.dot(x_ref, x_ref)) + + c_tot = np.zeros((6,6)) + result = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) + for i, c in enumerate([c_iso, c_hex, c_tet, c_ort, c_mon, c_tri]): + c_tot = c_tot + c + x_tot = _c2x(c_tot) + x_d = x_ref - x_tot + result[i] = 1 - (np.sqrt(np.dot(x_d, x_d)) / n) + + # Rework the fractions + result[1:6] = result[1:6] - result[0:5] + + return ElasticNorms(result[0], result[1], result[2], result[3], result[4], + result[5]) + + +def decompose(c): + + c_work = np.copy(c) # so we don't stuff up input + out = [] + for i in range(5): + x = _c2x(c_work) + m = _get_projector(i) + xh = np.dot(m, x) + ch = _x2c(xh) + out.append(ch) + c_work = c_work - ch + out.append(c_work) # triclinic + return out + + +def _c2x(c): + """Convert elastic constants tensor (Voigt notation) into vector 'x' + """ + x = np.zeros(21) + + x[0] = c[0,0] + x[1] = c[1,1] + x[2] = c[2,2] + x[3] = np.sqrt(2.0) * c[1,2] + x[4] = np.sqrt(2.0) * c[0,2] + x[5] = np.sqrt(2.0) * c[0,1] + x[6] = 2.0 * c[3,3] + x[7] = 2.0 * c[4,4] + x[8] = 2.0 * c[5,5] + x[9] = 2.0 * c[0,3] + x[10] = 2.0 * c[1,4] + x[11] = 2.0 * c[2,5] + x[12] = 2.0 * c[2,3] + x[13] = 2.0 * c[0,4] + x[14] = 2.0 * c[1,5] + x[15] = 2.0 * c[1,3] + x[16] = 2.0 * c[2,4] + x[17] = 2.0 * c[0,5] + x[18] = 2.0 * np.sqrt(2.0) * c[4,5] + x[19] = 2.0 * np.sqrt(2.0) * c[3,5] ; + x[20] = 2.0 * np.sqrt(2.0) * c[3,4] ; + + return x + + +def _x2c(x): + """Convert vector 'x' into elastic constants tensor (Voigt notation) + """ + c = np.zeros((6,6)) + + c[0,0] = x[0] + c[1,1] = x[1] + c[2,2] = x[2] + c[1,2] = 1.0 / np.sqrt(2.0) * x[3] + c[0,2] = 1.0 / np.sqrt(2.0) * x[4] + c[0,1] = 1.0 / np.sqrt(2.0) * x[5] + c[3,3] = 1.0 / 2.0 * x[6] + c[4,4] = 1.0 / 2.0 * x[7] + c[5,5] = 1.0 / 2.0 * x[8] + c[0,3] = 1.0 / 2.0 * x[9] + c[1,4] = 1.0 / 2.0 * x[10] + c[2,5] = 1.0 / 2.0 * x[11] + c[2,3] = 1.0 / 2.0 * x[12] + c[0,4] = 1.0 / 2.0 * x[13] + c[1,5] = 1.0 / 2.0 * x[14] + c[1,3] = 1.0 / 2.0 * x[15] + c[2,4] = 1.0 / 2.0 * x[16] + c[0,5] = 1.0 / 2.0 * x[17] + c[4,5] = 1.0 / 2.0 * np.sqrt(2.0) * x[18] + c[3,5] = 1.0 / 2.0 * np.sqrt(2.0) * x[19] + c[3,4] = 1.0 / 2.0 * np.sqrt(2.0) * x[20] + + for i in range(6): + for j in range(6): + c[j,i] = c[i,j] + + return c + + +def _get_projector(order): + """Return the projector, M, for a given symmetry. + + Projector is a 21 by 21 numpy array for a given symmetry given by the + order argument thus: + order = 0 -> isotropic + order = 1 -> hexagonal + order = 2 -> tetragonal + order = 3 -> orthorhombic + order = 4 -> monoclinic + there is no projector for triclinic (it is just whatever is left once all + other components have been removed). + + NB: real division (like python 3) is used below, hence the __future__ + import at the top of the file. + """ + + srt2 = np.sqrt(2.0) + + # Isotropic + if order == 0: + m = np.zeros((21,21)) + m[0:9,0:9] = [ + [3/15, 3/15, 3/15, srt2/15, srt2/15, srt2/15, 2/15, 2/15, 2/15], + [3/15, 3/15, 3/15, srt2/15, srt2/15, srt2/15, 2/15, 2/15, 2/15], + [3/15, 3/15, 3/15, srt2/15, srt2/15, srt2/15, 2/15, 2/15, 2/15], + [srt2/15, srt2/15, srt2/15, 4/15, 4/15, 4/15, -srt2/15, -srt2/15, + -srt2/15], + [srt2/15, srt2/15, srt2/15, 4/15, 4/15, 4/15, -srt2/15, -srt2/15, + -srt2/15], + [srt2/15, srt2/15, srt2/15, 4/15, 4/15, 4/15, -srt2/15, -srt2/15, + -srt2/15], + [2/15, 2/15, 2/15, -srt2/15, -srt2/15, -srt2/15, 1/5, 1/5, 1/5], + [2/15, 2/15, 2/15, -srt2/15, -srt2/15, -srt2/15, 1/5, 1/5, 1/5], + [2/15, 2/15, 2/15, -srt2/15, -srt2/15, -srt2/15, 1/5, 1/5, 1/5]] + + # hexagonal + elif order == 1: + m = np.zeros((21,21)) + m[0:9,0:9] = [[3/8, 3/8, 0, 0, 0, 1/(4*srt2), 0, 0, 1/4], + [3/8, 3/8, 0, 0, 0, 1/(4*srt2), 0, 0, 1/4], + [0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1/2, 1/2, 0, 0, 0, 0], + [0, 0, 0, 1/2, 1/2, 0, 0, 0, 0], + [1/(4*srt2), 1/(4*srt2), 0, 0, 0, 3/4, 0, 0, -1/(2*srt2)], + [0, 0, 0, 0, 0, 0, 1/2, 1/2, 0], + [0, 0, 0, 0, 0, 0, 1/2, 1/2, 0], + [1/4, 1/4, 0, 0, 0, -1/(2*srt2), 0, 0, 1/2]] + + # tetragonal + elif order == 2: + m = np.zeros((21,21)) + m[0:9,0:9] = [[1/2, 1/2, 0, 0, 0, 0, 0, 0, 0], + [1/2, 1/2, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1/2, 1/2, 0, 0, 0, 0], + [0, 0, 0, 1/2, 1/2, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1/2, 1/2, 0], + [0, 0, 0, 0, 0, 0, 1/2, 1/2, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1]] + + # orthorhombic + elif order == 3: + m = np.zeros((21,21)) + for jj in range(9): + m[jj,jj] = 1 + + # monoclinic + elif order == 4: + m = np.eye(21) + for jj in [9, 10, 12, 13, 15, 16, 18, 19]: + m[jj,jj] = 0 + + else: + raise ValueError("Order must be 0, 1, 2, 3 or 4") + + return m + + +def _isortho(vec1, vec2, vec3, thresh=10.0*np.sqrt(np.spacing(1))): + "Confirms three vectors are mutually orthogonal, returns True or False" + dot_prods = np.absolute([np.dot(vec1, vec2), np.dot(vec1, vec3), + np.dot(vec2, vec3)]) + return not(np.any(dot_prods > thresh)) + + +def _ndistinct(vec, thresh=10.0*np.sqrt(np.spacing(1))): + """Returns the number and indicies of distinct values in a length-3 array + + Ignores differences of thresh. + """ + if ((np.abs(vec[0] - vec[1]) < thresh) and + (np.abs(vec[0] - vec[2]) < thresh)): + # All are the same + nd = 1 + i = None + elif ((np.abs(vec[0] - vec[1]) > thresh) and + (np.abs(vec[0] - vec[2]) > thresh) and + (np.abs(vec[1] - vec[2]) > thresh)): + # All are the different + nd = 3 + i = None + else: + # One is different + nd = 2 + if np.abs(vec[0] - vec[1]) < thresh: + i = 2 + elif np.abs(vec[0] - vec[2]) < thresh: + i = 1 + elif np.abs(vec[1] - vec[2]) < thresh: + i - 0 + else: + assert False, "Logic error in ndistinct!" + return nd, i + + +def _veceq(v1, v2, thresh): + """Return 1 if two 3D eigen-vectors are equal, otherwise 0 + + Ignores differences smaller than thresh and switch + sign of the vector if needed. + """ + # eignvector sign is arbitary, so, for each vector compare + # the sign of the largest element, and switch such that it is + # positive + i = np.absolute(v1).argmax() + if v1[i] < 0: + v1 = v1 * -1 + if v2[i] < 0: + v2 = v2 * -1 + + if (np.any((v1 - v2) > thresh)): + return 0 + else: + return 1 + + +def _estimate_basis(a, b): + """Estimate the best basis vectors for decomposition. + + This is the bisectrices of vectors a and b.""" + c = np.zeros((3,3)) + for j in [0, 1, 2]: + i = np.argmax(np.absolute([np.dot(a[:,j], b[:,0]), + np.dot(a[:,j], b[:,1]), + np.dot(a[:,j], b[:,2])])) + if np.dot(a[:,j], b[:,i]) < 0.0: + c[:,j] = _bisectrix(a[:,j], -b[:,i]) + else: + c[:,j] = _bisectrix(a[:,j], b[:,i]) + + # Enforce orthoganality, and renormalise + c[:,2] = np.cross(c[:,0], c[:,1]) + c[:,1] = np.cross(c[:,0], c[:,2]) + c[:,0] = c[:,0] / np.linalg.norm(c[:,0]) + c[:,1] = c[:,1] / np.linalg.norm(c[:,1]) + c[:,2] = c[:,2] / np.linalg.norm(c[:,2]) + + return c + + +def _bisectrix(a,b): + """return the unit length bisectrix of 3-vectors A and B""" + c = a + b + c = c / np.linalg.norm(c) + return c + + +def _tryrotations(cr, rr): + """try all combinations of 180 rotations around principle axes + select the one that gives the highest Vp in the [1 1 1] direction + """ + x1r = [000, 180, 000, 000, 180, 000] + x2r = [000, 000, 180, 000, 180, 180] + x3r = [000, 000, 000, 180, 000, 180] + vp = [] + for r1, r2, r3 in zip(x1r, x2r, x3r): + _, _, _, _, p = fundamental.phasevels(rotate.rot_3(cr, r1, r2, r3), + 2000, 45, 45) + vp.append(p[0]) + inds = np.argsort(vp) + + # Rotate to fastest direction + a = x1r[inds[-1]] + b = x2r[inds[-1]] + g = x3r[inds[-1]] + crr = rotate.rot_3(cr, a, b, g) + + # Add rotation to rr rotation matrix + a = a * np.pi/180.0 + b = b * np.pi/180.0 + g = g * np.pi/180.0 + r1 = np.array([[1.0, 0.0, 0.0], + [0.0, np.cos(a), np.sin(a)], + [0.0, -np.sin(a), np.cos(a)]]) + r2 = np.array([[np.cos(b), 0.0, -np.sin(b)], + [0.0, 1.0, 0.0], + [np.sin(b), 0.0, np.cos(b)]]) + r3 = np.array([[np.cos(g), np.sin(g), 0.0], + [-np.sin(g), np.cos(g), 0.0], + [0.0, 0.0, 1.0]]) + rrr = np.dot(r3, np.dot(r2, np.dot(r1, rr))) + return crr, rrr diff --git a/pytasa/rotate.py b/pytasa/rotate.py new file mode 100644 index 0000000..355b1a0 --- /dev/null +++ b/pytasa/rotate.py @@ -0,0 +1,85 @@ +# -*- coding: utf-8 -*- +""" +pytasa.rotate - rotations and related transforms for elasticity + +""" + +import numpy as np + +def rot_c(c, g): + """Rotate a elastic constants matrix given a rotation matrix + + This implementation makes use of the Bond transform, which avoids + expanding the 6x6 matrix to a 3x3x3x3 tensor. + + See Bowers 'Applied Mechanics of Solids', Chapter 3 + """ + k = np.empty((6,6)) + # k1 + k[0:3,0:3] = g[0:3,0:3]**2 + + # k2 + k[0:3,3] = g[0:3,1]*g[0:3,2]*2.0 + k[0:3,4] = g[0:3,2]*g[0:3,0]*2.0 + k[0:3,5] = g[0:3,0]*g[0:3,1]*2.0 + + # k3 + k[3,0:3] = g[1,0:3]*g[2,0:3] + k[4,0:3] = g[2,0:3]*g[0,0:3] + k[5,0:3] = g[0,0:3]*g[1,0:3] + + # k4 + k[3,3] = g[1,1]*g[2,2] + g[1,2]*g[2,1] + k[3,4] = g[1,2]*g[2,0] + g[1,0]*g[2,2] + k[3,5] = g[1,0]*g[2,1] + g[1,1]*g[2,0] + k[4,3] = g[2,1]*g[0,2] + g[2,2]*g[0,1] + k[4,4] = g[2,2]*g[0,0] + g[2,0]*g[0,2] + k[4,5] = g[2,0]*g[0,1] + g[2,1]*g[0,0] + k[5,3] = g[0,1]*g[1,2] + g[0,2]*g[1,1] + k[5,4] = g[0,2]*g[1,0] + g[0,0]*g[1,2] + k[5,5] = g[0,0]*g[1,1] + g[0,1]*g[1,0] + + cr = np.dot(np.dot(k, c), k.T) + + return cr + + +def rot_3(c, alpha, beta, gamma, order=None): + """Rotate a elastic constants matrix around the three axes. + + Rotate a 6x6 numpy array representing an elastic constants matrix + (c) by alpha degrees about the 1-axis, beta degrees about the 2-axis, + and gamma degress about the 3-axies. All rotations are clockwise when + looking at the origin from the positive axis. The order of the rotations + matter and the default is to rotate about the 1-axis first, the 2-axis + second and the 3-axis third. Passing a three-tuple of integers to the + optional "order" argument can be used to change this. For example: + + rot_3(c, 45, 90, 30, order=(2, 3, 1) + + will result in a rotation of 90 degrees about the 2-axis, 30 degrees + about the 3-axis then 45 degrees about the 1-axis. + """ + if order is None: + order = (1, 2, 3) + + alpha = np.radians(alpha) + beta = np.radians(beta) + gamma = np.radians(gamma) + + # Three rotation matrices - in a list so we can order them + # given "order" + r_list = [np.array([[1.0, 0.0, 0.0], + [0.0, np.cos(alpha), np.sin(alpha)], + [0.0, -1.0*np.sin(alpha), np.cos(alpha)]]), + np.array([[np.cos(beta), 0.0, -1.0*np.sin(beta)], + [0.0, 1.0, 0.0], + [np.sin(beta), 0.0, np.cos(beta)]]), + np.array([[np.cos(gamma), np.sin(gamma), 0.0], + [-1.0*np.sin(gamma), np.cos(gamma), 0.0], + [0.0, 0.0, 1.0]])] + + rot_matrix = np.dot(np.dot(r_list[order[2]-1], r_list[order[1]-1]), + r_list[order[0]-1]) + + return rot_c(c, rot_matrix) diff --git a/pytasa/validate.py b/pytasa/validate.py new file mode 100644 index 0000000..df4874c --- /dev/null +++ b/pytasa/validate.py @@ -0,0 +1,28 @@ +# -*- coding: utf-8 -*- +""" +pytasa.validate - tools for checking elaatic constants + +""" + +import numpy as np + +def cij_stability(cij): + """Check that the elastic constants matrix is positive + definite - i.e. that the structure is stable to small + strains. This is done by finding the eigenvalues by + diagonalization and checking that they are all positive. + See Born & Huang, "Dynamical Theory of Crystal Lattices" + (1954) page 141.""" + + stable = False + try: + L = np.linalg.cholesky(cij) + stable = True + + except np.linalg.LinAlgError: + # Swallow this exception and return False + print("Crystal not stable to small strains") + print("(Cij not positive definite)") + print("Matrix: " + str(cij)) + + return stable diff --git a/pytasa/vti.py b/pytasa/vti.py new file mode 100644 index 0000000..a491252 --- /dev/null +++ b/pytasa/vti.py @@ -0,0 +1,83 @@ +# -*- coding: utf-8 -*- +""" +pytasa.vti - tools for working with VTI symmetry elastic constants + +""" + +import numpy as np + +from .validate import cij_stability +from .rotate import rot_c + +def planar_ave(this_cij): + + # Find VTI elastic constants + if (np.amax(this_cij) <= 0.0) and (np.amin(this_cij) >= 0.0): + # no data at this point. Fill with NaN + xi = float('nan') + phi = float('nan') + vti_cij = None + + else: + # Find VRH mean by rotating around X3 axis + voigt_cij = np.zeros((6,6)) + reuss_sij = np.zeros((6,6)) + num_rots = 0.0 + for theta in range(-180, 180, 5): + theta = np.radians(theta) + g = np.array([[np.cos(theta), np.sin(theta), 0.0], + [-np.sin(theta), np.cos(theta), 0.0], + [0.0, 0.0, 1.0]]) + rot_cij = rot_c(this_cij, g) + assert cij_stability(rot_cij) + voigt_cij = voigt_cij + rot_cij + reuss_sij = reuss_sij + np.linalg.inv(rot_cij) + num_rots = num_rots + 1.0 + voigt_cij = voigt_cij / num_rots + reuss_cij = np.linalg.inv(reuss_sij / num_rots) + vrh_cij = (voigt_cij + reuss_cij) / 2.0 + + # Now build new matrix enforcing VTI symmetry + vti_cij = np.zeros((6,6)) + vti_cij[0,0] = (vrh_cij[0,0] + vrh_cij[1,1])/2.0 + vti_cij[1,1] = vti_cij[0,0] + vti_cij[2,2] = vrh_cij[2,2] + vti_cij[3,3] = (vrh_cij[3,3] + vrh_cij[4,4])/2.0 + vti_cij[4,4] = vti_cij[3,3] + vti_cij[5,5] = vrh_cij[5,5] + vti_cij[0,1] = vti_cij[0,0] - 2.0*vti_cij[5,5] + vti_cij[0,2] = (vrh_cij[0,2] + vrh_cij[1,2])/2.0 + vti_cij[1,2] = vti_cij[0,2] + # Lower half + vti_cij[1,0] = vti_cij[0,1] + vti_cij[2,1] = vti_cij[1,2] + vti_cij[2,0] = vti_cij[0,2] + + assert cij_stability(vti_cij) + + # Calculate velocities - for xi and phi the density does not matter + vph, vpv, vsh, vsv = vti_velocities(vti_cij, 20000.0) + # Calculate anisotropy measures - note that v and h swap for s and p + phi = (vpv**2/vph**2) + xi = (vsh**2/vsv**2) + + return xi, phi, vti_cij + + +def vti_velocities(c, r): + + c = c * r + + A = (3.0/8.0)*(c[0,0]+c[1,1])+0.25*c[0,1]+0.5*c[5,5] + C = c[2,2] + F = 0.5*(c[0,2]+c[1,2]) + N = (1.0/8.0)*(c[0,0]+c[1,1])-0.25*c[0,1]+0.5*c[5,5] + L = 0.5*(c[3,3]+c[4,4]) + + vph = np.sqrt(A/r) + vpv = np.sqrt(C/r) + + vsh = np.sqrt(N/r) + vsv = np.sqrt(L/r) + + return vph, vpv, vsh, vsv diff --git a/tests/rot_tests.py b/tests/rot_tests.py new file mode 100644 index 0000000..77f9fd1 --- /dev/null +++ b/tests/rot_tests.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Test the rotation routines + +These test cases are taken from MSAT +""" +import unittest +import numpy as np + +import pytasa.rotate + +olivine_cij_voigt = np.array([[320.5, 68.1, 71.6, 0.0, 0.0, 0.0], + [68.1, 196.5, 76.8, 0.0, 0.0, 0.0], + [71.6, 76.8, 233.5, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 64.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 77.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 78.7]]) + +albite_cij_voigt = np.array([[69.9, 34.0, 30.8, 5.1, -2.4, -0.9], + [34.0, 183.5, 5.5, -3.9, -7.7, -5.8], + [30.8, 5.5, 179.5, -8.7, 7.1, -9.8], + [ 5.1, -3.9, -8.7, 24.9, -2.4, -7.2], + [-2.4, -7.7, 7.1, -2.4, 26.8, 0.5], + [-0.9, -5.8, -9.8, -7.2, 0.5, 33.5]]) + +class PytasaRotTestCase(unittest.TestCase): + """ + Test case for rot functions + """ + def test_identity(self): + """Check that roation by the identity matrix does not rotate + + We should end up with what we started with. This is from MSAT + test_MS_rotR_simple. + """ + np.testing.assert_almost_equal(pytasa.rotate.rot_c( + olivine_cij_voigt, np.eye(3)), olivine_cij_voigt, decimal=6) + + def test_a_90(self): + """Check that roation around X axis works + + Result can be calculated by hand by permuting values. + This is from MSAT test_MS_rotR_simple. + """ + c_ol_rot_a_90 = np.array([[320.5, 71.6, 68.1, 0.0, 0.0, 0.0], + [71.6, 233.5, 76.8, 0.0, 0.0, 0.0], + [68.1, 76.8, 196.5, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 64.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 78.7, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 77.0]]) + g = np.array([[1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, -1.0, 0.0]]) + np.testing.assert_almost_equal(pytasa.rotate.rot_c( + olivine_cij_voigt, g), c_ol_rot_a_90, decimal=6) + + def test_b_90(self): + """Check that roation around Y axis works + + Result can be calculated by hand by permuting values. + This based on MSAT test_MS_rotR_simple. + """ + c_ol_rot_b_90 = np.array([[233.5, 76.8, 71.6, 0.0, 0.0, 0.0], + [76.8, 196.5, 68.1, 0.0, 0.0, 0.0], + [71.6, 68.1, 320.5, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 78.7, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 77.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 64.0]]) + g = np.array([[0.0, 0.0, 1.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0]]) + print(pytasa.rotate.rot_c(olivine_cij_voigt, g)) + print(c_ol_rot_b_90) + np.testing.assert_almost_equal(pytasa.rotate.rot_c( + olivine_cij_voigt, g), c_ol_rot_b_90, decimal=6) + + def test_c_90(self): + """Check that roation around Z axis works + + Result can be calculated by hand by permuting values. + This based on MSAT test_MS_rotR_simple. + """ + c_ol_rot_c_90 = np.array([[196.5, 68.1, 76.8, 0.0, 0.0, 0.0], + [68.1, 320.5, 71.6, 0.0, 0.0, 0.0], + [76.8, 71.6, 233.5, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 77.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 64.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 78.7]]) + g = np.array([[0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + print(pytasa.rotate.rot_c(olivine_cij_voigt, g)) + print(c_ol_rot_c_90) + np.testing.assert_almost_equal(pytasa.rotate.rot_c( + olivine_cij_voigt, g), c_ol_rot_c_90, decimal=6) + + + def test_rot_3_triclinic(self): + """Check that rotating albite forwards and backwards works + + This is based on MSAT test_MS_rot_3_triclinic + """ + c1 = pytasa.rotate.rot_3(albite_cij_voigt, 30.0, 45.0, 60.0) + c2 = pytasa.rotate.rot_3(c1, -30.0, -45.0, -60.0, order=(3, 2, 1)) + np.testing.assert_almost_equal(c2, albite_cij_voigt, decimal=6) + + +def suite(): + return unittest.makeSuite(PytasaRotTestCase, 'test') + + +if __name__ == '__main__': + unittest.main(defaultTest='suite') + diff --git a/tests/test_decompose.py b/tests/test_decompose.py new file mode 100644 index 0000000..0f2a979 --- /dev/null +++ b/tests/test_decompose.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Test the routines in pytasa.decompose + +Most of these tests confirm that our functions return the values given in the +paper by Browaeys and Chevrot (2004) and match test cases found in MSAT +(Walker and Wookey 2012). + +References: + + Browaeys, J. T. and S. Chevrot (2004) Decomposition of the elastic tensor + and geophysical applications. Geophysical Journal international, + 159:667-678. + Walker, A. M. and Wookey, J. (2012) MSAT -- a new toolkit for the analysis + of elastic and seismic anisotropy. Computers and Geosciences, + 49:81-90. +""" +import numpy as np +import pytasa.decompose +import pytasa.rotate + +# Define input and output from Browaeys and Chevrot (olivine, page 671) +ol_c_ref = np.array([[ 192.0, 66.0, 60.0, 0.0, 0.0, 0.0 ], + [ 66.0, 160.0, 56.0, 0.0, 0.0, 0.0 ], + [ 60.0, 56.0, 272.0, 0.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 60.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 62.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 49.0 ]]) +ol_c_iso = np.array([[ 194.7, 67.3, 67.3, 0.0, 0.0, 0.0 ], + [ 67.3, 194.7, 67.3, 0.0, 0.0, 0.0 ], + [ 67.3, 67.3, 194.7, 0.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 63.7, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 63.7, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 63.7 ]]) +ol_norm_iso = 0.793 +ol_c_hex = np.array([[ -21.7, 1.7, -9.3, 0.0, 0.0, 0.0 ], + [ 1.7, -21.7, -9.3, 0.0, 0.0, 0.0 ], + [ -9.3, -9.3, 77.3, 0.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, -2.7, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, -2.7, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 0.0, -11.7 ]]) +ol_norm_hex = 0.152 +ol_c_tet = np.array([[ 3.0, -3.0, 0.0, 0.0, 0.0, 0.0 ], + [ -3.0, 3.0, 0.0, 0.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 0.0, -3.0 ]]) +# NB: we do not know tet and ort from the paper (they are reported together) +ol_c_ort = np.array([[ 16.0, 0.0, 2.0, 0.0, 0.0, 0.0 ], + [ 0.0, -16.0, -2.0, 0.0, 0.0, 0.0 ], + [ 2.0, -2.0, 0.0, 0.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, -1.0, 0.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]]) +ol_norm_tet_and_ort = 0.055 +ol_c_mon = np.zeros((6,6)) +ol_norm_mon = 0.0 +ol_c_tri = np.zeros((6,6)) +ol_norm_tri = 0.0 + +albite_c = np.array([[ 69.9, 34.0, 30.8, 5.1, -2.4, -0.9], + [ 34.0, 183.5, 5.5, -3.9, -7.7, -5.8], + [ 30.8, 5.5, 179.5, -8.7, 7.1, -9.8], + [ 5.1, -3.9, -8.7, 24.9, -2.4, -7.2], + [ -2.4, -7.7, 7.1, -2.4, 26.8, 0.5], + [ -0.9, -5.8, -9.8, -7.2, 0.5, 33.5]]) + +def test_olivine_norms(): + """Check that we get the correct norms if we feed in the decomposed + values for olivine""" + norms = pytasa.decompose.norms(ol_c_ref, ol_c_iso, ol_c_hex, ol_c_tet, + ol_c_ort, ol_c_mon, ol_c_tri) + np.testing.assert_almost_equal(norms.isotropic, ol_norm_iso, decimal=3) + np.testing.assert_almost_equal(norms.hexagonal, ol_norm_hex, decimal=3) + np.testing.assert_almost_equal(norms.orthorhombic + norms.tetragonal, + ol_norm_tet_and_ort, decimal=3) + np.testing.assert_almost_equal(norms.monoclinic, ol_norm_mon, decimal=3) + np.testing.assert_almost_equal(norms.triclinic, ol_norm_tri, decimal=3) + +def test_olivine_decompose(): + """Check that we get the correct decomposed elasticity for olivine""" + decomp = pytasa.decompose.decompose(ol_c_ref) + np.testing.assert_allclose(decomp[0], ol_c_iso, atol=0.1) + np.testing.assert_allclose(decomp[1], ol_c_hex, atol=0.1) + np.testing.assert_allclose(decomp[2], ol_c_tet, atol=0.1) + np.testing.assert_allclose(decomp[3], ol_c_ort, atol=0.1) + np.testing.assert_allclose(decomp[4], ol_c_mon, atol=0.1) + np.testing.assert_allclose(decomp[5], ol_c_tri, atol=0.1) + +def test_olivine_axes(): + """Check that we get the correct rotation for olivine + + This checks that the olivine example from page 671 od B&C always comes + out with the correct axes however we rotate the input""" + + # Rotate the correct result on to the x3>x2>x1 reference frame + # this won't matter for future decomposition but we need it for the test + # cases. + c_ref = pytasa.rotate.rot_3(ol_c_ref, 0, 0, 90) + + for a, b, g in zip([0.0, 0.0, 45.0, 33.234, 123.0, 0.0, 0.0, 77.5, 266.0], + [0.0, 2.0, 2.0, 44.7, 266.0, 33.3, 0.0, 77.5, 0.0], + [0.0, 0.0, 0.0, 345.0, 0.0, 10.0, 10.0, 10.0, 10.0]): + c_test = pytasa.rotate.rot_3(c_ref, a, b, g) + np.testing.assert_allclose(pytasa.decompose.axes(c_test)[0], c_ref, + atol=0.000001) + +def test_triclinic_axes(): + """Check that different rotations for albite always give us the same axes""" + + c1, _ = pytasa.decompose.axes(albite_c) + + for a, b, g in zip([0.0, 0.0, 45.0, 33.234, 123.0, 0.0, 0.0, 77.5, 266.0], + [0.0, 2.0, 2.0, 44.7, 266.0, 33.3, 0.0, 77.5, 0.0], + [0.0, 0.0, 0.0, 345.0, 0.0, 10.0, 10.0, 10.0, 10.0]): + c2 = pytasa.rotate.rot_3(albite_c, a, b, g) + c2, _ = pytasa.decompose.axes(c2) + + np.testing.assert_allclose(c1, c2) + +def test_isortho(): + "Check that out orthogonal check routine works" + # These should be true + assert pytasa.decompose._isortho([1, 0, 0], [0, 1, 0], [0, 0, 1]) + assert pytasa.decompose._isortho([0, 1, 0], [1, 0, 0], [0, 0, 1]) + assert pytasa.decompose._isortho([1, 0, 0], [0, 0, 1], [0, 1, 0]) + assert pytasa.decompose._isortho([0, 0, 1], [0, 1, 0], [1, 0, 0]) + assert pytasa.decompose._isortho([1, 0, 0], [0, 0, 1], [0, 1, 0]) + assert pytasa.decompose._isortho([0, 1, 0], [0, 0, 1], [1, 0, 0]) + assert pytasa.decompose._isortho([-1, 0, 0], [0, 1, 0], [0, 0, 1]) + assert pytasa.decompose._isortho([0, -1, 0], [1, 0, 0], [0, 0, 1]) + assert pytasa.decompose._isortho([-1, 0, 0], [0, 0, 1], [0, 1, 0]) + assert pytasa.decompose._isortho([0, 0, -1], [0, 1, 0], [1, 0, 0]) + assert pytasa.decompose._isortho([-1, 0, 0], [0, 0, 1], [0, 1, 0]) + assert pytasa.decompose._isortho([0, -1, 0], [0, 0, 1], [1, 0, 0]) + # These should be false + assert not(pytasa.decompose._isortho([1, 0, 0], [0, 1, 0.01], [0, 0, 1])) + assert not(pytasa.decompose._isortho([1, 0, 0], [0.01, 1, 0], [0, 0, 1])) + assert not(pytasa.decompose._isortho([1, 0, 0], [0, 1, 0], [0.01, 0, 1])) diff --git a/tests/vti_tests.py b/tests/vti_tests.py new file mode 100644 index 0000000..df4ff73 --- /dev/null +++ b/tests/vti_tests.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Test the vti routines - based on results from my MATLAB implementation +""" +import unittest +import numpy as np + +import pytasa.vti + +class PytasaVTITestCase(unittest.TestCase): + """ + Test case for vti functions + """ + def test_example_point(self): + """Check we get the same results as those given in old files + + This test is for lat=-60, lon=85, r=3555 with cij from: + TX2008.V2.T9.6.topo.P001.geog_cij.dat and Xi from: + TX2008.V2.T9.6.topo.P001.Xi.dat + """ + ln_xi_percent_file = 24.319 + cij_file = np.array([[1091.91560284, 458.538622806, 456.496942875, + -1.04865250962, 1.84037540811, -38.9969923147], + [458.538622806, 1101.84460008, 452.074983665, + 2.75366574376, 0.5841747463, 24.6917896573], + [456.496942875, 452.074983665, 1224.90181259, + -4.09920876836, -5.74073877119, 6.43661390972], + [-1.04865250962, 2.75366574376, -4.09920876836, + 262.58650102, 1.908492257, 2.66401429285], + [1.84037540811, 0.5841747463, -5.74073877119, + 1.908492257, 264.1072491, 0.433616876376], + [-38.9969923147, 24.6917896573, 6.43661390972, + 2.66401429285, 0.433616876376, 356.605598817]]) + xi, phi, vti_cij = pytasa.vti.planar_ave(cij_file) + + np.testing.assert_almost_equal(np.log(xi)*100.0, ln_xi_percent_file, decimal=3) + + +def suite(): + return unittest.makeSuite(PytasaVTITestCase, 'test') + + +if __name__ == '__main__': + unittest.main(defaultTest='suite') +