anisotropy/data/analyze.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter(action = 'ignore')\n",
    "\n",
    "from anisotropy.database import Database, tables\n",
    "import pathlib\n",
    "\n",
    "import peewee as pw\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "db = Database(pathlib.Path(\"anisotropy.db\").resolve())\n",
    "execution = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data(execution: int, field: str):\n",
    "    if not db.getExecution(execution):\n",
    "        print(\"Execution not found\")\n",
    "\n",
    "    for model in db.tables:\n",
    "        try:\n",
    "            column = getattr(model, field)\n",
    "        \n",
    "        except AttributeError:\n",
    "            pass\n",
    "\n",
    "        else:\n",
    "            break\n",
    "\n",
    "    query = model.select(tables.Shape.alpha, column, tables.Shape.direction, tables.Shape.label)\n",
    "    idn = db.tables.index(model)\n",
    "\n",
    "    for table in reversed(db.tables[ :idn]):\n",
    "        query = query.join(table, pw.JOIN.LEFT_OUTER)\n",
    "        \n",
    "    query = query.switch(tables.Shape)\n",
    "    query = query.where(\n",
    "        tables.Shape.exec_id == execution,\n",
    "    #    tables.Shape.label == structure,\n",
    "    )\n",
    "    query = query.order_by(tables.Shape.label, tables.Shape.direction, tables.Shape.alpha)\n",
    "\n",
    "    with db:\n",
    "        if query.exists():\n",
    "            table = []\n",
    "            for row in query.dicts():\n",
    "                for k in row.keys():\n",
    "                    if type(row[k]) == list:\n",
    "                        row[k] = str(row[k])\n",
    "\n",
    "                table.append(row)\n",
    "        \n",
    "        else:\n",
    "            table = None\n",
    "\n",
    "    if table is None:\n",
    "        print(\"Results not found\")\n",
    "\n",
    "    else:\n",
    "        return pd.DataFrame(table)\n",
    "\n",
    "def permeability(label, direction):\n",
    "    def concrete(execution, field):\n",
    "        df = load_data(execution, field)\n",
    "        return df[df.label == label][df.direction == direction][field].to_numpy()\n",
    "\n",
    "    alpha = concrete(7, \"alpha\")\n",
    "    viscosity = concrete(5, \"viscosity\")\n",
    "    viscosityKinematic = concrete(5, \"viscosityKinematic\")\n",
    "    length = concrete(7, \"length\")\n",
    "    flowRate = concrete(5, \"flowRate\")\n",
    "    areaCellOutlet = concrete(7, \"areaCellOutlet\")\n",
    "    pressureInlet = concrete(5, \"pressureInlet\")\n",
    "    pressureOutlet = concrete(5, \"pressureOutlet\")\n",
    "    density = concrete(5, \"density\")\n",
    "\n",
    "    return viscosity * length * flowRate / (areaCellOutlet * (pressureInlet - pressureOutlet))\n",
    "\n",
    "def nanmean(arr):\n",
    "    temp = arr.copy()\n",
    "\n",
    "    if np.isnan(temp[0]):\n",
    "        temp[0] = temp[1]\n",
    "\n",
    "    for n, item in enumerate(temp):\n",
    "        if np.all(np.isnan(item)):\n",
    "            \n",
    "            vals = temp[n - 1 : n + 2]\n",
    "\n",
    "            if np.sum(~np.isnan(vals)) <= 1:\n",
    "                vals = temp[n - 2 : n + 3]\n",
    "\n",
    "            temp[n] = vals[~np.isnan(vals)].mean()\n",
    "\n",
    "    return temp\n",
    "\n",
    "def filter_group(arr, nan = True, qhigh = True, quantile = 0.97):\n",
    "    temp = arr.copy()\n",
    "    check = True\n",
    "    quan = np.quantile(temp[~np.isnan(temp)], quantile)\n",
    "    limit = 1000\n",
    "\n",
    "    while check:\n",
    "        if nan and np.any(np.isnan(temp)):\n",
    "            temp = nanmean(temp)\n",
    "            check = True\n",
    "        \n",
    "        elif qhigh and np.any(quan < temp):\n",
    "            temp[quan < temp] = np.nan\n",
    "            check = True\n",
    "\n",
    "        else:\n",
    "            check = False \n",
    "        \n",
    "        if limit <= 0:\n",
    "            break\n",
    "\n",
    "        else:\n",
    "            limit -= 1\n",
    "\n",
    "    return temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({\n",
    "    \"font.size\": 18,\n",
    "    \"lines.markersize\": 9.0,\n",
    "    \"lines.linewidth\": 3\n",
    "})\n",
    "savefig = True\n",
    "anisotropy = pd.DataFrame({\"simple\": None, \"bodyCentered\": None, \"faceCentered\": None}, index = [\"A21\", \"A31\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = load_data(execution, \"volume\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>alpha</th>\n",
       "      <th>volume</th>\n",
       "      <th>direction</th>\n",
       "      <th>label</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.010</td>\n",
       "      <td>6.299808e-15</td>\n",
       "      <td>[0.0, 0.0, 1.0]</td>\n",
       "      <td>bodyCentered</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.015</td>\n",
       "      <td>6.101743e-15</td>\n",
       "      <td>[0.0, 0.0, 1.0]</td>\n",
       "      <td>bodyCentered</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.020</td>\n",
       "      <td>5.902559e-15</td>\n",
       "      <td>[0.0, 0.0, 1.0]</td>\n",
       "      <td>bodyCentered</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.025</td>\n",
       "      <td>5.702339e-15</td>\n",
       "      <td>[0.0, 0.0, 1.0]</td>\n",
       "      <td>bodyCentered</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.030</td>\n",
       "      <td>5.508176e-15</td>\n",
       "      <td>[0.0, 0.0, 1.0]</td>\n",
       "      <td>bodyCentered</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>316</th>\n",
       "      <td>0.250</td>\n",
       "      <td>1.841166e-15</td>\n",
       "      <td>[1.0, 1.0, 1.0]</td>\n",
       "      <td>simple</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>317</th>\n",
       "      <td>0.255</td>\n",
       "      <td>1.705807e-15</td>\n",
       "      <td>[1.0, 1.0, 1.0]</td>\n",
       "      <td>simple</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>318</th>\n",
       "      <td>0.260</td>\n",
       "      <td>1.574410e-15</td>\n",
       "      <td>[1.0, 1.0, 1.0]</td>\n",
       "      <td>simple</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>319</th>\n",
       "      <td>0.265</td>\n",
       "      <td>1.447357e-15</td>\n",
       "      <td>[1.0, 1.0, 1.0]</td>\n",
       "      <td>simple</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>320</th>\n",
       "      <td>0.270</td>\n",
       "      <td>1.325051e-15</td>\n",
       "      <td>[1.0, 1.0, 1.0]</td>\n",
       "      <td>simple</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>321 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     alpha        volume        direction         label\n",
       "0    0.010  6.299808e-15  [0.0, 0.0, 1.0]  bodyCentered\n",
       "1    0.015  6.101743e-15  [0.0, 0.0, 1.0]  bodyCentered\n",
       "2    0.020  5.902559e-15  [0.0, 0.0, 1.0]  bodyCentered\n",
       "3    0.025  5.702339e-15  [0.0, 0.0, 1.0]  bodyCentered\n",
       "4    0.030  5.508176e-15  [0.0, 0.0, 1.0]  bodyCentered\n",
       "..     ...           ...              ...           ...\n",
       "316  0.250  1.841166e-15  [1.0, 1.0, 1.0]        simple\n",
       "317  0.255  1.705807e-15  [1.0, 1.0, 1.0]        simple\n",
       "318  0.260  1.574410e-15  [1.0, 1.0, 1.0]        simple\n",
       "319  0.265  1.447357e-15  [1.0, 1.0, 1.0]        simple\n",
       "320  0.270  1.325051e-15  [1.0, 1.0, 1.0]        simple\n",
       "\n",
       "[321 rows x 4 columns]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Porosity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABBUAAAIhCAYAAAARuoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADDnUlEQVR4nOzdd3hU1dbH8e9Meif0FiD0jjSVTmiCwlVABARpAiqK135BvBh87YoFG4JKERBBpIiIYOhFFAhSpJMECFID6XVm3j9yicTMhBDClOT3eZ48zjl7nXNWYBsya3YxWCwWCyIiIiIiIiIiN8jo6ARERERERERExDWpqCAiIiIiIiIihaKigoiIiIiIiIgUiooKIiIiIiIiIlIoKiqIiIiIiIiISKGoqCAiIiIiIiIihaKigoiIiIiIiIgUiooKIiIiIiIiIlIo7o5OQPJnNps5c+YMAQEBGAwGR6cjIiIiIiIixZzFYiExMZHKlStjNOY/FkFFBSd35swZQkJCHJ2GiIiIiIiIlDCnTp2iatWq+caoqODkAgICgOy/zMDAQJtxmZmZrFmzhh49euDh4WGv9ERyUT8UZ6G+KM5CfVGcgfqhOAv1RdeRkJBASEhIzvvR/Kio4OSuTnkIDAy8blHB19eXwMBA/Q8qDqN+KM5CfVGchfqiOAP1Q3EW6ouupyBT8LVQo4iIiIiIiIgUiooKIiIiIiIiIlIoKiqIiIiIiIiISKGoqCAiIiIiIiIihaKigoiIiIiIiIgUiooKIiIiIiIiIlIoKiqIiIiIiIiISKGoqCAiIiIiIiIiheLu6ARuVGJiIlOnTmXJkiVERUXh5uZG3bp1GTRoEOPHj8fT07PInvXoo4/y+eefA1C9enWio6PzjT937hxvv/02K1eu5OTJk/j4+NCoUSOGDx/Oww8/jMFgKLLcRERERESkZLFYLGRmZmI2mx2dSqFkZmbi7u5OWloaJpPJ0ekUe0ajEQ8Pj1v+PtSligoxMTF07tw55829r68v6enp7Ny5k507dzJ//nwiIiIIDg6+6WetX7+eGTNmFDh+165d3HXXXVy6dAkAf39/EhMT2bJlC1u2bOG7775jxYoVRVr0EBERERGR4s9kMnHx4kUSExPJzMx0dDqFZrFYqFixIqdOndIHrnbi4eFBQEAAZcuWxc3N7ZY8w2WKCllZWfTp04fo6GgqVarE3Llz6datG2azmcWLFzNmzBgiIyMZOnQoP/744009KyUlhTFjxuDu7k6zZs3YuXNnvvHx8fH07t2bS5cuUb9+fb7++mtatWpFRkYGM2fO5Omnn+bnn3/mqaee4tNPP72p3EREREREpOQwmUycOnWK9PR0goKC8Pf3x83NzSXflJvNZpKSkvD398do1Ez8W8lisWAymUhKSuLKlSukpqYSEhJySwoLLlNUmDNnDvv27QNgyZIltGnTBsge0jFw4EDMZjMPPvggq1atIiIigq5duxb6WZMmTeL48eNMmjSJ06dPX7eo8O6773L27Fl8fHxYtWoVoaGhAHh6evL444+TkJDAiy++yIwZM3jqqaeoW7duoXMTEREREZGS4+LFi6Snp1OtWjV8fHwcnc5NMZvNZGRk4O3traKCnfj7+xMUFMTJkye5ePEiFSpUKPJnuMzf5Jw5cwAICwvLKShca9CgQTlv5ufOnVvo5/z6669MmzaNunXr8tJLLxXomqvPuzaHa40fPx5/f39MJhPz588vdG4iIiIiIlJyWCwWEhMTCQoKcvmCgjiOj48PgYGBJCYmYrFYivz+LlFUSElJYevWrQD06tXLaozBYKBnz54ArFmzplDPSU9PZ9SoUVgsFmbMmIG3t/d1rzl8+DAnT57MNzd/f386dOhwU7mJiIiIiEjJkpmZSWZmJv7+/o5ORVxcQEBATn8qai5RVDh48GDOCqeNGze2GXe17ezZs8TFxd3wc1555RUOHjzIww8/TKdOnQp0zf79+/M8P7/c/vzzzxvOS0RERERESp6r74Fu1QJ7UnJc7UO3YucQl1hT4cyZMzmvq1SpYjPu2rYzZ85QunTpAj8jMjKSt99+mwoVKvDOO+/cstwSEhJyFiexJj09nfT09JzjhIQEgOtWla62OXo12MS0TAK8PRyagziOs/RDEfVFcRbqi+IM1A9dV2ZmJhaLBYvF4rLbSF7r6tD74vL9uJKr/SgzM7NARaob+XnhEkWFxMTEnNe+vr42465tu/aa68nKymLUqFFkZWUxbdo0SpUqdctzs1VUeOONN5gyZUqe82vWrMn3/letXbv2ujG3yqkk+HC/G83LWuhQ0Uw1jdIqsRzZD0Wupb4ozkJ9UZyB+qHrcXd3p2LFiiQlJZGRkeHodIrMjbxXk6KRkZFBamoqmzZtIisr67rxKSkpBb63SxQVbrU333yTPXv20Lt3bx544AGH5jJx4kSeeeaZnOOEhARCQkLo0aMHgYGBNq/LzMxk7dq1dO/eHQ8Px4wUeHHZATItsfx2wcBvF4x0b1CeTx+8zSG5iGM4Qz8UAfVFcR7qi+IM1A9dV1paGqdOncLf379A6705u6sLTwYEBLjklpiuLC0tDR8fHzp27FigvnR1xHxBuERRISAgIOd1fhWTa9uuvSY/f/75J//3f/+Hv78/n3766U3nZuuNf0Fz8/LywsvLK895Dw+PAv0jUNC4ohafkskPe//Kda5J1VL6h6uEclQ/FPkn9UVxFuqL4gzUD12PyWTCYDBgNBqLxRaMV6c8XP2exH6MRiMGg+GG3lcW+N43k5i9VK5cOed1bGyszbhr2669Jj+PP/44GRkZTJo0ieDgYJKSknJ9XR0aYrFYcs5dO7/kRnMLDAwslqu37j55GZP57+1J3I0GBrUOcWBGIiIiIiIicqu5RFGhQYMGOZWsa3db+KerbRUrVizwIo1RUVFA9rSDgICAPF/z588H4OTJkznnPvnkk5zrr93xoSC5NWzYsEB5uZqw+uXZNqErz99Vj8pB3vRsXJHygdaH1WSZzDzw+XY+33icy8nFZ26YiIiIiIg4r/DwcAwGw3WnXmRlZfHAAw/kxL7xxhs5bSNGjMg5fyNfGzZsuMXfneO4RFHB19eXdu3aAbB69WqrMRaLhZ9//hmAHj162C23unXrUq1atXxzS05OZvPmzXbPzd7KBXjxeFhtNr0Qxv/da3t7zV8Onue3qDje+OkQd7wRwbOL/iAt02THTEVERERERPIymUwMGTKExYsXA/Dqq68yceLEPHFGo5EKFSrk+xUcHGzv9B3CJYoKAMOHDwdg/fr17NixI0/74sWLOXHiBADDhg0r8H2jo6Nzttew9nX1udWrV88599RTT+VcbzAYcp63cOFCoqOj8zzjk08+ISkpCTc3N4YMGVLg3FyVu5uRYD9Pm+3zfo3JeZ2RZebYhSS8PbT3roiIiIiIOI7JZGLo0KEsWrQIgClTpjBp0iSrsSEhIZw9ezbfr++//96e6TuMSxUVmjRpgsVioX///kRERADZi30sXryYMWPGANCrVy+6du2a69prh7lYe9N/s5577jkqVqxISkoK99xzD7t27QKyt+347LPP+O9//wvA2LFjqVu3bpE/35WcvJTClmMXc50bdmd1B2UjIiIiIiKSXVAYNmwYCxcuBGDy5MlMnjzZwVm5BpfY/QGy92hdsWIFYWFhREdH061bN3x9fTGbzaSlpQHQvHnznDUQ7CkoKIiVK1dy11138eeff9KqVSsCAgJIS0vLWdSxR48evP/++3bPzdmElPZhyWNtmLs9hlX7/sLfy517mlayGf/97tPUKudPs5BS9ktSRERERERKDLPZzIgRI1iwYAEAkyZNYsqUKQ7OynW4TFEBoEaNGuzdu5d3332X77//nqioKDw8PGjUqBGDBw9m/PjxeHraHnZ/K7Vs2ZIDBw7w1ltvsXLlSk6dOoWfnx+NGzdm+PDhjBo1StumkD1dpGX10rSsXpr/9m7IkbOJNqc+JKRl8t9l+0nOMNGyejAj29WgZ6OKuLvpz1FERERERG6e2Wxm5MiRzJs3D4AJEybw6quvOjgr1+JSRQWAgIAApky
      "text/plain": [
       "<Figure size 1200x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = load_data(execution, \"porosity\")\n",
    "alpha = df[\"alpha\"].unique()\n",
    "simple = df[df.label == \"simple\"][df.direction == \"[0.0, 0.0, 1.0]\"][\"porosity\"].to_numpy()\n",
    "bodyCentered = df[df.label == \"bodyCentered\"][df.direction == \"[0.0, 0.0, 1.0]\"][\"porosity\"].to_numpy()\n",
    "faceCentered = df[df.label == \"faceCentered\"][df.direction == \"[0.0, 0.0, 1.0]\"][\"porosity\"].to_numpy()\n",
    "\n",
    "fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize = (12, 6))\n",
    "\n",
    "ax.plot(alpha, np.pad(simple, (0, alpha.size - simple.size), 'constant', constant_values = np.nan), \":\", label = \"КП\")\n",
    "ax.plot(alpha, np.pad(bodyCentered, (0, alpha.size - bodyCentered.size), 'constant', constant_values = np.nan), \"--\", label = \"КОП\")\n",
    "ax.plot(alpha, np.pad(faceCentered, (0, alpha.size - faceCentered.size), 'constant', constant_values = np.nan), \"-.\", label = \"КГЦ\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$m$\")\n",
    "plt.show()\n",
    "\n",
    "if savefig:\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(\"porosity-rounded.tiff\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+cAAAIhCAYAAAAsOMuhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABkAElEQVR4nO3deXxU1f3/8fdMMmQPi6BhiQRRsApWxKUWQRFBoVoR1yIV1J9aqn5dAHcltAgVgn6/LdraUlErahVbES2WyFLqghShVZAHVIEAASlhyUJCMsnc3x9pxkAyySx35p6ZvJ6PRx4kM3PPnDvzYZL3Peee67IsyxIAAAAAAHCM2+kOAAAAAADQ1hHOAQAAAABwGOEcAAAAAACHEc4BAAAAAHAY4RwAAAAAAIcRzgEAAAAAcBjhHAAAAAAAhxHOAQAAAABwWLLTHYgln8+n3bt3KysrSy6Xy+nuAAAAAAASnGVZKi8vV7du3eR2Bx4fb1PhfPfu3crNzXW6GwAAAACANmbnzp3q0aNHwPvbVDjPysqSVP+iZGdnN/sYr9erpUuXasSIEfJ4PLHsHnAUahGmoBZhAuoQpqAWYQLqML6UlZUpNzfXn0cDaVPhvGEqe3Z2dovhPD09XdnZ2RQ6HEUtwhTUIkxAHcIU1CJMQB3Gp9ZOrWZBOAAAAAAAHEY4BwAAAADAYYRzAAAAAAAcRjgHAAAAAMBhhHMAAAAAABxGOAcAAAAAwGGEcwAAAAAAHEY4BwAAAADAYYRzAAAAAAAcRjgHAAAAAMBhhHMAAAAAAByW7HQHAAAAkHiKD1Xp4OGaoB/fMaOdundIi2KPAMBshHMAAADYqvhQlS4uWKnqWl/Q26Qku7V88kUEdABtFtPaAQAAYKuDh2tCCuaSVF3rC2mkHQASDeEcAAAAAACHEc4BAAAAAHAY4RwAAAAAAIcRzgEAAAAAcFhYq7Xv379f77zzjpYtW6Z169apqKhItbW16tKli84++2yNHz9eV111VVgdWrdunRYvXqzPPvtMW7Zs0b59+1RWVqbs7GydeuqpGjVqlCZOnKhOnTqF1T4AAABa5/Sl0Jx+fjskwj4AiJ2wwnlOTo5qa2v9P6empsrj8ai4uFjFxcVatGiRRo4cqYULFyo9PT2ktl944QU9++yzR7WdlpamAwcO6OOPP9bHH3+s//3f/9U777yj888/P5zuAwAAoAVOXwrN6ee3QyLsA4DYCmtae21trc4991w999xz+vrrr1VVVaWKigpt27ZNt956qyRpyZIluuOOO0Ju+9xzz9Xs2bP1ySef6ODBg6qqqlJZWZnKy8v10ksvqUuXLiopKdHo0aNVWloaTvcBAADQAqcvheb089shEfbBBMWHqrShuDTor+JDVU53GQhbWCPny5cv19ChQ5vcnpeXp3nz5ik5OVnPP/+8XnnlFc2YMUO5ublBt33TTTc1e3tmZqZuuukm5eTk6NJLL9V//vMfvfvuu7rxxhvD2QUAAACgRc1NS6+trdXOCmnj7jIlJx/9pzTT0u3F7AO0NWGF8+aCeWO33nqrnn/+eUnS2rVrQwrnrfne977n/37Xrl22tQsAAOpxniwSRSS13HIwTFbBF6ub3EowtFcksw94DxCPwgrnrUlNTfV/X1dXZ2vbf//73/3f9+7d29a2AQBo6xipQqKItJYJhomDA46IF1EJ5ytXrvR/379//4jbq66u1p49e/Tuu+/qiSeekCSdfPLJuuKKKyJuGwAAfItAAjt0zGinlGR3yMG4Y0Y72/pALUPigCPii+3h/NChQ5o5c6YkafDgwerbt2/YbaWmpqq6urrJ7YMGDdKrr76qlJSUsNsGAABAdHTvkKblky9itBKOs+MgTaQj74zcI1i2hnOfz6cf//jH2rNnj1JTUzV37tyI2svJydGRI0dUUVGhw4cPS6o/333WrFk68cQTW92+urr6qHBfVlYmSfJ6vfJ6vc1u03B7oPuBWKEWYQpqsW1pfKnUULeLZo1Qh7FlRx0cn5Gs4zNC+1OzYVs7nj/SNpzeHmbUwe5DVRrxfx+FPPK+9J5B6tYhLeLtJWn3oSodrDy6JhoWJvzXjgNNFyZM9/i3hRmC/T9tazi/55579O6770qSnn32WZ1xxhkRtbd9+3b/9//5z3/0hz/8QU8++aTOPfdcPfbYY/rZz37W4vYzZ87UtGnTmty+dOnSVq+/XlhYGFafAbtRizAFtdg27KyQwvnz4MMPP1RRpu3daYI6jA2n68CO54+0Dae3hzl1UF0b2vbVtT4tXrpCuTZsf6BaenJ9kmotVzOPTJa+WNv0VpelRwfUqROTjI1RWVkZ1ONclmVZdjzh5MmTNWfOHEnSM888o3vvvdeOZptYs2aNzj//fPl8Pi1evFiXX355wMc2N3Kem5urkpISZWdnN7uN1+tVYWGhhg8fLo/HY3v/gWBRizAFtdi2bNxdptG/broKdWvenvg9nd6t+d+tdqAOQ9fcaFtLGo+2OV0Hdjx/pG04vT2oA+oocZSVlalz584qLS0NmEMlm0bOH3jgAX8wLygoiFowl6Rzzz1XF1xwgVatWqXf/va3LYbzlJSUZs9L93g8rf5iD+YxQCxQi/Ej0c8poxbbhmOnR4ayXSzqI1Z1GO//n4vDnErbsAiW03Vgx/NH2obT29slnmuZOjCnjhCZYN+LiMP5lClTVFBQIEmaNWuWJk2aFGmTrerevbsk6auvvor6cwFAMBJlNdjm/ohrOK9t4+6ypue1GfRHHGCXRPj/zErlkBKjloG2JKJw3ngq+6xZszRlyhRbOtWarVu3SpKysrJi8nwA0BoTVoONVMt/xCWr4Ium0+pM+yPO6dcQiYFg6/yl0Jx+fjuYsA/UMhBfwg7njYN5QUGBLSPmdXV1crvdcrmaW/Cg3rJly7RmzRpJ0kUXXRTxcwJIDPEeykwY3Yj3AwwmvIZ2iPdaRmJw+lJoTj+/HRJhHwDEVljhvPE55k8//bTuu+++oLd98cUXdfPNN0uSVqxYcVTA3rlzp0aPHq2JEydq+PDh6tWrlz+o79y5UwsWLND06dNlWZY6deoU0vMCSFx2hTInQ1EijG7Y8T5E8h7wGpqxDyZIhAMcJuxD9w5pjr4uTj+/HRJhH5ysRRNmHwCxFHI437Fjh2bPni1Jcrvdeuqpp/TUU08FfPzkyZM1efLkoNv/17/+pZ/85CeSpHbt2ik7O1tVVVX+65xLUq9evfTWW28pJycn1O4DSEB2jfgSiiIT6fvAe5AYBxjsEMnaB4lQR4mwD4nAhGDo9EEap2uR2Qdoa0IO5z6f76jv9+7d2+LjKyoqgm67W7duevPNN7Vy5Up9+umn2r17t0pKSpSUlKQTTzxR3/3ud3XllVdq7NixSkvjPx0A+xCKnMd7YAanA0mkax8kQh0lwj6YINJaDhQMa2tr9eGHH+qCCy6I6iKZTgdjyYxaTITZB0CwQg7neXl5iuTS6BMmTNCECROava9du3a65pprdM0114TdPgC0VU6PsMAMkdaBHSNVbf30BJjBjlpuLhh6vV4VZUqnd8uO6qWqEuX/gtO/m5w+4AiEwpbrnAMAnGXCCEuicPoPyUjYVQeRjFRRizAJo67OMuHzgKnxiCeEcwC2iOdAkwgSZYTFaSb8IRkJE+rAhD4AMIMpnweRHKSJdOSdkXuEgnAOIGLxHmiABqb8IQkAMEOkI++Rbk+4b1sI54AN2vqoMYEGAAAkqkhPj4hke6cXJkRsEc6BCCXKqHFbP8AAAImE0TYgcTi5MCFii3AORCgRRo0T5QCDk/hDGEgcdvx/dvqAJ4tgQeJ3ExBvCOcAEuIAg9MS4Q9h/oizh9OhDJGL9P+zKQc8WakcifC7CWhLCOdAAkiEMJAI++DkarB2iPfz2kx4DU0JZYhcJP+fOeAJk3CQBogfhHMgziVCGIh0H0wIZZEyZXQjkvPanH4fTHgNIw1lTr+GAADAOYRzwADNjRrX1tZqZ4W0cXdZi6OViTBCE+k+2BHKTAhF8T66YcLlYtr6awgz/i9HKhH2AYmBWgRii3AOyNkp1S2PGier4IvVTW4
      "text/plain": [
       "<Figure size 1200x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = load_data(execution, \"flowRate\")\n",
    "simple = df[df.label == \"simple\"].groupby(df.direction)\n",
    "alpha = simple.get_group(\"[0.0, 0.0, 1.0]\")[\"alpha\"].to_numpy()\n",
    "\n",
    "k1 = (permeability(\"simple\", \"[0.0, 0.0, 1.0]\"))\n",
    "k2 = (permeability(\"simple\", \"[1.0, 0.0, 0.0]\"))\n",
    "k3 = (permeability(\"simple\", \"[1.0, 1.0, 1.0]\"))\n",
    "\n",
    "#poly = np.polynomial.Polynomial.fit(alpha, anisotropy_21, 1)\n",
    "\n",
    "fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize = (12, 6))\n",
    "ax.plot(alpha, k2 / k1, \"s\", label = \"$k_2$ / $k_1$\")\n",
    "ax.plot(alpha, k3 / k1, \"^\", label = \"$k_3$ / $k_1$\")\n",
    "#ax.plot(alpha, poly(alpha), \"-\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "#plt.ylabel(\"Анизотропия проницаемости\")\n",
    "#plt.title(\"Простая кубическая\")\n",
    "plt.show()\n",
    "\n",
    "if savefig:\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(\"anisotropy-simple.tiff\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "anisotropy[\"simple\"] = [(k2 / k1).mean(), (k3 / k1).mean()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "81.01844712765529"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k1.mean() / (0.9869 * 10**-12) * 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07995710547028301"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k1.mean() * (10**6)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Body-centered structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0cAAAGXCAYAAAB4GyuFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABEMUlEQVR4nO3de5xdZXno8d9DSAiQQIBAuBNFueihJBJQENvhGKxFIiKXigds0EM42ojctIJIbRWwlSiFepToaVPvhEvSBpRiwAEkaAAJgkBAIMpFKQkkZEgmN97zx9pDJ8OemTUze+21957f9/OZz2bWetdazzwZ9qxnv5cVKSUkSZIkabjbouwAJEmSJKkRWBxJkiRJEhZHkiRJkgRYHEmSJEkSYHEkSZIkSYDFkSRJkiQBsGXZARRl/PjxaeLEiWWH0RJeeeUVtt1227LDGHbMeznMe3nMfTnMezmaNe8bX90IwJZbNO8tZLPmvtk1Ut7vu+++5Smlnavta97f7H5MnDiRe++9t+wwWkJ7ezttbW1lhzHsmPdymPfymPtymPdymPfymPtyNFLeI+J3ve1zWJ0kSZJymbNkDnOWzCk7DKkwFkeSJEnKxeJIrc7iSJIkSZKwOJIkSZIkwOJIkiRJkgCLI0mSJEkCLI4kSZIkCWjh5xzl1dnZyQsvvEBnZycbN24sO5yGtP322/PII4+UHUZuW265JaNHj2bnnXdm9OjRZYcjSZKkJjGsi6NVq1bx/PPPs/POO7Prrruy5ZZbEhFlh9VwVq9ezdixY8sOI5eUEhs3bqSjo4Pf//73TJgwge23377ssCRJktQEhvWwuuXLl7Pnnnuyww47MHLkSAujFhARjBw5kh122IE999yTFStWlB2SJDWujheYdP+FsOqZsiORpIYwrHuO1q9fz9Zbb112GCrI1ltvzbp168oOQ5Ia1+LZbL/qYbjlIjhpTtnRqAlcd/J1ZYcgFWpY9xwB9ha1MP9tJakPGzrhl98gSLD0ZvjDA2VHpCYwfpvxjN9mfNlhSIUZ9sWRJEnD0oPXQno1+++NnbDgbEip1JDU+OYsmcOcJXPKDkMqjMWRJEnDTUpw+z/A+le6NsALj8Ljt5QalhqfxZFa3bCecyRJ0rD0xG2w9sXNt21YAzedC/s+ACO8PVB17dPbyw5BKpQ9R5IkDTd3/GO3XqNu1rwEv/q3+scjSQ3C4mgYmjJlCgcddFDZYbym0eKRpJb2X4/Ac0uq79vwCiz8AnS+XM+I1EQuX3Q5ly+6vOwwpMLkKo4i4oKIuDYinoyIFBHLBnOxiPhIRNwfEWsj4vmI+HZE7Fyl3ZzKdap9nTiYayuzceNGHnroISZNmlTodWbOnMkee+xB6mdyb73ikSRV3PlV2LSh9/2bNsAd3vyquhsfu5EbH7ux7DCkwuQdVHwp8CLwK2DcYC4UEecAXwVuBz4F7AmcCxweEYellKr073NalW2LB3N9ZR5++GHWrVtXaDGSUmL+/Pkcd9xx/S6nXY94JEndPP8QpE2979+4Fp5qr1s4ktRI8hZH+6aUngSIiIeAMQO5SESMB74E3AO8O6XsXTki7gH+g6xYurTncSml7w3kOurfkiVLAJg8eXJh17jnnnt49tln+cAHPtAQ8UiSuvnE3Zt9297eTltbWzmxSFKDyVUcdRVGQ/ABYBvgqq7CqHLeBRHxJHAqVYqjyLodxgIdKXU9jKGxTfnST1nesb7PNuPHjOLei46uU0Sbu//++wE266lZuXIlp59+Oj/5yU+48sormTFjxpCuMW/ePMaNG8dRRx3VEPFIkiRJedRrrc5DK693V9n3C+CUiBiTUurosW8VWXG0PiLuAC5KKf2ywDiHrL/CKG+boixZsoS99tqLHXfcEYD77ruPk046iZQSd911F4cccsiQrzFv3jyOOeYYRo4c2RDxSJIkSXnUqzjavfL6bJV9zwJRafNYZdsfga8B9wGvAAcDZwN3RsQxKaWF1S4SETOAGQATJkygvb29z6C23357Vq9ePZCfo2bKuu6SJUs4/PDDWb16Nd/+9re54IIL+LM/+zNmz57NjjvuWDWuTZs25Y536dKlLF26lM997nO5jhlMPAPR2dnZ7+9Bo+ro6Gja2JuZeS+PuS+HeS9Hs+Z95cqVAE0Ze5dmzX2za5a816s42qbyuq7Kvs4ebUgpfbZHm/kR8QNgCfAN4M3VLpJSmg3MBpgyZUrqbwz1I488wtixY/sJvRhlXHfZsmWsXLmSAw88kDPPPJNrrrmGL3zhC1x00UV9LpywevXq3PH+9Kc/ZfTo0Rx//PGMGdP31LTBxjMQo0ePbtr5TM4DKId5L4+5L4d5L0ez5n3csnEATRl7l2bNfbNrlrzX6zlHayqvW1XZN7pHm6pSSo8Dc4E3RcR+NYxt2Oia33PVVVdxww03cPPNN/P5z39+s0Jk3bp1nHHGGbzxjW9k7Nix7Lfffnzzm9/MfY158+YxderUfgujvPEAzJ07lyOPPJIxY8YwceLE3LFIkiRJA1Gv4ui5yuseVfbtAaRubfqyrPI6vgYxDTtdxcgZZ5zBunXrePDBB1/XZuPGjey6667ccsstrFq1irlz5/KVr3yFuXPn9nv+Z555hnvvvTfXKnV54wHYYYcdmDlzJpdcckmu80qSJEmDUa9hdfeQzQU6HPhtj33vAJZWWYyhmq7hdM/XMLZhY8mSJYwfP56vf/3rrF27lk9/+tNMnDiRD37wg6+12XbbbfniF7/42veTJk3imGOO4ec//zknn3xyn+efP38+EcH73//+msUDcPTRR792fkmSVJ726e1lhyAVquY9RxGxd0QcEBHdlyr7d2AtMDMiRnRrOw14I/D9btu2jYjR9BARk4GTgEdSSk/UOu7h4P77739t/s3VV19NW1sbp556KosX9/5c3Q0bNrBo0SL+5E/+pN/zz5s3jyOPPJKdd965sHgkSaqLjhfgX/8CVj1TdiSS6ihXcRQRp0XERRFxEbAzsH3X9xFxWo/m3wEeodsQupTSC8DngcOAhRExIyL+Dvgh8ChwRbfj3ww8FRHfiIhzI+LMiPi/ZMuAb6KyGp0GZsWKFTzzzDOvFSMjR47k+uuvZ+LEiUybNo2nnnqq6nEzZ85kzJgxfOQjH+nz/C+++CJ33HFH7iF1g41HkqS6WDwbfv8LuOWisiNpKJcvupzLF11edhhSYfL2HH0M+GLlaxdgXLfvP5bnBCmlWcDpwI7AlcDHyRZY+LMeQ+r+CCwEjgK+APwzcCxwDfC2lNLPc8ZcivFjRtWkTa11ze/pvnLbuHHjuOmmmwA45phjeOmllzY75txzz+Xuu+/mhhtuYNSovmO+8cYb2bhxI8cff3xh8UiSVBcbOuGX34D0Kiy9Gf7wQNkRNYy7n7mbu5+p9thKqTXkmnOUUmrLe8K+2qaU5gBz+jn+j0DP3qimce9FR5cdQlVTp04lpfS67W94wxt4/vnXT+E6++yzufXWW7ntttsYPfp1oxxfZ968eRx88MG5V5MbaDySJNXNg9dmhRHAxk5YcDaccRvU6DETzez6k68vOwSpUPVarU5N5KyzzmLhwoXcdtttuecPHX744Vx66aWFxLNp0yY6OzvZsGEDKSU6OztZt67aI7MkSRqilOD2f4D1r3RtgBcehcdvKTUsSfVRr9Xq1CR+97vfcdVVV7HVVlvxhje84bXt73rXu/jJT37S63Gf+cxnCovpu9/9Lqeffvpr32+99dbss88+LFu2rLBrSpKGqSdug7Uvbr5twxq46VzY9wEYMbxvnS5YeAEAl029rORIpGLYc6TN7LPPPq/1znR0dNDR0cEf/vCHPgujok2fPp2U0mZfFkaSGoarmrWWO/6xW69RN2tegl/9W/3jaTDOOVKrsziSJGkoXNWsdfzXI/Dckur7NrwCC78AnS/XMyJJdWZxJEnSYBW9qpm9UvV151dh04be92/aAHe4jLXUyiyOJEkarGqrmlVZiXPQ7JWqr+cfgrSp9/0b18JT7XULR1L9De9ZhZIkDVZfq5rt9+dDP3+1XqndDh76edW7TziXRhru7DmSJGk
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = load_data(execution, \"flowRate\")\n",
    "bodyCentered = df[df.label == \"bodyCentered\"].groupby(df.direction)\n",
    "alpha = bodyCentered.get_group(\"[0.0, 0.0, 1.0]\")[\"alpha\"].to_numpy()\n",
    "\n",
    "k1 = filter_group(permeability(\"bodyCentered\", \"[0.0, 0.0, 1.0]\"))\n",
    "k2 = filter_group(permeability(\"bodyCentered\", \"[1.0, 0.0, 0.0]\"))\n",
    "k3 = filter_group(permeability(\"bodyCentered\", \"[1.0, 1.0, 1.0]\"))\n",
    "\n",
    "#poly = np.polynomial.Polynomial.fit(alpha, anisotropy_21, 10)\n",
    "\n",
    "fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize = (12, 6))\n",
    "ax.plot(alpha, k2 / k1, \"s\", label = r\"$k_2$ / $k_1$\")\n",
    "ax.plot(alpha, k3 / k1, \"^\", label = r\"$k_3$ / $k_1$\")\n",
    "ax.axvline(0.13, linestyle = \"-.\", color = \"green\")\n",
    "#ax.plot(alpha, poly(alpha), \"-\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "#plt.ylabel(\"Анизотропия проницаемости\")\n",
    "#plt.title(\"Кубическая объемноцентрированная\")\n",
    "\n",
    "if savefig:\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(\"anisotropy-bodycentered.tiff\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "anisotropy[\"bodyCentered\"] = [(k2 / k1).mean(), (k3 / k1).mean()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Face-centered structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtoAAAGDCAYAAAAVh7eRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6i0lEQVR4nO3de5hcVZno/+9LyAXoAMFgDBEIAgJ6FBiCguJMZ0QHceSig0cZQGYc0GGCBh1vI3qYEdDxJ8qQ4ajR3xHFC5IjiQOMiFwaGAETIkEuAbkFBlAkhFwa0kmnWeePvTtWmuru6k7tXV1V38/z1FOptdfe661a6e631l577UgpIUmSJKm+tml0AJIkSVIrMtGWJEmSCmCiLUmSJBXARFuSJEkqgIm2JEmSVAATbUmSJKkA2zY6gKJMnTo1zZw5s9FhtJznn3+eHXbYodFhaBTsu+ZkvzUv+6552XfNqZH9tnTp0pUppV0Hlrdsoj1z5kzuuOOORofRcrq6uujs7Gx0GBoF+6452W/Ny75rXvZdc2pkv0XEY9XKnToiSZIkFaCmRDsiPhMRCyLikYhIEbFipA1FxIci4gcRcX9E9EXEoLekjIhz8naqPf5xpG1LkiRJZat16sj5wCrg18DOo2zrM8DLgDuBHYBX1rDPWcDKAWVLR9m+JEmSVJpaE+29U0qPAETEPUDHKNrqBB5PKb0YEVdRW6K9KKW0YhRtSZIkSQ1V09SR/iR7a6SUVqSUXhzpfhGxY0S07EWbkiRJak1j/WLI3wBrgJ6IuDUi3tHogCRJkqRaREqDXpNYfYd86khKaeaoG82mjrwzpRSDbJ8LHADcCjwH7AfMBaYDf5tSumSQ/U4HTgeYNm3aIZdddtloQ9Qguru76egYzcwhNZp915zst+Zl3zUv+645NbLfZs+evTSlNGtg+ZhMtAfZ52XAPcAkYPeUUvdQ9WfNmpVcR7v+XFu0edl3zcl+a172XfOy75pTg9fRrppoN83c55TSsxHxDeAc4E3AtY2N6I9mnfsLVnZvHLLO1I4J3HH225qqLUmSJI1e0yTauRX589RGBjHQcIlvrXXGWluSWp9f3iWpOM2WaO+bPz/d0CgkbcFkrXn55V2SilP3VUciYo+I2D8ixo9y/20jYqcq5bsDfw88S3aRpKQxwmRNkqSXqmlEOyJOBvbMX+4KTIiIs/PXj6WULq2o/j3gz4C9+ONUDyLiXcCB+ct98rL+Y6xOKf17/u8O4NGIWAQs54+rjvxdvu39KaX1Nb4/SdIY4tkPSe2k1qkjHyRLnit9IX++CbiU4b0H+MAgx3gM6E+01wM/Ad4IHEeWXK8ErgO+nFJaXGPMkqQxxrMfktpJTYl2Sqmz1gMOVjeldCpwag37byAbvZYkSZKa1li/M6QkSZLUlJpt1RFJKp3ziiVJo2GiXQdTOybU9Ee42dqSlClzXnHVpP6aq7d4aVIvSc3BRLsOyvyD5x9XqbWVfbGgX94lqTgm2pK2msla8/LLuyQVx0Rb0lYrO1lzzrQkqRmYaEtqOq7F3LzKPPvhFzJJjWaiLUkqTZlJrV/IJDWaibaG5AoIkiRJo2OirSGVPSLkqd768UuSJEmNZaKtMaWVT/WW/SWilT/LsrmqiiRpNEy0pZKY+DavMkf9Teqbl2eRJA1koi1JY8jAJKyrq4vOzs7GBKMR8cu0pIFMtNW2nA8uSQLPRqg4JtpqW44+NS+nV0gOFtSTfw9UFBNtSU3HxEG1aPUvZCaH0thnoi1Jakl+IaufskfPW3m0vpXfm15qm0YHoLGtltGeZh4RkiQNr+zR81YerW/l96aXckRbQyp7BYRWP9VbJj9LSZIay0RbY0ornyorO/F1mThJkhrLRFsqSSt/iZDkWSSNTc4JbywTbUmS6sCzSBqLnBPeWF4MqbblhZ6SJPDvgYrjiLbalqfJJDUzp6rUj2cjVBQTbUmSmpCDBdLY59QRSZI0pLKnVrTyVI5Wfm96KUe0JUnSkMoePW/l0fpWfm96KUe0JUmSpAKYaEuSJEkFMNGWJElqUc4JbyznaEuSJLUo54Q3Vk0j2hHxmYhYEBGPRESKiBUjbSgiPhQRP4iI+yOiLyLSMPX3i4hFEfFcRDwfEbdExJ+PtF1JkiSpEWod0T4fWAX8Gth5lG19BngZcCewA/DKwSpGxN7ArcAm4MvAGuA04OcR8Y6U0nWjjEGSJEkqRa2J9t4ppUcAIuIeoGMUbXUCj6eUXoyIqxgi0Qa+SJbQH5JSWpa3+z3gXuDiiNg/pTTkiLgkSZJa16xzf/HSu6Nec/UWL6d2TGjo9Jmapo70J9lbI6W0IqX04nD1ImIH4Bigqz/JzvfvBr4NvBo4dGvjkSRJUvN6SZI9yjpFGourjrwemAjcVmXb7fmzibYkSZLGtLG46shu+fOTVbb1l82otmNEnA6cDjBt2jS6urrqHly76+7u9nNtUvZdc7Lfmpd917zsu9bSyL4ci4n29vnzhirbegbU2UJKaT4wH2DWrFmps7Oz7sG1u66uLvxcm5N915zst+Zl3zUv+65JDJiPPZhG9uVYnDryQv48scq2SQPqSJIkSWPSWEy0n8qfq00P6S+rNq1EkiRJGjPGYqJ9N9m0kcOrbDssf76jvHAkSZKkkat7oh0Re0TE/hExfjT758v4XQl0RsSBFcftAP4OeBBYXJdgJUmSpILUdDFkRJwM7Jm/3BWYEBFn568fSyldWlH9e8CfAXsBKyqO8S6gP3HeJy/rP8bqlNK/VxzjM8BbgWsj4mvAWrI7Q84A3unNaiRJktrb1I4Jw66TPbVjQknRVFfrqiMfJEueK30hf74JuJThvQf4wCDHeAzYnGinlB6KiDcDXwI+DUwgu/37Ud5+XZIkSQPv+DgWV4upKdFOKXXWesDB6qaUTgVOHcFxlgPH1lpfkiRJGkvG4sWQkiRJUtMz0ZYkSZIKYKItSZIkFcBEW5IkSSqAibYkSZJUABNtSZIkqQAm2pIkSVIBTLQlSZKkAphoS5IkSQUw0ZYkSZIKYKItSZIkFcBEW5IkSSqAibYkSZJUABNtSZIkqQAm2pIkSVIBTLQlSZKkAphoS5IkSQXYttEBSJIkqTXMOvcXrOzeOGSdqR0TuOPst5UUUWM5oi1JkqS6GC7JrrVOqzDRliRJkgpgoi1JkiQVwERbkiRJKoCJtiRJklQAE21JkiSpACbakiRJUgFMtCVJkqQCmGhLkiSpLqZ2TKhLnVbhnSElSZJUF+1yx8daOaItSZIkFcBEW5IkSSqAibYkSZJUgJoS7Yj4TEQsiIhHIiJFxIrRNBYRp0TEnRGxPiKejohvR8SuVepdkrdT7fFXo2lbkiRJKlOtF0OeD6wCfg3sPJqGIuIs4KvATcBHgVcCHwMOj4g3pJSer7LbyVXKFo+mfUmSJKlMtSbae6eUHgGIiHuAjpE0EhFTgXOBJcBbU0p9efkS4D/IEu/zB+6XUvr+SNqRJEmSxoqapo70J9lb4Thge2Bef5KdH/dK4BHgpGo7RWbHiHAuuSRJkppKWQnsofnzbVW23Q7sHxHVRsnX5I/1EfGLiHhjUQFKkiRJ9RQppZHtkE8dSSnNHME+VwJ/CWyfUlo/YNuXgU8A+6WUfpuXfQmYACwFngcOBOYCOwBHp5SuG6Sd04HTAaZNm3bIZZddNqL3puF1d3fT0TGimUMaI+y75mS/NS/7rnnZd82pkf02e/bspSmlWQPLy7oz5Pb584Yq23oG1CGl9OkBdRZFxA+BZcDXgX2rNZJSmg/MB5g1a1bq7OwcfcSqqqurCz/X5mTfNSf7rXnZd83LvmtOY7Hfypo68kL+PLHKtkkD6lSVUnoQuBzYJyJeXcfYJEmSpLorK9F+Kn+eUWXbDCBV1BnKivx5ah1ikiRJkgpTVqK9JH8+vMq2w4AHUkrdNRynf8rI03WJSpIkSSpI3RPtiNgjIvaPiPEVxT8F1gNzImJcRd13Aa8CflBRtkNETGKAiDgYOAFYnlJ6uN5xS5IkSfVU08WQEXEysGf+cldgQkScnb9+LKV0aUX17wF/BuxFPtUjpfRMRHw
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = load_data(execution, \"flowRate\")\n",
    "faceCentered = df[df.label == \"faceCentered\"].groupby(df.direction)\n",
    "alpha = faceCentered.get_group(\"[0.0, 0.0, 1.0]\")[\"alpha\"].to_numpy()\n",
    "\n",
    "k1 = filter_group(permeability(\"faceCentered\", \"[0.0, 0.0, 1.0]\"))\n",
    "k2 = filter_group(permeability(\"faceCentered\", \"[1.0, 0.0, 0.0]\"))\n",
    "k3 = filter_group(permeability(\"faceCentered\", \"[1.0, 1.0, 1.0]\"))\n",
    "\n",
    "#poly = np.polynomial.Polynomial.fit(alpha, anisotropy_21, 10)\n",
    "\n",
    "fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize = (12, 6))\n",
    "ax.plot(alpha, filter_group(k2 / k1, quantile = 0.97), \"s\", label = r\"$k_2$ / $k_1$\")\n",
    "ax.plot(alpha, filter_group(k3 / k1, quantile = 0.8), \"^\", label = r\"$k_3$ / $k_1$\")\n",
    "#ax.plot(alpha, poly(alpha), \"-\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "#plt.ylabel(\"Анизотропия проницаемости\")\n",
    "#plt.title(\"Кубическая гранецентрированная\")\n",
    "plt.show()\n",
    "\n",
    "if savefig:\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(\"anisotropy-facecentered.tiff\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "anisotropy[\"faceCentered\"] = [(k2 / k1).mean(), (k3 / k1).mean()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Anisotropy mean values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>simple</th>\n",
       "      <th>bodyCentered</th>\n",
       "      <th>faceCentered</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A21</th>\n",
       "      <td>2.2591</td>\n",
       "      <td>0.7571</td>\n",
       "      <td>1.1312</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A31</th>\n",
       "      <td>1.5660</td>\n",
       "      <td>1.0036</td>\n",
       "      <td>0.9148</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     simple  bodyCentered  faceCentered\n",
       "A21  2.2591        0.7571        1.1312\n",
       "A31  1.5660        1.0036        0.9148"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "anisotropy.round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analytical porosity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def porosity_a(a, label):\n",
    "    if label == \"simple\":\n",
    "        val = (1 - np.pi * (2 - 3 * a ** 2 * (3 - a)) / (12 * (1 - a) ** 3))\n",
    "\n",
    "    elif label == \"bodyCentered\":\n",
    "        a13 = a[a <= 0.134]\n",
    "        a13g = a[a > 0.134]\n",
    "        a = a13\n",
    "        val = (1 - np.pi * np.sqrt(3) * (1 - 2 * a ** 2 * (3 - a)) / (8 * (1 - a) ** 3)\n",
    "        )\n",
    "        a = a13g\n",
    "        val2 = (1 - np.pi * np.sqrt(3) * (1 - 2 * a ** 2 * (3 - a)) / (8 * (1 - a) ** 3) +\n",
    "            np.pi * (np.sqrt(3) - 2 * (1 - a)) ** 2 * \n",
    "            (np.sqrt(3) - (1 - a)) / (8 * (1 - a)**3)\n",
    "        )\n",
    "        val = np.append(val, val2)\n",
    "\n",
    "    elif label == \"faceCentered\":\n",
    "        val = (1 - np.pi * (1 - 3 * a ** 2 * (3 - a)) / (3 * np.sqrt(2) * (1 - a) ** 3))\n",
    "\n",
    "    return val\n",
    "\n",
    "def pad(x, y):\n",
    "    return np.pad(y, (0, x.size - y.size), 'constant', constant_values = np.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = np.linspace(0, 0.3, 50)\n",
    "porosity_simple = porosity_a(alpha[alpha <= 0.2929], \"simple\")\n",
    "porosity_bodyCentered = porosity_a(alpha[alpha <= 0.1835], \"bodyCentered\")\n",
    "porosity_faceCentered = porosity_a(alpha[alpha <= 0.134], \"faceCentered\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAGKCAYAAAC/yWNIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABedElEQVR4nO3dd3gUVdvH8e9JJ4QWepUuSBcEGwKKCoKNoqKgVLEgIpZH1NdHfeyKWBBEUFFQQBQVBbGgAQsiKF1ApEnvLYTUPe8fs5CEbEICm51N8vtc116wZ9qdk8nmzsyZ+xhrLSIiIiL+EuJ2ACIiIlK4KLkQERERv1JyISIiIn6l5EJERET8SsmFiIiI+JWSCxEREfErJRciIiLiV64mF8aYEGPMfcaYNcaYRGPMFmPMSGNM8Vxub7N5xed37CIiIuJbmMvHHwUMBT4DRgINve9bGGM6Wms9udjHT8DbJ7Wl+DVKERERyTXXkgtjTCPgHmCGtbZ7hvaNwOvATcBHudjVBmvt5PyJUkRERPLKzdsivQADvHpS+3ggAeid2x0ZYyKMMTH+C01EREROl5vJxXmAB/g9Y6O1NhFY6l2eGz1wkpEjxpjdxpg3jDGl/BmoiIiI5J6bYy6qAHuttUk+lm0DLjTGRFhrk3PYx+/AdOAfoCRwFTAEaGeMudBaq4GdIiIiAeZmchEN+EosABIzrJNtcmGtbXNS0wfGmOXAM8C93n99MsbcDtwOUKxYsZbVq1fPZdin5vF4CAnRU76+qG98U79kT33jm/ole+ob3zL2y99//73XWls+v45l3Jpy3RizAqhgra3oY9nHQE8g8hRXLnztNxyIB/6w1l6Ym21atWplFy9enJfD5CguLo727dv7bX+FifrGN/VL9tQ3vqlfsqe+8S1jvxhj/rDWtsqvY7mZ2m0HyhljIn0sq4pzyyRPiQWAtTbl+L7PMD4RERE5DW4mF4u8x2+dsdEYEwU0B07rUoJ3+2rArjOMT0RERE6Dm8nFNMACw05qH4Qz1uLD4w3GmDrGmAYZVzLGlM1mv//DGUvypd8iFRERkVxzbUCntXaFMeZNYIgxZgYwm/QKnfPIXEBrLnAWTl2M4x4zxpwP/Aj8C8TgPC3SAVgIvJHvX4SIiIhk4Xb572HAJpynNroAe3GSgsdzUfo7DjgHuA0oC6QB64BHgVe89TJEREQkwFxNLqy1aThziow8xXo1fbR9AXyRP5GJiIjI6dKDwCIiIuJXbt8WERGRAiQxMZE9e/aQmJhIamqq2+H4VKpUKVavXu12GAEVFhZGVFQU5cuXJyoqyu1wlFyIiEjuHDp0iF27dlG+fHkqVapEWFgYxphTbxhgR44coUSJEm6HETDWWlJTU4mPj+fff/+lYsWKlCrl7hRbSi5ERCRX9u7dS7Vq1YiOjnY7FMnAGEN4eDhlypQhMjKSnTt3up5caMyFiIjkSnJyMsWKFXM7DMlBsWLFSErKbtquwFFykQ/+2JXKN6t2uh2GiIjfBeNtEEkXLN8fJRd+tvNQIu+sSGLwpD94cPoy4pOCc8CTiIhIflFy4WcvzllDgjefWLBhH27NOisiIuIWJRd+9t+rG3F+5VCMgVduaE6JqHC3QxIREQkoJRd+Vio6nDuaRTF7aFta14rNsvxwYooLUYmISF7FxcVhjOHll1/Osiw1NZXu3btjjOHuu+8G4IknnsAYc8qXr/0VNnoUNZ80rFwyS9vfu47QfeyvDL20HgMurkVISHAMvBERkdxLTU2lV69ezJgxgzvuuIPRo0dnWv7UU09Rq1atLNutXr2aZ599NlBhukrJRYAkp3oYNnUpRxJTeWb2arYcSOCpaxu7HZaIiORBWloavXv35pNPPmHQoEGMGTMmyxManTt3plWrVlm2jYuLKzLJhW6LBMiBhGRCvVcqIsJC6H3+WS5HJCIieZGWlkafPn2YNm0a/fv3Z9y4cUHz6Gew0ZWLAKlYMopP77yQ1+euo3yJSOpXLDqlaUVECjqPx0Pfvn2ZMmUKffv2Zfz48UoscqDkIoAiwkJ44MqzfS77Y/N+ElM8XFS3XICjEhGRnHg8Hvr168fkyZMB6NOnDyEhuvCfEyUXQeBwYgpDpyxl28Fj9L+oFg91Opuo8FC3wxIRyZNR3/3Na3PXAXDvZfW47/L6mZY//dVfTPh5IwCPXtWQQZfUzrR8xIzlTPl9CwDPXt+Em9vUyLR86JQlzFy2HYDXbmrOtc2rZlo+YOIi5q7ZDcCEW1vR8ZyKfvm6Ro4cye7du+nRowfz5s1j4MCBLF++nJiYGL/svzBS6hUERn6zlm0HjwHw6Z9bOZigx1VFRILF7t27ufHGG5k6dSpjx45l48aNDB8+3O2wgpqSiyAw5NJ6XNqgAgDPXN+YSqWiXI5IRESO69q1K5MnTyY0NJTu3bvTq1cvxo8fz5w5c9wOLWjptkgQKF8ikndua8Uv/+zj4npZx1ykpnkIC1UeKCLB7b7L62e5FZLRY13P4bGu52S7/LluTXmuW9Nsl7/eqwWv92qR7fJ3+p4HwJEjRyhRwn+D5tu1a0dYWPqvy9GjRxMXF8fAgQNZuXIlpUuX9tuxCgv9xgoSxhificWuw4lc9so8vli6zYWoRETkZLGxsYwfP55t27YxZMgQt8MJSkougpjHY3lg+jI270vg3qlLeWbWX26HJCIiQJcuXejfvz8ffvghn332mdvhBB0lF0FsT3wSm/clAGAMdPCOyxAREfeNGjWKs846i8GDB7N79263wwkqGnMRxCqWjGL2vW15cuYqYotHcGEd1cAQEQmU9u3bY63NdnnJkiXZtGnTifdPPPEETzzxxGnvrzBRchHkYiLDeKlnMzyerCfkxr1HSU3zUE/VPkVEJIjotkgBcfIMqsmpHoZOWULXN37mgwWbikw2LCIiwU/JRQH1+tx1rNh2iKRUD09/tfrE2AwRERG3KbkooK5uVoUGlZzbIQ91Opua5Yq7HJGIiIhDYy4KqLMrleCLIRfx8eKt3NK6xqk3EBERCRBduSjAIsNC6XP+WVnGY8QnpdLvvd9ZsfWQS5GJiEhRpuSiEHpy5ip+XLuH68f8wnu/bHQ7HBERKWKUXBQy2w8eY/aKHQCkeiyxxSNcjkhERIoaJReFTJXSxZh9b1ta1CjNdc2rcG3zqm6HJCIiRYwGdBZCZ5UtzseDLyAlzZNl2f6jyUSEhRATqW+9iIjkD125KKTCQ0OIjsicQFhruf/jpVz12k/8+e8BlyITEZHCTslFETLpt838uHYP/+5P4Ia3FrBlvwpviYiI/ym5KEJii0dQIsq5mtH3wppUj412OSIRESmMdOO9COnatArNq5fmzR//4cFOZ7sdjoiIFFK6clHEVCsTzXPdmhIZFpqpPSXNw6vf/83hxBSXIhMRCS5xcXEYY3j55ZezLEtNTaV79+4YY7j77rszLdu4cSN33XUX9erVo1ixYpQuXZoLLriAV199laSkpCz7mjhxIsYYPvnkE59xbNq0CWMMQ4YM8c8XFgBKLgSAN+au49Xv13HVaz/xx+b9bocjIhK0UlNT6dWrFzNmzOCOO+5g9OjRJ5Z98cUXNGrUiEmTJnHllVfy+uuv8+STT1KyZEnuu+8+LrjgAnbt2uVi9IGh2yLChj3xjP7xHwC2HjjGwo37aXlWrMtRiYgEn7S0NHr37s0nn3zCoEGDGDNmDMY4UzCsXLmSm266iXLlyhEXF0edOnVObHfvvfcyfvx4br/9dm644YYTV0UKK125EGqXj+H1Xi0oGRXG+bVjGXxJnVNvJCJSxKSlpdGnTx+mTZtG//79GTduXKYE4b///S+JiYmMGzcuU2Jx3KBBg+jZsyfz589n1qxZgQw94JRcCOAM9vx62CW8emMLQk+aCC3VRzEuEZGixOPx0LdvX6ZMmULfvn0ZP358psQiMTGRWbNmUa1aNa666qps9zNo0CAAPv3003yP2U26LSInVC1dzGf74zNXcSw5jaeubUSJqPAARyUiBcaPz8G853O37rm3wTWvZ26bORT+fD9327d7GDqMyNz20Y3w9xxK5LTOafB4PPT
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize = (8, 6))\n",
    "\n",
    "ax.plot(alpha, pad(alpha, porosity_simple), \":\", label = \"КП\")\n",
    "ax.plot(alpha, pad(alpha, porosity_bodyCentered), \"--\", label = \"КОЦ\")\n",
    "ax.plot(alpha, pad(alpha, porosity_faceCentered), \"-.\", label = \"КГЦ\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$m$\")\n",
    "plt.show()\n",
    "\n",
    "if savefig:\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(\"porosity-analytical.tiff\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = load_data(execution, \"porosity\")\n",
    "alpha = df[\"alpha\"].unique()\n",
    "simple = df[df.label == \"simple\"][df.direction == \"[0.0, 0.0, 1.0]\"][\"porosity\"].to_numpy()\n",
    "bodyCentered = df[df.label == \"bodyCentered\"][df.direction == \"[0.0, 0.0, 1.0]\"][\"porosity\"].to_numpy()\n",
    "faceCentered = df[df.label == \"faceCentered\"][df.direction == \"[0.0, 0.0, 1.0]\"][\"porosity\"].to_numpy()\n",
    "\n",
    "part_simple = pad(alpha, simple)[alpha <= 0.29] / porosity_a(alpha[alpha <= 0.29], \"simple\")\n",
    "part_bodyCentered = pad(alpha, bodyCentered)[alpha <= 0.18] / porosity_a(alpha[alpha <= 0.18], \"bodyCentered\")\n",
    "part_faceCentered = pad(alpha, faceCentered)[alpha <= 0.13] / porosity_a(alpha[alpha <= 0.13], \"faceCentered\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABoAAAAQCAYAAAAI0W+oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABJ0AAASdAHeZh94AAAB00lEQVR4nL3UPWzNYRTH8U8bCUm9JRoaA+KlIiIiMTQIg8TgZTBZLAYGE2kXluNEBEkHkQ5MEkwWJiQkommMHQwGgyqLl6ANDUtTw30u1/W/tyTit5w85/md55vz/5/n6ZiZmfE/NKdxkZkXsRW96MZXjOMOhiLiQ5P/JVa2OPttRPRUgnASo3iAd+hCH87gWGb2RcTrpppJXKoAfWnZERZGxLfmisw8h9M4heNN2xMRcaYC9Is6GxdVkKJbJa6b7cBWau6olQ6U+LRib25mHsYKTBXPcERMzwrKzAHMxyK14dhRDrhQYe/BjabcWGYeiYjH9USnag0gcKJA7mNPRLxv8l3D7gLrwiZcxSrcy8zNdWNHu3uUmcuwrXSyAPsjYrRlwc+6QfTjTkQcbNcRiIi3EXEbe7AE12eDFF0pcWc90RbUABzHM2zMzO4/KKl/4q6/AhUtL3G6raumvhJf1BM/pi4ze9WejcnGiszsxFksxZOI+FTyG/AqIqaa/KswVJY3fwNhL85n5gjG8AHLsAur8QZHG/yH0J+Zw2rv4WeswT7Mw10MVoEeYq3aOG/BYrUL+FztnlyOiI8N/kdYX7zb1f7HBEaK/0ZE/BjptuP9L/UdP8eNfEo28QwAAAAASUVORK5CYII=",
      "text/latex": [
       "$\\displaystyle 35$"
      ],
      "text/plain": [
       "35"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "porosity_a(alpha[alpha <= 0.18], \"bodyCentered\").size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAGDCAYAAAAs+rl+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3M0lEQVR4nO3de5xcVZnv/8+TJuRKiDEXkjgQRSQORogEBDmSRlHkkJ8CDUKYZGhyEogSkKCZEUEnHh3hDAnJkQjBMAyDQRISQxgPEA4XGxhBIdxRQG6tTudAujWBNN2dS/f6/bF3JdXVddldt71r7+/79apX0WvfVq0U/fRae61nm3MOERERia8BYVdAREREKkvBXkREJOYU7EVERGJOwV5ERCTmFOxFRERiTsFeREQk5vYLuwKVMnr0aDdp0qReZe+//z7Dhg0Lp0I1TO1WHLVbcdRuxVG7FSdu7fb000+3OefGZJbHNthPmjSJzZs39ypramqivr4+nArVMLVbcdRuxVG7FUftVpy4tZuZ/TFbuYbxRUREYk7BXkREJOYU7EVERGJOwV5ERCTmFOxFRERiTsFeREQk5hTsRUREYi626+yD6urqorW1la6uLvbs2RN2dSLpwAMP5OWXX67qNQcOHMjYsWMZMWJEVa8rIhJHiQ727777Lu+88w5jxozhoIMOYr/99sPMwq5W5OzYsYMDDjigatdzztHZ2UlLSwuAAr6ISIkSHezb2tr40Ic+xNChQ8OuiqQxM4YOHcrEiRPZsmWLgr0Udu1h8P7WvuXDxsKi16pfH5GISfQ9+127djFkyJCwqyE5DBkyhN27d4ddDakF2QJ9vnKRhEl0sAc0bB9h+rcRESmPxAd7ERGRuFOwFxERiTkFexERkZhTsE+QpqYmzIwlS5b02bZnzx4aGhowMy6++GIAFi9ejJkxYsQIzCznK9v5RKpq2Nj+lYskTKKX3olnz549zJw5kw0bNjB//nxWrFjRa/uVV17J5MmT+xz38ssv86Mf/aha1RTJTcvrRPJSsE+47u5uZs2axfr165k3bx433HBDn1nwX/jCF5g+fXqfY5uamhTsRURqgIbxE6y7u5vZs2ezdu1a5syZw0033aTlbiIiMaSefYVsfLaFa+9/lS3bO5kwcgiLTjmc06dODLtae/X09NDY2Mgdd9xBY2Mjq1atUqAXEYkpBfsK2PhsC1dseJHO3d0AtGzv5IoNLwJEIuD39PRwwQUXsHr1agBmz57NgAEa5BERiSv9hq+Aa+9/dW+gT+nc3c21978aUo16W7p0KbfddhtnnXUWY8aMYe7cubS3t4ddLRERqRAF+wrYsr2zX+XVtnXrVs455xzWrFnDjTfeyFtvvcXll18edrVERKRCFOwrYMLI7A/XyVVebTNmzGD16tXU1dXR0NDAzJkzWbVqFZs2bQq7aiIiUgEK9hWw6JTDGTKwrlfZkIF1LDrl8JBq1Nv06dPZb7990zVWrFjB+PHjmTt3Ltu3bw+vYiIiUhEK9hVw+tSJXH3mFCaOHIIBE0cO4eozp0Ricl42o0aNYtWqVbS0tLBgwYKwqyMiImWm2fgVcvrUiZEN7tmcdtppzJkzh1tuuYWGhgbOOOOMsKskIiJlop697LVs2TIOOeQQLrroIrZu3Rp2dUREpEzUs0+Q+vp6nHM5t48YMYLm5ua9Py9evJjFixezY8eOos4nIiLRoJ69iIhIzCnYi4iIxJyCvYiISMwp2IuIiMScgr2IiEjMKdiLiIjEnIK9iIhIzCnYi4iIxJyCvYiISMwp2IuIiMRc4GBvZleY2Toze9PMnJk1F3NBM/t7M3vWzDrN7B0zu9nMxuTY99Nm9qCZ7TCz98xsk5kdVcx1RUREkqo/PfsfAZ8D3gC2FXMxM1sI/DvwLvAN4CbgXKDJzIZl7Hsc8AjwYeB7wD8BhwGPmdmUYq4vIiKSRP15EM6hzrk3AczsJWB4fy5kZqOBHwJPAZ93znX75U8B/4EX/H+UdsiPgV3Aic65Fn/fO4GXgaXAF/tzfRERkaQK3LNPBfoSnA4MBa5PBXr/vL8E3gRmpcrM7KPAMcC6VKD3920B1gEnm9lBJdYncZqamjAzlixZ0mfbnj17aGhowMy4+OKLe21rbm7m61//OocddhhDhgxh5MiRHH/88SxfvpydO3f2Odett96KmRV8LViwoGKfVURE9qnmI26P8d+fyLLtN8BMMxvunGsPsO8c4GjgnrLXMoH27NnDzJkz2bBhA/Pnz2fFihV7t919993MnDmTuro6zj//fI488kg6Ojq49957WbhwIbfddhv33Xcf48aN63PeSy+9lGOOOaZPeVtbGwsXLqzoZxIRkX2qGewn+O8tWba1AObv84cA+wJMLGvtEqq7u5tZs2axfv165s2bxw033ICZAfDSSy9x7rnn8sEPfpBHHnmEQw89dO9x3/jGN1i1ahUXXnghX/3qV/eOGqT77Gc/y1lnndXnms3NzQr2IiJVVM1gP9R/7zvuC10Z+/Rn373M7ELgQoBx48bR1NTUa3t7e3uvsgMPPJAdO3YUrnlMdHR0ALBz50527NhBd3c38+bNY/369cyePZslS5bQ3t6+d/8rr7ySrq4urrvuOsaOHdunrc4991zuu+8+7rrrLtatW8epp54KQFeX90/U2dmZtX1T19i1a1fB9u/q6urz71grMr9vEozarThqt+Ikpd2qGew7/PdBQGfGtsEZ+6Tvmylz372ccz8Ffgowbdo0V19f32t7U1MT6WUvv/wyBxxwQKDK98u1h8H7W/uWDxsLi14r//UCGjrU+/to0KBBDBs2jPPPP5/169fT2NjIv/7rvzJgwL4pHF1dXdx///186EMf4tRTT83ZTl/72te46667uO+++/jqV78KwODB3j/RkCFDsh43fLg3t3P//fcv2P6DBw9m6tSp/f+wEZD5fZNg1G7FUbsVJyntVs2kOlv892zD7xMBl7ZPoX0h+xB/NGQL9PnKq6ynp4cLLriA1atXAzB79uxegR7gtddeY+fOnQUD7dFHHw3Aiy++WJnKiohIyaoZ7J/y34/Psu044FV/cl6QfR3wdHmrlxxLly7ltttu46yzzmLMmDHMnTu31/A9wHvvvQd4tzryGTFiBADvvvtuZSorIiIlq0iwN7ODzWyymQ1MK74bb/h+gZnVpe37/wEfAW5PlTnnXgc2A2eb2YS0fScAZwMPO+ferkTdk2Dr1q2cc845rFmzhhtvvJG33nqLyy+/vNc+QYN40D8KREQkPIHv2ZvZbOAQ/8cxwP5mdpX/8x+dcz9L2/02YDpe9rtmAOdcq5l9F1gCPGhmd+ANyX8TeAVYnnHJbwC/wsuYd71fdgneHyjfDFpv6WvGjBmsXr2auro6GhoamDlzJqtWreLMM8/kS1/6EgCHHXYYgwYN4tlnn817rmeeeQaAKVOU1FBEJKr607P/H8AP/NdYYGTaz/8jyAmcc0uBC4BReBnyvgbcCUxPG8JP7fs4UI/3x8IP/eu8jpdR7/l+1FsyTJ8+nf322/d33ooVKxg/fjxz585l+/btgDcx7tRTT+W//uu/eOCBB3Ke6+abbwbgzDPPrGidRUSkeP3JoFfvnLMcr/oc+zZnOc+tzrkjnXODnXNjnXNznHNZZ645555wzn3eOTfcOXeAc+4U59wz/f2QVTdsbP/KQzZq1ChWrVpFS0tLr6x23//+9xk0aBCXXXYZzc3NfY675ZZbWLt2LSeeeCIzZsyoYo1FRKQ/qrn0LjlCXF5XrNNOO405c+Zwyy230NDQwBlnnMEnP/lJbr/9dmbNmsWUKVNobGzcm0HvvvvuY9OmTRx11FHceeedfRLqiIhIdOh59rLXsmXLOOSQQ7jooovYutUbbGloaOCJJ57g7/7u77j33ntZsGABV111Fdu2beO6667jiSeeyJoqV0REokM9+wSpr6/HOZdz+4gRI7IO1x966KGsXLky8HUaGxtpbGzMuX3SpEl56yEiIuWlnr2IiEjMKdiLxExrRyuNmxpp62wLuyoiksXGZ1s44ZqH+fC37+GEax5m47OVTwirYC8
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams.update({\n",
    "    \"lines.markersize\": 6.0\n",
    "})\n",
    "fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize = (8, 6))\n",
    "\n",
    "ax.plot(alpha, pad(alpha, part_simple), \"o\", label = \"КП\")\n",
    "ax.plot(alpha, pad(alpha, part_bodyCentered), \"s\", label = \"КОП\")\n",
    "ax.plot(alpha, pad(alpha, part_faceCentered), \"^\", label = \"КГЦ\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "#plt.ylabel(r\"\")\n",
    "plt.show()\n",
    "\n",
    "if savefig:\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(\"porosity-proportion.tiff\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANIAAAAQCAYAAABjuSH9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABJ0AAASdAHeZh94AAAISklEQVR4nO2bfbBWVRXGf1dBNFQySKmpVBgpIPXaB0IqQhqUgIFhNQ1MOgPkqCMK+JFlD8vGlAoKtQ9QRqyYmlIkFSWUmAi1nEkcbSAlPsWKBLuGApFAf6x94Nx9z7nvu897++8+M+/se9bZX+vZ66y999r7Nh08eJBOdKITjaFLLDCz9wG3Ap8GegJ/B5YAJulf9VZsZk3ApPAbCDQB64B7gfmSDkT5xwPnA83AmcBxwCJJExLanAD8NDxOlnRvSb5RwFRgAId1/BMwR9IzJWUuAK4GhgAnADuBF4G5kh5rRPdcuSTuzWwzcHJRXcB2Sb0b1d3MZgEfA/oBvYA9wJbQr7sl7SwoU3ks6+W5KsdRW+3ai5ldBtxXo5oDko48IirYFyf1cuBZ4HvARpz4Z8ysZ63O5fAzYD5wCvBzXMF3AD8CFhbk/zpOYDPwakI7Wd/fD9wNvFkj3yzgUeAjwDJgLvAc8FngqUBuXObbwJO4QT0MzAaWAu8GhhU0k6p7I9y/AVjB77sdoTtwHdAdeCLkXwS8DcwEXgi8x6g0lok8J3MctVWPvTxPMbcG/DbkeRzazkg/BE4ErpF0V67ROTihtwFX1NHJccCXgE3AIEk7gvwo4EFgopktkbQ4V+w6YBvwV9ybrazVTq69Jtxz7AQWAzNK8vUO77YDZ0j6Z+7dcJycW/FByuSTgeuB+4EpkvZFdXbtAN2hOvctkmYW6duo7gHHS9pbUN9twM3AV4Ero9fJY5nCcwMcZ+XrshdJz+MfU1Ed2ew9H+CI3Iu+wAhgM/CDuE7grdDB7kUVRxgX0tmZkqFj+4BbwuPVUadXSlovqcqm7Rrgk7g3f6udfCfjOv8xb0hZ+8Au3PsBYGbdcAPeSsHghnL/jUTJuncw92VI0j33rs1HFPDLkJ5WUCZpLCvwnMxxhHrtpay/pwOD8dl2KeQ+JGB4SJfH60tJu4Cn8KlzcB1tZWvzjQXvMtl5wYM0BDPrD9yBr6FX1ci+HtgHDDKzXlE9Q/G1/JM58adw41oMHDCzUWZ2o5lNNbMhJW1U0b0R7ruZ2QQzuzn0a7iZHVmQL1X3WhgT0hcSypQhlefK9pVoL2WYEtIFkvZD66XdB0P6cknh9bjX7AesqNFQ5iVOLXjXJ6Rdwt9/qVFXKcysC75Z3IovM9qFpNfN7EZgDrDWzJbg03tf4GJ8H/CVXJGPh3QvsAb4cNT+KmC8pNdy4iq6N8J9bw5vmDNsMrPLJf0uE1TQvRXMbAZwLNAD38Oci39Ed5SVSUAqz5XsK9VeimBmxwATgP34vgxoPSP1COkbJXVk8nfW0d7SkE4zs3flOtEV36hlOKGOutrDN4CzgMsk7amngKTvA5fgRE8GbgIuBV4BFkbLnhNDej1wEDgP99xnAMuBocCvoiaq6F6V+/uAC/CPqTtwOjAP34A/bmZn5jMn6h5jBr7MvBb/iJYBIyInUhWpPFe1r2R7KcDn8XFYJumVTNgm/N1B+AUwERiJe79f497mQuA9uEf4AFAzRFkGMzsb9yqzy0LWJeVuAL4F3IlHbf4BfAi4HVhkZs2SbgjZM0fzNnCxpM3h+cWw4X0JON/MhuT68H/XPYMki0R/Bq4wszeB6XhkLdtPpOoet9U71HES8Al8JlpjZqMlPdegKqk8J3Nc1V4KkC3r5hUpAIe9Xg+KkclbarUU1o1jcI/3GvDl8FuPD8KukLU9D1iKMEX/BF8K3VIje77cMGAW8LCkaZI2StodDGEcvnmcbmbZ8qAlpGtygwuApN3Ab8LjoJy8iu4dxn3Aj0M6NBNU0L0QkrZLeghfavbEx6FRtIS0Lp5TOa5qLzHMbGCofxvQ6uwwPyO9FNJ+JfVk0ZmydXwrhCjLrPDLd+boUNcOSZvqqasAx+b6udcsdswA3GNm9+CbymuDbHRI24RjJe02s2dxozoL37RmnLSU9CM7JD0mqitV9w7lHjcu8OVehlTd24WkLWa2Fmg2s1756FkFJPOcyHFVe4nRJsiQIT8jZQSPMLP4oPY44BxgN/CHkkbqxReBo/BDtKr4D7Cg5Lcm5FkdnvPTeLeQtgnzRvIs/LoCX7MPiDkJyDbF9TqEMt07mvssupf/IFJ1rwfvDen+dnPVRkfyXMRxVXs5hPCBTsR1XRC/PzQjSdpgZsvxKfsq4K5cPsO92zxJh+Lu4fyjK7AhPk8xs+Ml/TuSNQPfwT1M5WhP2ChOKnpnZjNxr3p/wRWh3+PnC1PMbJ6kV3PlPoMb7F7g6dDOFjN7BI9qTcVvG2T5R+Br9BZ8453vQ5LuFbnvD2zNy4L8FHz/A60PV5N0D/J++FWjVkGQYOzfxIMETxddX0pBFZ5TOG7AXvK4FA9ePJoPMmSIgw1X4kTeaX7naR1wNn7O8TLwtSj/Cvyg71T8MDGPJ8xsD74B3gX0B0bhd7XGSPpbpNBYYGx4zM4JhpjZwvD3DkmFJ9AJeAA/K7kQWGdmD+Eb7v740qcJuEmt749dhRM9x/ye2hpc37G4d5oUGxqJugekcv8FfE+zCr/7tgsPZY8CjsbX8PlrQlV0vwi43cxW47PBTuAk/LZCn1B+cqxIxbFM5bkKx40gW9bNL3rZahqVtAE/I1iID+J0fHDmAoNVcEGxHTyAhzAnANPwUOZ8YED+fCOHZg5vGkcGWZ+cbHxC24UIh50X4VdY1uJ7gun4UugxYKSkuVGZbcBHcS9/Gu4xhwGPAOdIerCgqVTdq3C/Er831xe/LjMNN/DVOF+j8zcEquiOf3gL8GXfJXh4+nPA6/hMOVDS2gJ1mkkcywo8J3NcFWH2P5eCIEOGps5/o+hEJxrH/wDJobz0mEzU2wAAAABJRU5ErkJggg==",
      "text/latex": [
       "$\\displaystyle 0.914869058316947$"
      ],
      "text/plain": [
       "0.9148690583169465"
      ]
     },
     "execution_count": 152,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "part_simple.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANIAAAAQCAYAAABjuSH9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABJ0AAASdAHeZh94AAAIPklEQVR4nO2bf7BVVRXHPw9RM1QysJh+TAkTjZD2sEJIUV4WlKSJYdM0kDQD1BgDCJhE0ddl4wQ1oK/oByDjM2NqTJERoSc/ZDLUdEafo82jJBASShIIQ4GIH/2x9oHt4dx7zzn3Tn+978yd/c7ee+211177x1pr79d0/PhxutCFLtSH7ukMM3sfcDvwWaAX8A9gBWCS/pW3YTMbBUwFBkTtPAsskPRUFbqrgMnAUOA8YA/wItAqaXU9NGY2BrgSaAY+CpwDLJM0tkK7vYDRwCjgIuC9wOHQ9j3APZKONVL+QDsWuC98TpR0d6p8fOBfDccknZaiKyR/oJkHfBzoD/QGDgLb8TmxUNKeVP0mYEL4DQSagE3A3cDiSuMVaAvrPqKtOmZl5c/Lo1uqYj9c2V8DngHuBLbiE+KpMLHyMJwHPAJcArQDrcBzwBeAJ0KHsuh+CKzDFfcwMB9YBZwPDG8AzXdxRTUDO3OIcgOwBLgUeBq4C3gQ+Ag+Me4PEyfdp1LyB9r3AwuBN6r063nAKvweC3V+l0FXVH6Am4EewNogxzLgCHAb8ELob4xfAYuBDwK/xsfp7cDPgbZKTMroPqLNM2ZQTv5cPNIn0s+AdwFTJP0kamQBPqB3AN+owbAPMBPYBVws6Z9RWQuu6NvxAY/pJgK3APcCkyQdTpWfnsGrKM3NwA7gr/jOtKGaLMBLwLXAqngnNbPZ+EbzReB6fHHVJX8ob8JPmj3A8tDOKZD0PL6YToGZJafd4oziovIDnCvpUAafO4DZwLeBm0LeaOArwMvAYEm7Q/4Z+BiNM7MVkpan2iqs+6gs15gFlJE/F49uUeV+wAhgG/DTVD0Bb+ID0aMG3w+Edp+OJxGApA3AfnyXiTt6Jr5I/0bGQAba/9ZLI2mDpM2ScjmGkh6TtDJtjkh6FfhF+ByeIissf4QpwKdwi+DNPH2MYWYXAUPw3XZVuryo/IHmlEUUcH9IPxTljQ7p/GQRhTYOA3PC5+RUnwvrMYXcY1ZG/rw8YtOuJaRrMibOfuAJ/IgeUoPpZtyPGGxmveMCM7sCt0vXpWg+g0+u5cAxMxtlZrea2VQzG1qBTxmaRiJR7pFUfhn5MbMLgbm4P/B4yT5NCulSSUdLtpEX14T0hSivT0i3ZtRP8oaFEypBaT02aMyqIi+P2LT7cEhfqlB3M35i9QfWV2pQ0l4zuxVYAHSa2Qr8SOyHm0lrga+nyD4R0kNAB+6DxMI8DoyR9FqdNA2BmXUHvho+2+OyMvKH9u7Dd+XZJft0FjAWOIr7JQ2Fmc0EzgZ64n7M5fgimhtVS06hCzKa6BvS7uHvP4fvUnpsxJjVQhEe8YnUM6SvV6ib5L+jVgck3YX7Dt2BicAs3HF/BWhLmzy4XwZuJx8HhuE798XAGuAK4LcNoGkU5uIKXy3p0XRhCfm/BwwCxks6WLJPX8J10y7plZJtVMNM3MSfhi+idmBEaoIn5uR0M3tnkhl8HIvqnRf9XVaPjRizWsjNo1u1wrIws28BD+BRmn541Odj+PG+LERosvpxBLhW0kZJb0h6Ebe7dwBXpo76MjSNkG0KMAPfUcdVqJNbfjO7FN/t5tcKi9dAYtYtqqONipDUR1ITbr5dj58qHWZ2SVTtN8CjuMydZrbIzFrxwMgwfGcHiF2Hwnps4JhVRFEe8UJKTpyeWRWj/H01OjAcmAc8LGm6pK2SDkh6Dh+YncAMM+sbkSVtdkjaFrcn6QCuHIDBddLUBTObjIeAO4EWSXsz6gwnp/zBdPglbk7PSbdVoF8DgU/ik67qfUu9kLRL0kO4md8L739SdhT3nWYBrwE3ht/m0L/9oWp8Iu8LaS49NmrMqqEMj9hH+ktI+1eom0RnKvlQCT4f0lNCi5IOmNkz+IQaxEkHNOG9r0KbyUXwWRn9LUJTGmY2Db9X+xNwVYZ5lqCI/Hs5Od6HzCxNArDEzJbgzu60Cjz/n0EGACRtN7NOoNnMeidRuhBhmxd+J2Bmb8Pn0G5JL0dFRfV4No0Zs2oozCM+kRLFjzCz9EXtOcBlwAHgjzU6cWZIK4V4k/w4zLket48HpHkHJA5orIAyNKUQggd34iZKS5VFBMXk/w+wtMKvI9TbGL4zzYswQcfhQYaltaVpKN4T0jyL98vAGfglbYyieqx7zHKgMI8TJ5KkLWa2Bj+yvwmcuJDFHcUewCJJJ+Lo4e7pdGBLFOv/A35XMMnMFknaGdX/HL4gDwFPRry3m9lKPKo1FZ+0Cc0IYCS+Y7XXQ1MGZjYHv0B9FneuTzHnUsgtf3BgJ1Tgext+at2b9dwlwg248/5Io4MMZtYf2CXp9VR+N+D7eKDgSUVPx8zsXEn/TtVvBn6Eny5xlK+wHhs0ZlVRhkf6ZcNN+AT/sfm7p03485gW3KT7Tqr+evwC8gL8IhfcyV4HfBrYZGYPAa8CF+JmTxMwS6k3WvjiHQQsMH+n1hHavQ7f8SakFVqUxsyuC2Vw8s5jqJm1hb93S5oZ1b8RX0RH8QUyJeOY3yapLfouK39ZJGZd1kuGt6Co/MDVwA/MbCN+IuwB3o2/CuiLyzUxxWatmR3ETeD9uNyj8Dd610j6e0bXyui+MErInxtvOUolbcHvCNrwBTQDj8C0AkPyKD9c5l6NP8foxP2BGfhF7mpgpKTWDLodeGRrIW5LT8VfDawELpP0YANomjnpAI8MeX2jvDGp+sl9yGl42FcZv/GNkL8MwmXh5eQPMjRTTP51uPlyPh6puwV/FrUXt1IGSupM0TyAh6/HAtPxMPZiYICk32d1qozuS6KZYvLnRlPXv1F0oQv143/2dYuALA3BlgAAAABJRU5ErkJggg==",
      "text/latex": [
       "$\\displaystyle 0.866128471396414$"
      ],
      "text/plain": [
       "0.8661284713964136"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "part_bodyCentered[~np.isnan(part_bodyCentered)].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMYAAAAQCAYAAABN/ABvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABJ0AAASdAHeZh94AAAHyUlEQVR4nO2ae7BXVRXHPxdFNERSUakpXyQGpMFUCKEgqWgBKpnVOOBjRsgxR5SXZmOLr44KliRKD1DGF46OD8RHaCgxkdljJq+TDaQMgoglCgZdA8KE/lj7wLn7nsPvnPO7P/+635k7+/72a33PXvux1tq7adeuXXSgAx1ojX3jDEmfAW4AzgIOBf4JLAJkZv8q2rGkJuDS8NcPaAJWAncD88xsZ72yJa0FjsqhsMHMehbgORZ4IPwcb2Z3R+UzgS8DvYEewDbgzcBrjpltiupfDNxTQ+xOM9snald6vEK7kcBEoC97xuwvwCwz+0O9vELbSnNC0mnAFcBg4GBgE/AqMNvMFtfDrdF6abUwJPUCXgIOB54E/g4MxAf+LElDYoF7wQLgAuBd4CFgK3AG8Avgq8CF7SR7C3B7Rv4HtQhK+iwwJ9Q9MKfa1cDLwPPhW7oCg4DpwARJg8zsrVT9VwDl9HUK8DXg2YyyUuMV+M8EpuETbhGwEfgccA5wnqQLzWxBPbyq6kXSrcBUYD3wVOB2GPAl4FRgcap6FW4N1Ut8YvwcH4ArzezO1EfOCkRuAi7L6Xw3JI3BlbwGGGhmG0P+fsDjwDhJi8xsYTvI3mxm02txyuDYhO8gm4CFwJScqgeZ2faM9jcB1wE/AC5P8s3sFVwJWTKTHXxelF96vCT1DJw3ACea2bupsuHAb/BdfkFVXgGl9SJpPL4o7gMmmNmOqLxz+ndFbg3VS6dUYS9gBLAW+FnU1oD/4ArqmtV5hDEhvS1RciC3A7g+/LyiQbKL4kp8l7gk9J+JrMEPeCSkxxURJukEfEd7G/hVVFxqvAKOwvX3p/SiCO2WAS34Dl2ZVxW9SOqCL5Z1ZCyKwO/DWrxqcWu0Xjql6gwP6ZLYnjWzFuD3wCdCJ7WQ2PZvZJQleaeEHbFe2V0kjZV0naSJkoZLamMnpyGpDzADt3WXF/ieLIwO6V8L1p8Q0vlm9lFUVna8AFYBO4CBknqkG0gaCnQDXqiTVxW9nIEvyIXATkkjJV0TdDO4AJ+i3PLQLnpJm1LHh/T1nA5W4btHb2BpDWHJrndMRtmxKdnH4jZrPbJ7ssd5TrBG0iVm9tu4I0n7hvrr8CO3ECRNwf2Q7rjTdzI++DMKtD0AGAt8hDvTMcqOF2b2vqRrgFnACkmLcLOwF3A2bnt/r05eVfTylZBuB5qBL0QylwPfMrP36uSW1GuIXtInRveQbsnpJ8n/ZC2B7DmSJkk6JEWkM60doIPrlH0PcBq+OLoCJwBzgaOBZyV9MaOvHwEDgIvNbFutD0lhCm4+XIUP/nPAiFoKDvh24P5c5BAmKDteAJjZ7cA38UUzHrgWOB94C7g3NrEq8Kqil8NDOhXYhTu23YATgSXAUODRGryKcEvQEL20Cde2Ex4GxgFn4rvZk/gOcjrwKXy3PhLIDEEWhZnFUYa/AZdJ+gCYjEcoEvsdSSfhp8Rt6VBmQVk9Qx9H4FGiGUCzpFFm9nKN5slxPTenvNJ4SZoG3AzcgUfX3gE+D9wCPCipv5lNq4NXFSSb7f+As81sbfj9aggyvAYMkzS4hg4KcWuUXtInRrL6u8eVovzNNYQRbLXR+A72HnBR+FsVyLeEqsmO1m6yA34Z0qFJRjCh7sfNguuzGhWBmW0wsydwE+LQ0GcuJPXDv3k9rUOU6T7LjheSTgVmAk+Z2SQze8PMtobJMAZ3JidLSkyx0ryoppfk/+bUoki+cyvw6/BzYE6fRbm1QnvrJX1ivBbS3jl9JV5+nr0ZE/0QV9zMiND+oa+NZramEbLxyQVuXiU4MNX/dikzpH2XpLtwp/yqvQkwszclrQD6S+qRjiZFKORAlhwvgFEhXZbR11ZJf8YXyACynfoivKroJWmzOadNciF4QE55UW6ZaC+9pE+MZIBHSErnI6kbMAS/dPpjGaIZ+C6wH36J1SjZSZQkPSH+C8zP+WsOdV4Mv4uaWZ8OaabywqQeF8rnF+wzRtZ4AXQJaV5INslvEy4twauKXpbivkXfuE1A4oyvyShrrzGrWy+7TwwzWy1pCX4UfR+4M1VP+O4718x2x/xDnLszsDqOTUs6yMz+HeX1B36M7xq7owYVZfcB1qXzQv7RuL0N4XIryNiGP7doA0nT8Z31vvSTEEm98aclW6L6nYAbcUfzpb08izgfd5ifqeFAlhqvgN/hdxsTJM01s7dT7b6OT9rt+K11JV5V9BJ27KfxyNhE4KcpXiNwP2oz7iRnoSa3j0MvsfN9OT6Qd8jfuawETsLj2a8DP4zqL8Uvmo7BL4HSeF7SNtwhbgH6ACPxNy2jzewfdcr+Dm5DL8ffyLTgocqRwP643fiTrI8ugW8At0h6Ed/hNgFHAMPw0Ok7eDQoD8lxnXWjHKPseD2G31OcDqyU9ETg0wc3s5qAa3Oe0ZThVVYv4ItoADBL/parGZ8j5+K79KXxpC7JreF6aXXUmdlqPBZ8L/7xk/HJNhsYVOKdFLjiuuFx4kl4uG4e0DfrfqGC7GXAM6HOBUHGMNwcuggYlXXrWhIv4EftYXhYdCpwHvA+vmP2M7MVWQ3DiXYyxR3IsuO1E58gVwMrcH9iMm5GLgbONLPZ9fKqMifMbD3+JmoO7odMxN9HPQ0MMbPHs2SV4NZwvTR1PDvvQAfa4v+O5k7D0lMsqwAAAABJRU5ErkJggg==",
      "text/latex": [
       "$\\displaystyle 0.90543798706737$"
      ],
      "text/plain": [
       "0.90543798706737"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "part_faceCentered[~np.isnan(part_faceCentered)].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "19649669bd52b0be75e091dcf60d2128e4a347083ff474cfec5ff9275df3ceed"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}