diff --git a/src/test_christoffel_stepbystep.ipynb b/src/test_christoffel_stepbystep.ipynb index 074dfdb..f3b09b1 100644 --- a/src/test_christoffel_stepbystep.ipynb +++ b/src/test_christoffel_stepbystep.ipynb @@ -159,6 +159,16 @@ "Cijkl.shape" ] }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# normalize the stiffness tensor to density\n", + "Cijkl_norm = Cijkl /density" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -182,7 +192,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "metadata": {}, "outputs": [ { @@ -191,63 +201,17 @@ "(450, 3, 3)" ] }, - "execution_count": 7, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "Mil = ch._christoffel_matrix(wavevectors, Cijkl)\n", + "Mil = ch._christoffel_matrix(wavevectors, Cijkl_norm)\n", "\n", "Mil.shape" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Step 3: Normalize the Christoffel matrices\n", - "\n", - "$$\n", - "M_{il(norm)} = \\frac{M_{il}}{\\rho}\n", - "$$" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "norm_Mil = Mil / density" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(array([[ 59.4, 0. , 0. ],\n", - " [ 0. , 59.4, 0. ],\n", - " [ 0. , 0. , 118.8]]),\n", - " array([[11.19276427, 0. , 0. ],\n", - " [ 0. , 11.19276427, 0. ],\n", - " [ 0. , 0. , 22.38552855]]))" - ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# check\n", - "Mil[0], norm_Mil[0]" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -259,7 +223,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -268,13 +232,13 @@ "((450, 3), (450, 3, 3))" ] }, - "execution_count": 10, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "eigenvalues, eigenvectors = ch._calc_eigen(norm_Mil)\n", + "eigenvalues, eigenvectors = ch._calc_eigen(Mil)\n", "\n", "eigenvalues.shape, eigenvectors.shape" ] @@ -293,7 +257,7 @@ "V_p = \\sqrt{ \\frac{\\lambda_1}{\\rho} }, \\quad V_{s1} = \\sqrt{ \\frac{\\lambda_2}{\\rho} }, \\quad V_{s2} = \\sqrt{ \\frac{\\lambda_3}{\\rho} }\n", "$$\n", "\n", - "where $\\rho$ is the material density. If the Christoffel tensor has already been normalised with respect to the density, the general formulae are as follows:\n", + "where $\\rho$ is the material density. If the stiffness of the Christoffel tensor has already been normalised with respect to the density, the general formulae are as follows:\n", "\n", "$$\n", "v_p = \\sqrt{\\lambda_n}, \\quad \\lambda_n = v^2_p\n", @@ -304,7 +268,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": {}, "outputs": [ { @@ -313,7 +277,7 @@ "(450, 3)" ] }, - "execution_count": 11, + "execution_count": 10, "metadata": {}, "output_type": "execute_result" } @@ -326,7 +290,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "metadata": {}, "outputs": [ { @@ -344,7 +308,7 @@ " [3.02647863, 3.31846894, 4.95976291]])" ] }, - "execution_count": 12, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } @@ -356,7 +320,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": {}, "outputs": [ { @@ -410,7 +374,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 13, "metadata": {}, "outputs": [ { @@ -419,13 +383,13 @@ "(450, 3, 3, 3)" ] }, - "execution_count": 14, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "dMilk = ch._christoffel_gradient_matrix(wavevectors, Cijkl)\n", + "dMilk = ch._christoffel_gradient_matrix(wavevectors, Cijkl_norm)\n", "\n", "dMilk.shape" ] @@ -447,7 +411,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 14, "metadata": {}, "outputs": [ { @@ -456,7 +420,7 @@ "(450, 3, 3)" ] }, - "execution_count": 15, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -483,7 +447,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 15, "metadata": {}, "outputs": [ { @@ -492,7 +456,7 @@ "(450, 3, 3)" ] }, - "execution_count": 16, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -513,7 +477,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 16, "metadata": {}, "outputs": [ { @@ -522,7 +486,7 @@ "(450, 3)" ] }, - "execution_count": 17, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -535,25 +499,25 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[17.75488102, 17.75488102, 25.10919354],\n", - " [26.55988968, 17.75488102, 22.53753526],\n", - " [17.55089769, 20.07852322, 25.52255753],\n", - " [20.7581618 , 22.72141005, 23.46868979],\n", - " [27.76084197, 21.48549935, 21.54844796],\n", - " [32.82474072, 17.90867688, 20.14843909],\n", - " [27.53796823, 26.21150498, 19.78026365],\n", - " [18.69947453, 28.05020057, 22.31781173],\n", - " [17.93107493, 31.21136785, 21.28368726],\n", - " [19.24616973, 17.33475982, 26.77048673]])" + "array([[3.34555889, 3.34555889, 4.73133475],\n", + " [5.00468997, 3.34555889, 4.24675622],\n", + " [3.30712223, 3.78340366, 4.80922509],\n", + " [3.91146821, 4.28140382, 4.42221402],\n", + " [5.23098586, 4.0485207 , 4.06038213],\n", + " [6.1851782 , 3.3745387 , 3.79657793],\n", + " [5.18898968, 4.93904371, 3.7272025 ],\n", + " [3.523549 , 5.28550981, 4.20535363],\n", + " [3.37875917, 5.88116975, 4.01049317],\n", + " [3.62656298, 3.26639529, 5.04437285]])" ] }, - "execution_count": 18, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -579,7 +543,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 18, "metadata": {}, "outputs": [ { @@ -588,7 +552,7 @@ "(450, 3, 3)" ] }, - "execution_count": 19, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -617,7 +581,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 19, "metadata": {}, "outputs": [ { @@ -626,7 +590,7 @@ "(450, 3)" ] }, - "execution_count": 20, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } @@ -663,7 +627,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 20, "metadata": {}, "outputs": [], "source": [ @@ -679,7 +643,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 21, "metadata": {}, "outputs": [], "source": [ @@ -695,7 +659,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 22, "metadata": {}, "outputs": [], "source": [ @@ -711,7 +675,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 23, "metadata": {}, "outputs": [], "source": [ @@ -729,7 +693,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 24, "metadata": {}, "outputs": [ { @@ -771,9 +735,9 @@ " 4.731335\n", " 3.345559\n", " 3.345559\n", - " 25.109194\n", - " 17.754881\n", - " 17.754881\n", + " 4.731335\n", + " 3.345559\n", + " 3.345559\n", " \n", " \n", " 1\n", @@ -782,9 +746,9 @@ " 4.812982\n", " 3.345559\n", " 3.226995\n", - " 22.537535\n", - " 17.754881\n", - " 26.559890\n", + " 4.246756\n", + " 3.345559\n", + " 5.004690\n", " \n", " \n", " 2\n", @@ -793,9 +757,9 @@ " 4.836447\n", " 3.317070\n", " 3.221318\n", - " 25.522558\n", - " 20.078523\n", - " 17.550898\n", + " 4.809225\n", + " 3.783404\n", + " 3.307122\n", " \n", " \n", " 3\n", @@ -804,9 +768,9 @@ " 4.854540\n", " 3.344671\n", " 3.165071\n", - " 23.468690\n", - " 22.721410\n", - " 20.758162\n", + " 4.422214\n", + " 4.281404\n", + " 3.911468\n", " \n", " \n", " 4\n", @@ -815,9 +779,9 @@ " 4.877352\n", " 3.307152\n", " 3.169422\n", - " 21.548448\n", - " 21.485499\n", - " 27.760842\n", + " 4.060382\n", + " 4.048521\n", + " 5.230986\n", " \n", " \n", " ...\n", @@ -837,9 +801,9 @@ " 4.876995\n", " 3.311479\n", " 3.165452\n", - " 26.085941\n", - " 17.028408\n", - " 20.282967\n", + " 4.915384\n", + " 3.208669\n", + " 3.821927\n", " \n", " \n", " 446\n", @@ -848,9 +812,9 @@ " 4.854745\n", " 3.342711\n", " 3.166828\n", - " 26.198323\n", - " 17.960696\n", - " 18.433161\n", + " 4.936560\n", + " 3.384341\n", + " 3.473367\n", " \n", " \n", " 447\n", @@ -859,9 +823,9 @@ " 4.836353\n", " 3.318637\n", " 3.219843\n", - " 25.428975\n", - " 20.628592\n", - " 17.562545\n", + " 4.791591\n", + " 3.887053\n", + " 3.309317\n", " \n", " \n", " 448\n", @@ -870,9 +834,9 @@ " 4.813008\n", " 3.345199\n", " 3.227329\n", - " 24.088264\n", - " 21.434537\n", - " 19.575681\n", + " 4.538961\n", + " 4.038918\n", + " 3.688653\n", " \n", " \n", " 449\n", @@ -881,9 +845,9 @@ " 4.731335\n", " 3.345559\n", " 3.345559\n", - " 25.109194\n", - " 17.754881\n", - " 17.754881\n", + " 4.731335\n", + " 3.345559\n", + " 3.345559\n", " \n", " \n", "\n", @@ -891,36 +855,36 @@ "" ], "text/plain": [ - " polar_ang azimuthal_ang phase_Vp phase_Vs1 phase_Vs2 group_Vp \\\n", - "0 0.0 0.0 4.731335 3.345559 3.345559 25.109194 \n", - "1 9.8 0.0 4.812982 3.345559 3.226995 22.537535 \n", - "2 11.2 222.5 4.836447 3.317070 3.221318 25.522558 \n", - "3 12.4 85.0 4.854540 3.344671 3.165071 23.468690 \n", - "4 13.6 307.5 4.877352 3.307152 3.169422 21.548448 \n", - ".. ... ... ... ... ... ... \n", - "445 166.4 146.6 4.876995 3.311479 3.165452 26.085941 \n", - "446 167.6 9.0 4.854745 3.342711 3.166828 26.198323 \n", - "447 168.8 231.5 4.836353 3.318637 3.219843 25.428975 \n", - "448 170.2 94.0 4.813008 3.345199 3.227329 24.088264 \n", - "449 180.0 0.0 4.731335 3.345559 3.345559 25.109194 \n", + " polar_ang azimuthal_ang phase_Vp phase_Vs1 phase_Vs2 group_Vp \\\n", + "0 0.0 0.0 4.731335 3.345559 3.345559 4.731335 \n", + "1 9.8 0.0 4.812982 3.345559 3.226995 4.246756 \n", + "2 11.2 222.5 4.836447 3.317070 3.221318 4.809225 \n", + "3 12.4 85.0 4.854540 3.344671 3.165071 4.422214 \n", + "4 13.6 307.5 4.877352 3.307152 3.169422 4.060382 \n", + ".. ... ... ... ... ... ... \n", + "445 166.4 146.6 4.876995 3.311479 3.165452 4.915384 \n", + "446 167.6 9.0 4.854745 3.342711 3.166828 4.936560 \n", + "447 168.8 231.5 4.836353 3.318637 3.219843 4.791591 \n", + "448 170.2 94.0 4.813008 3.345199 3.227329 4.538961 \n", + "449 180.0 0.0 4.731335 3.345559 3.345559 4.731335 \n", "\n", " group_Vs1 group_Vs2 \n", - "0 17.754881 17.754881 \n", - "1 17.754881 26.559890 \n", - "2 20.078523 17.550898 \n", - "3 22.721410 20.758162 \n", - "4 21.485499 27.760842 \n", + "0 3.345559 3.345559 \n", + "1 3.345559 5.004690 \n", + "2 3.783404 3.307122 \n", + "3 4.281404 3.911468 \n", + "4 4.048521 5.230986 \n", ".. ... ... \n", - "445 17.028408 20.282967 \n", - "446 17.960696 18.433161 \n", - "447 20.628592 17.562545 \n", - "448 21.434537 19.575681 \n", - "449 17.754881 17.754881 \n", + "445 3.208669 3.821927 \n", + "446 3.384341 3.473367 \n", + "447 3.887053 3.309317 \n", + "448 4.038918 3.688653 \n", + "449 3.345559 3.345559 \n", "\n", "[450 rows x 8 columns]" ] }, - "execution_count": 25, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" } @@ -929,87 +893,102 @@ "data" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**I think I need to normalize the ``eigenvalue_gradients`` or ``group_velocities`` before calculating the group velocity magnitudes**" - ] - }, { "cell_type": "code", - "execution_count": 26, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "0 1.000000\n", - "1 0.882355\n", - "2 0.994372\n", - "3 0.910944\n", - "4 0.832497\n", - " ... \n", - "445 1.007871\n", - "446 1.016853\n", - "447 0.990745\n", - "448 0.943061\n", - "449 1.000000\n", - "Length: 450, dtype: float64" - ] - }, - "execution_count": 26, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "(data['group_Vp'] / density) / data['phase_Vp']" - ] - }, - { - "cell_type": "code", - "execution_count": 27, + "execution_count": 44, "metadata": {}, "outputs": [ { "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
polar_angazimuthal_angphase_Vpphase_Vs1phase_Vs2group_Vpgroup_Vs1group_Vs2
514.6170.04.8921643.3407853.1107783.7965783.3745396.185178
615.632.54.9140493.3024923.1171663.7272024.9390445.188990
\n", + "
" + ], "text/plain": [ - "0 4.731335\n", - "1 4.246756\n", - "2 4.809225\n", - "3 4.422214\n", - "4 4.060382\n", - " ... \n", - "445 4.915384\n", - "446 4.936560\n", - "447 4.791591\n", - "448 4.538961\n", - "449 4.731335\n", - "Name: group_Vp, Length: 450, dtype: float64" + " polar_ang azimuthal_ang phase_Vp phase_Vs1 phase_Vs2 group_Vp \\\n", + "5 14.6 170.0 4.892164 3.340785 3.110778 3.796578 \n", + "6 15.6 32.5 4.914049 3.302492 3.117166 3.727202 \n", + "\n", + " group_Vs1 group_Vs2 \n", + "5 3.374539 6.185178 \n", + "6 4.939044 5.188990 " ] }, - "execution_count": 27, + "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "data['group_Vp'] / density" + "data.iloc[5:7, :]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "Something is wrong here Vs2 > Vp. I think I skipped the step where the volocitoies are ordered. TO FOLLOW UP!\n", "\n", "---\n", + "\n", "# TEST" ] }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 27, "metadata": {}, "outputs": [], "source": [ @@ -1062,7 +1041,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 28, "metadata": {}, "outputs": [], "source": [ @@ -1078,18 +1057,18 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[-133.42561669, 21.32574585, 153.13073556],\n", - " [ -30.35018512, -10.78150255, 115.24183953],\n", - " [ 45.78019137, 10.32751876, 191.4713609 ]])" + "array([[-25.14143899, 4.01841829, 28.85448192],\n", + " [ -5.71889676, -2.03156257, 21.71506303],\n", + " [ 8.62637863, 1.94601823, 36.07902033]])" ] }, - "execution_count": 30, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } @@ -1101,18 +1080,18 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[-133.42561669, 21.32574585, 153.13073556],\n", - " [ -30.35018512, -10.78150255, 115.24183953],\n", - " [ 45.78019137, 10.32751876, 191.4713609 ]])" + "array([[-25.14143899, 4.01841829, 28.85448192],\n", + " [ -5.71889676, -2.03156257, 21.71506303],\n", + " [ 8.62637863, 1.94601823, 36.07902033]])" ] }, - "execution_count": 31, + "execution_count": 30, "metadata": {}, "output_type": "execute_result" } @@ -1132,18 +1111,18 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[-21.44569807, 3.42771889, 24.61293117],\n", - " [ -4.54237322, -1.61361812, 17.24771839],\n", - " [ 4.67893025, 1.05551634, 19.56918738]])" + "array([[-4.04102093, 0.64588636, 4.63782385],\n", + " [-0.85592109, -0.30405467, 3.24999404],\n", + " [ 0.88165258, 0.19889134, 3.68742932]])" ] }, - "execution_count": 32, + "execution_count": 31, "metadata": {}, "output_type": "execute_result" } @@ -1155,18 +1134,18 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([[-21.44569807, 3.42771889, 24.61293117],\n", - " [ -4.54237322, -1.61361812, 17.24771839],\n", - " [ 4.67893025, 1.05551634, 19.56918738]])" + "array([[-4.04102093, 0.64588636, 4.63782385],\n", + " [-0.85592109, -0.30405467, 3.24999404],\n", + " [ 0.88165258, 0.19889134, 3.68742932]])" ] }, - "execution_count": 33, + "execution_count": 32, "metadata": {}, "output_type": "execute_result" } @@ -1186,16 +1165,16 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([32.82474072, 17.90867688, 20.14843909])" + "array([6.1851782 , 3.3745387 , 3.79657793])" ] }, - "execution_count": 34, + "execution_count": 33, "metadata": {}, "output_type": "execute_result" } @@ -1207,16 +1186,16 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([32.82474072, 17.90867688, 20.14843909])" + "array([6.1851782 , 3.3745387 , 3.79657793])" ] }, - "execution_count": 35, + "execution_count": 34, "metadata": {}, "output_type": "execute_result" } @@ -1236,7 +1215,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 35, "metadata": {}, "outputs": [ { @@ -1247,7 +1226,7 @@ " [ 0.23222296, 0.052387 , 0.97125079]])" ] }, - "execution_count": 36, + "execution_count": 35, "metadata": {}, "output_type": "execute_result" } @@ -1259,7 +1238,7 @@ }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 36, "metadata": {}, "outputs": [ { @@ -1270,7 +1249,7 @@ " [ 0.23222296, 0.052387 , 0.97125079]])" ] }, - "execution_count": 37, + "execution_count": 36, "metadata": {}, "output_type": "execute_result" } @@ -1290,7 +1269,7 @@ }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 37, "metadata": {}, "outputs": [ { @@ -1299,7 +1278,7 @@ "array([0.46813036, 0.13431025, 0.48536834])" ] }, - "execution_count": 38, + "execution_count": 37, "metadata": {}, "output_type": "execute_result" } @@ -1311,7 +1290,7 @@ }, { "cell_type": "code", - "execution_count": 39, + "execution_count": 38, "metadata": {}, "outputs": [ { @@ -1320,7 +1299,7 @@ "array([2.41859939, 2.86906211, 2.90122615])" ] }, - "execution_count": 39, + "execution_count": 38, "metadata": {}, "output_type": "execute_result" } @@ -1339,7 +1318,7 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 39, "metadata": {}, "outputs": [ { @@ -1348,7 +1327,7 @@ "array([0.72299326, 0.27253055, 0.24036651])" ] }, - "execution_count": 40, + "execution_count": 39, "metadata": {}, "output_type": "execute_result" } @@ -1359,7 +1338,7 @@ }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 40, "metadata": {}, "outputs": [ { @@ -1368,7 +1347,7 @@ "array([2.98310074, 3.48292513, 0.22187522])" ] }, - "execution_count": 41, + "execution_count": 40, "metadata": {}, "output_type": "execute_result" }