Flow_curves_models_fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0747f869",
   "metadata": {},
   "source": [
    "# Theoretical models fitted to flow curves\n",
    "\n",
    "In this notebook I use the curve fitting function from Python SK-Learn package to apply rheological models to experimental data. The input data consist of a csv file with two columns corresponding to stress and strain rate measurements. \n",
    "\n",
    "This software calculates viscosity and fits several models to the flow curves. \n",
    "\n",
    "The output consists of two spreadsheets: \n",
    "- experimental data and stress values obtained from each fit\n",
    "- scores for each fit and relative parameters values\n",
    "\n",
    "Those shown here as an example are synthetic data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c2be77d",
   "metadata": {},
   "source": [
    "## 1. import libraries and data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b866c6b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Import curve fitting package from scipy\n",
    "from scipy.optimize import curve_fit\n",
    "from sklearn.metrics import r2_score, explained_variance_score, mean_absolute_error, mean_squared_error\n",
    "from pandas import ExcelWriter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5842d1e3",
   "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>Shear rate</th>\n",
       "      <th>Shear stress</th>\n",
       "      <th>Viscosity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.000000</td>\n",
       "      <td>66.044951</td>\n",
       "      <td>33.022476</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.472529</td>\n",
       "      <td>71.833355</td>\n",
       "      <td>29.052583</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.054742</td>\n",
       "      <td>79.110979</td>\n",
       "      <td>25.897763</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.978503</td>\n",
       "      <td>88.388858</td>\n",
       "      <td>22.216613</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.998865</td>\n",
       "      <td>98.365375</td>\n",
       "      <td>19.677543</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Shear rate  Shear stress  Viscosity\n",
       "0    2.000000     66.044951  33.022476\n",
       "1    2.472529     71.833355  29.052583\n",
       "2    3.054742     79.110979  25.897763\n",
       "3    3.978503     88.388858  22.216613\n",
       "4    4.998865     98.365375  19.677543"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# import from csv containing stres and strain rate and delete eventual empty cells\n",
    "df_ares = pd.read_csv('input_data.csv').dropna() \n",
    "#define viscosity. We calculate it rather than importing it from the instrument software\n",
    "df_ares['Viscosity'] =  df_ares['Shear stress']/df_ares['Shear rate'] \n",
    "df_ares.head() #visualize input table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "842f4f8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract numpy arrays from the dataframe columns\n",
    "x1 = df_ares['Shear rate'].values\n",
    "y_stress = df_ares['Shear stress'].values\n",
    "y_visc = df_ares['Viscosity'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a68976e",
   "metadata": {},
   "source": [
    "## 2. Modelling : equation coefficients calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8a1c3d5",
   "metadata": {},
   "source": [
    "equations from:\n",
    "https://blog.rheosense.com/modeling-non-newtonian-fluids"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59f691a4",
   "metadata": {},
   "source": [
    "### 2.1 Models for non-newtonian fluids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1410437e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to calculate the Ostwald de Waele (power law) model \n",
    "def power(x,  K, n):\n",
    "    return K*x**(n-1)\n",
    "\n",
    "# Function to calculate the Carreau model\n",
    "# from https://en.wikipedia.org/wiki/Carreau_fluid\n",
    "def Carreau(x ,eta0, eta_inf, L, n):\n",
    "    return eta_inf + (eta0 - eta_inf)*(1 + (L*x)**2)**((n-1)/2)\n",
    "\n",
    "# Function to calculate the Yasuda model\n",
    "def Yasuda(x ,eta0, eta_inf, L, n,a):\n",
    "    return (eta0-eta_inf) * ((1 + (L*x)**a)**((n-1)/a))\n",
    "\n",
    "# Function to calculate the Cross model\n",
    "def Cross(x ,eta0, eta_inf, k, n):\n",
    "    return eta_inf + (eta0-eta_inf)/(1 +(k*x)**n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fc437e81",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([41.55,  0.57]),\n",
       " array([72.59,  7.06, -1.44,  0.17]),\n",
       " array([ 20.68, -35.49,   2.02,   0.57,  11.11]),\n",
       " array([211.08,   7.15,   4.57,   0.87])]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def non_newtonian_fits(x,y):\n",
    "    non_newtonian_models = [power, Carreau, Yasuda, Cross]\n",
    "    model_param_list = []\n",
    "    for model in non_newtonian_models:\n",
    "        param, param_cov = curve_fit(model, x,y , maxfev=5000)\n",
    "        model_param_list.append(np.round(param,2))\n",
    "    return model_param_list\n",
    "\n",
    "NonNewtonian_param_list = non_newtonian_fits(x1, y_visc)\n",
    "\n",
    "NonNewtonian_param_list "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5da91e0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#param1, param_cov = curve_fit(Carreau, x1, y_visc , maxfev=5000)\n",
    "#NonNewtonian_param_list[1] = param1\n",
    "#NonNewtonian_param_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa59434d",
   "metadata": {},
   "source": [
    "### 2.2 Models that assume that the fluid has a yield stress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "48fa67a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to calculate the Bingham model\n",
    "def Bingham(x ,y0, pv):\n",
    "    return y0 + pv*x\n",
    "# https://glossary.oilfield.slb.com/en/terms/b/bingham_plastic_model\n",
    " \n",
    "# Function to calculate the HB model\n",
    "def Herschel_Bulkley(x ,y0, K, n):\n",
    "    return (y0 + K*x**n)\n",
    "\n",
    "# Function to calculate the HB model\n",
    "def Casson(x ,yc, eta):\n",
    "    return ((np.sqrt(yc) + np.sqrt(eta*x))**2)/1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8c4c7c1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def yield_stress_fits(x,y):\n",
    "    yield_stress_models = [Bingham, Herschel_Bulkley, Casson]\n",
    "    model_param_list = []\n",
    "    for model in yield_stress_models:\n",
    "        param, param_cov = curve_fit(model, x,y , maxfev=5000)\n",
    "        model_param_list.append(np.round(param,2))\n",
    "    return model_param_list\n",
    "\n",
    "yield_stress_param_list = yield_stress_fits(x1, y_stress)\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4601cc97",
   "metadata": {},
   "source": [
    "### 3. Modelling : viscosity values from each fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "27166902",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now let's calculate the viscosity values from the fit\n",
    "eta_pow = NonNewtonian_param_list[0][0] * (x1**(NonNewtonian_param_list[0][1]-1))\n",
    "\n",
    "# Carreau parameters\n",
    "\n",
    "eta0 = NonNewtonian_param_list[1][0]\n",
    "eta_inf = NonNewtonian_param_list[1][1]\n",
    "L = NonNewtonian_param_list[1][2]\n",
    "n = NonNewtonian_param_list[1][3]\n",
    "\n",
    "\n",
    "eta_car = eta_inf + (eta0 - eta_inf)*(1 + (L*x1)**2)**((n-1)/2)\n",
    "\n",
    "# Yasuda parameters\n",
    "eta0 = NonNewtonian_param_list[2][0]\n",
    "eta_inf = NonNewtonian_param_list[2][1]\n",
    "L = NonNewtonian_param_list[2][2]\n",
    "n = NonNewtonian_param_list[2][3]\n",
    "a = NonNewtonian_param_list[2][4]\n",
    "\n",
    "eta_yas = (eta0-eta_inf) * ((1 + (L*x1)**a)**((n-1)/a))\n",
    "#cross parameters\n",
    "c_eta0 = NonNewtonian_param_list[3][0]\n",
    "c_eta_inf = NonNewtonian_param_list[3][1]\n",
    "c_k = NonNewtonian_param_list[3][2]\n",
    "c_n = NonNewtonian_param_list[3][3]\n",
    "\n",
    "eta_cross = c_eta_inf + (c_eta0-c_eta_inf)/(1 +(c_k*x1)**c_n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4f24279a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now let's calculate the stress and viscosity values from the fit\n",
    "\n",
    "y_bin = yield_stress_param_list[0][0] + yield_stress_param_list[0][1]*x1\n",
    "eta_bin = y_bin/x1\n",
    "\n",
    "y_hb =   yield_stress_param_list[1][0] + yield_stress_param_list[1][1]*x1**yield_stress_param_list[1][2]\n",
    "eta_hb = y_hb/x1\n",
    "\n",
    "y_cas = (np.sqrt(yield_stress_param_list[2][0]) + np.sqrt(yield_stress_param_list[2][1]*x1))**2\n",
    "eta_cas = y_cas/x1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2475312",
   "metadata": {},
   "source": [
    "### 4. Models scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "403bdab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_models(eta_list):    \n",
    "    MLA_names = []\n",
    "    ExplVar_scores = [] # explained variance list\n",
    "    MAE_scores = [] # mean average error list\n",
    "    RMSE_scores = [] # root mean square error list\n",
    "    R2_scores = [] # regression coefficients list\n",
    "    MLA_names = ['Ostwald – de Waele power law', 'Carreau', 'Carreau-Yasuda', 'Cross', 'Bingham',\"Herschel-Bulkley\", 'Casson']\n",
    "       \n",
    "    for y_pred in eta_list:\n",
    "                \n",
    "        #model scores  \n",
    "        R2 = r2_score(y_visc, y_pred)\n",
    "        EV = explained_variance_score(y_visc, y_pred)\n",
    "        MAE = mean_absolute_error(y_visc, y_pred)\n",
    "\n",
    "        MSE = mean_squared_error(y_visc, y_pred)\n",
    "        RMSE = np.sqrt(MSE)\n",
    "\n",
    "        # add results to lists for a final comparison\n",
    "               \n",
    "        MAE_scores.append(round(MAE, 2))\n",
    "        ExplVar_scores.append(round(EV, 2))\n",
    "        RMSE_scores.append(round(RMSE, 2))\n",
    "        R2_scores.append(round(R2, 4))\n",
    "        \n",
    "    \n",
    "    #create table to compare MLA metrics\n",
    "    MLA_columns = ['Model Name', 'Explained Variance','MAE', 'RMSE', 'R2']\n",
    "    zippedList =  list(zip(MLA_names, ExplVar_scores, MAE_scores, RMSE_scores, R2_scores))\n",
    "    df = pd.DataFrame(zippedList, columns = MLA_columns)#.sort_values(by=['R2'],ascending=False)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c6ce9f03",
   "metadata": {},
   "outputs": [],
   "source": [
    "eta_list = [eta_pow, eta_car, eta_yas, eta_cross ,eta_bin, eta_hb, eta_cas]\n",
    "df_results = compare_models(eta_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "dae285e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "list_power =['K','n']\n",
    "list_carreau = ['eta0',  'eta_inf','L', 'n']\n",
    "list_yasuda = ['eta0', 'eta_inf', 'L', 'n','a']\n",
    "list_cross = ['eta0', 'eta_inf', 'k', 'n']\n",
    "list_bingham = ['yield_stress', 'plastic viscosity']\n",
    "list_hb = ['yield_stress', 'K', 'n']\n",
    "list_casson = ['yield_stress','casson_viscosity']\n",
    "\n",
    "models_name_list = ['Ostwald – de Waele power law', 'Carreau', 'Carreau-Yasuda', 'Cross', 'Bingham',\"Herschel-Bulkley\", 'Casson']\n",
    "coefficient_list = [list_power, list_carreau, list_yasuda, list_cross, list_bingham, list_hb, list_casson]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "332a3129",
   "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>Shear rate</th>\n",
       "      <th>Shear stress</th>\n",
       "      <th>Viscosity</th>\n",
       "      <th>Ostwald – de Waele power law_viscosity</th>\n",
       "      <th>Carreau_viscosity</th>\n",
       "      <th>Carreau-Yasuda_viscosity</th>\n",
       "      <th>Cross_viscosity</th>\n",
       "      <th>Bingham_viscosity</th>\n",
       "      <th>Herschel-Bulkley_viscosity</th>\n",
       "      <th>Casson_viscosity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.000000</td>\n",
       "      <td>66.044951</td>\n",
       "      <td>33.022476</td>\n",
       "      <td>30.840977</td>\n",
       "      <td>33.039351</td>\n",
       "      <td>30.814877</td>\n",
       "      <td>33.110987</td>\n",
       "      <td>39.035000</td>\n",
       "      <td>35.145331</td>\n",
       "      <td>29.594894</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.472529</td>\n",
       "      <td>71.833355</td>\n",
       "      <td>29.052583</td>\n",
       "      <td>28.152711</td>\n",
       "      <td>29.191320</td>\n",
       "      <td>28.128886</td>\n",
       "      <td>29.209794</td>\n",
       "      <td>33.021673</td>\n",
       "      <td>30.201274</td>\n",
       "      <td>26.358153</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.054742</td>\n",
       "      <td>79.110979</td>\n",
       "      <td>25.897763</td>\n",
       "      <td>25.705851</td>\n",
       "      <td>25.826734</td>\n",
       "      <td>25.684097</td>\n",
       "      <td>25.842758</td>\n",
       "      <td>28.170759</td>\n",
       "      <td>26.186923</td>\n",
       "      <td>23.614809</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.978503</td>\n",
       "      <td>88.388858</td>\n",
       "      <td>22.216613</td>\n",
       "      <td>22.945185</td>\n",
       "      <td>22.259271</td>\n",
       "      <td>22.925767</td>\n",
       "      <td>22.289029</td>\n",
       "      <td>23.387508</td>\n",
       "      <td>22.193237</td>\n",
       "      <td>20.746951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.998865</td>\n",
       "      <td>98.365375</td>\n",
       "      <td>19.677543</td>\n",
       "      <td>20.799665</td>\n",
       "      <td>19.692546</td>\n",
       "      <td>20.782062</td>\n",
       "      <td>19.730071</td>\n",
       "      <td>20.158858</td>\n",
       "      <td>19.466367</td>\n",
       "      <td>18.681347</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6.246547</td>\n",
       "      <td>109.289774</td>\n",
       "      <td>17.496029</td>\n",
       "      <td>18.899304</td>\n",
       "      <td>17.589413</td>\n",
       "      <td>18.883310</td>\n",
       "      <td>17.627105</td>\n",
       "      <td>17.644365</td>\n",
       "      <td>17.315045</td>\n",
       "      <td>16.968808</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>7.945852</td>\n",
       "      <td>124.144791</td>\n",
       "      <td>15.623849</td>\n",
       "      <td>17.041610</td>\n",
       "      <td>15.699951</td>\n",
       "      <td>17.027188</td>\n",
       "      <td>15.731454</td>\n",
       "      <td>15.489855</td>\n",
       "      <td>15.441708</td>\n",
       "      <td>15.399278</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>9.967245</td>\n",
       "      <td>141.157305</td>\n",
       "      <td>14.162118</td>\n",
       "      <td>15.459083</td>\n",
       "      <td>14.226534</td>\n",
       "      <td>15.446000</td>\n",
       "      <td>14.249140</td>\n",
       "      <td>13.883680</td>\n",
       "      <td>14.017062</td>\n",
       "      <td>14.142594</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>12.609521</td>\n",
       "      <td>163.332623</td>\n",
       "      <td>12.953119</td>\n",
       "      <td>13.972388</td>\n",
       "      <td>12.960294</td>\n",
       "      <td>13.960563</td>\n",
       "      <td>12.973237</td>\n",
       "      <td>12.560673</td>\n",
       "      <td>12.815265</td>\n",
       "      <td>13.028496</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>15.820988</td>\n",
       "      <td>189.006329</td>\n",
       "      <td>11.946557</td>\n",
       "      <td>12.673617</td>\n",
       "      <td>11.949771</td>\n",
       "      <td>12.662892</td>\n",
       "      <td>11.954698</td>\n",
       "      <td>11.547628</td>\n",
       "      <td>11.868100</td>\n",
       "      <td>12.107372</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>19.944853</td>\n",
       "      <td>221.115647</td>\n",
       "      <td>11.086351</td>\n",
       "      <td>11.472122</td>\n",
       "      <td>11.095723</td>\n",
       "      <td>11.462414</td>\n",
       "      <td>11.094737</td>\n",
       "      <td>10.725200</td>\n",
       "      <td>11.072311</td>\n",
       "      <td>11.298162</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>25.081726</td>\n",
       "      <td>262.393593</td>\n",
       "      <td>10.461545</td>\n",
       "      <td>10.395547</td>\n",
       "      <td>10.397295</td>\n",
       "      <td>10.386749</td>\n",
       "      <td>10.393031</td>\n",
       "      <td>10.078998</td>\n",
       "      <td>10.420774</td>\n",
       "      <td>10.607904</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>31.643086</td>\n",
       "      <td>309.887889</td>\n",
       "      <td>9.793226</td>\n",
       "      <td>9.407000</td>\n",
       "      <td>9.812207</td>\n",
       "      <td>9.399039</td>\n",
       "      <td>9.807147</td>\n",
       "      <td>9.558744</td>\n",
       "      <td>9.869992</td>\n",
       "      <td>10.002948</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>39.909527</td>\n",
       "      <td>370.705635</td>\n",
       "      <td>9.288650</td>\n",
       "      <td>8.513499</td>\n",
       "      <td>9.330132</td>\n",
       "      <td>8.506294</td>\n",
       "      <td>9.326493</td>\n",
       "      <td>9.146816</td>\n",
       "      <td>9.407960</td>\n",
       "      <td>9.479894</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>50.143038</td>\n",
       "      <td>446.863687</td>\n",
       "      <td>8.911779</td>\n",
       "      <td>7.717568</td>\n",
       "      <td>8.938408</td>\n",
       "      <td>7.711036</td>\n",
       "      <td>8.937913</td>\n",
       "      <td>8.825010</td>\n",
       "      <td>9.022106</td>\n",
       "      <td>9.032958</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>63.146415</td>\n",
       "      <td>541.836191</td>\n",
       "      <td>8.580633</td>\n",
       "      <td>6.989097</td>\n",
       "      <td>8.611272</td>\n",
       "      <td>6.983183</td>\n",
       "      <td>8.615274</td>\n",
       "      <td>8.566573</td>\n",
       "      <td>8.687564</td>\n",
       "      <td>8.639655</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>79.558602</td>\n",
       "      <td>661.748097</td>\n",
       "      <td>8.317744</td>\n",
       "      <td>6.328132</td>\n",
       "      <td>8.340605</td>\n",
       "      <td>6.322777</td>\n",
       "      <td>8.350033</td>\n",
       "      <td>8.360989</td>\n",
       "      <td>8.396980</td>\n",
       "      <td>8.295914</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>100.000000</td>\n",
       "      <td>811.304600</td>\n",
       "      <td>8.113046</td>\n",
       "      <td>5.735497</td>\n",
       "      <td>8.119231</td>\n",
       "      <td>5.730643</td>\n",
       "      <td>8.134583</td>\n",
       "      <td>8.199300</td>\n",
       "      <td>8.144699</td>\n",
       "      <td>7.998423</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Shear rate  Shear stress  Viscosity  \\\n",
       "0     2.000000     66.044951  33.022476   \n",
       "1     2.472529     71.833355  29.052583   \n",
       "2     3.054742     79.110979  25.897763   \n",
       "3     3.978503     88.388858  22.216613   \n",
       "4     4.998865     98.365375  19.677543   \n",
       "5     6.246547    109.289774  17.496029   \n",
       "6     7.945852    124.144791  15.623849   \n",
       "7     9.967245    141.157305  14.162118   \n",
       "8    12.609521    163.332623  12.953119   \n",
       "9    15.820988    189.006329  11.946557   \n",
       "10   19.944853    221.115647  11.086351   \n",
       "11   25.081726    262.393593  10.461545   \n",
       "12   31.643086    309.887889   9.793226   \n",
       "13   39.909527    370.705635   9.288650   \n",
       "14   50.143038    446.863687   8.911779   \n",
       "15   63.146415    541.836191   8.580633   \n",
       "16   79.558602    661.748097   8.317744   \n",
       "17  100.000000    811.304600   8.113046   \n",
       "\n",
       "    Ostwald – de Waele power law_viscosity  Carreau_viscosity  \\\n",
       "0                                30.840977          33.039351   \n",
       "1                                28.152711          29.191320   \n",
       "2                                25.705851          25.826734   \n",
       "3                                22.945185          22.259271   \n",
       "4                                20.799665          19.692546   \n",
       "5                                18.899304          17.589413   \n",
       "6                                17.041610          15.699951   \n",
       "7                                15.459083          14.226534   \n",
       "8                                13.972388          12.960294   \n",
       "9                                12.673617          11.949771   \n",
       "10                               11.472122          11.095723   \n",
       "11                               10.395547          10.397295   \n",
       "12                                9.407000           9.812207   \n",
       "13                                8.513499           9.330132   \n",
       "14                                7.717568           8.938408   \n",
       "15                                6.989097           8.611272   \n",
       "16                                6.328132           8.340605   \n",
       "17                                5.735497           8.119231   \n",
       "\n",
       "    Carreau-Yasuda_viscosity  Cross_viscosity  Bingham_viscosity  \\\n",
       "0                  30.814877        33.110987          39.035000   \n",
       "1                  28.128886        29.209794          33.021673   \n",
       "2                  25.684097        25.842758          28.170759   \n",
       "3                  22.925767        22.289029          23.387508   \n",
       "4                  20.782062        19.730071          20.158858   \n",
       "5                  18.883310        17.627105          17.644365   \n",
       "6                  17.027188        15.731454          15.489855   \n",
       "7                  15.446000        14.249140          13.883680   \n",
       "8                  13.960563        12.973237          12.560673   \n",
       "9                  12.662892        11.954698          11.547628   \n",
       "10                 11.462414        11.094737          10.725200   \n",
       "11                 10.386749        10.393031          10.078998   \n",
       "12                  9.399039         9.807147           9.558744   \n",
       "13                  8.506294         9.326493           9.146816   \n",
       "14                  7.711036         8.937913           8.825010   \n",
       "15                  6.983183         8.615274           8.566573   \n",
       "16                  6.322777         8.350033           8.360989   \n",
       "17                  5.730643         8.134583           8.199300   \n",
       "\n",
       "    Herschel-Bulkley_viscosity  Casson_viscosity  \n",
       "0                    35.145331         29.594894  \n",
       "1                    30.201274         26.358153  \n",
       "2                    26.186923         23.614809  \n",
       "3                    22.193237         20.746951  \n",
       "4                    19.466367         18.681347  \n",
       "5                    17.315045         16.968808  \n",
       "6                    15.441708         15.399278  \n",
       "7                    14.017062         14.142594  \n",
       "8                    12.815265         13.028496  \n",
       "9                    11.868100         12.107372  \n",
       "10                   11.072311         11.298162  \n",
       "11                   10.420774         10.607904  \n",
       "12                    9.869992         10.002948  \n",
       "13                    9.407960          9.479894  \n",
       "14                    9.022106          9.032958  \n",
       "15                    8.687564          8.639655  \n",
       "16                    8.396980          8.295914  \n",
       "17                    8.144699          7.998423  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_final = df_ares\n",
    "\n",
    "for i in range(len(models_name_list)):\n",
    "    name = \"%s_viscosity\" % models_name_list[i]\n",
    "    df_final[name] = eta_list[i]\n",
    "\n",
    "df_final"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48eb4e6b",
   "metadata": {},
   "source": [
    "## 5. Output files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "cdca09b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter output file name: try\n",
      "try\n"
     ]
    }
   ],
   "source": [
    "output_name = input(\"Enter output file name: \")\n",
    "print(output_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "08d20b12",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([41.55,  0.57]),\n",
       " array([72.59,  7.06, -1.44,  0.17]),\n",
       " array([ 20.68, -35.49,   2.02,   0.57,  11.11]),\n",
       " array([211.08,   7.15,   4.57,   0.87]),\n",
       " array([62.93,  7.57]),\n",
       " array([50.19, 10.55,  0.93]),\n",
       " array([18.51,  5.75])]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l1 = NonNewtonian_param_list\n",
    "l2 = yield_stress_param_list\n",
    "\n",
    "parameters_values = [*l1, *l2]#.append()\n",
    "parameters_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a13dd1e1",
   "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>model</th>\n",
       "      <th>coefficient</th>\n",
       "      <th>value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Ostwald – de Waele power law</td>\n",
       "      <td>[K, n]</td>\n",
       "      <td>[41.55, 0.57]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Carreau</td>\n",
       "      <td>[eta0, eta_inf, L, n]</td>\n",
       "      <td>[72.59, 7.06, -1.44, 0.17]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Carreau-Yasuda</td>\n",
       "      <td>[eta0, eta_inf, L, n, a]</td>\n",
       "      <td>[20.68, -35.49, 2.02, 0.57, 11.11]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Cross</td>\n",
       "      <td>[eta0, eta_inf, k, n]</td>\n",
       "      <td>[211.08, 7.15, 4.57, 0.87]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Bingham</td>\n",
       "      <td>[yield_stress, plastic viscosity]</td>\n",
       "      <td>[62.93, 7.57]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Herschel-Bulkley</td>\n",
       "      <td>[yield_stress, K, n]</td>\n",
       "      <td>[50.19, 10.55, 0.93]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Casson</td>\n",
       "      <td>[yield_stress, casson_viscosity]</td>\n",
       "      <td>[18.51, 5.75]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          model                        coefficient  \\\n",
       "0  Ostwald – de Waele power law                             [K, n]   \n",
       "1                       Carreau              [eta0, eta_inf, L, n]   \n",
       "2                Carreau-Yasuda           [eta0, eta_inf, L, n, a]   \n",
       "3                         Cross              [eta0, eta_inf, k, n]   \n",
       "4                       Bingham  [yield_stress, plastic viscosity]   \n",
       "5              Herschel-Bulkley               [yield_stress, K, n]   \n",
       "6                        Casson   [yield_stress, casson_viscosity]   \n",
       "\n",
       "                                value  \n",
       "0                       [41.55, 0.57]  \n",
       "1          [72.59, 7.06, -1.44, 0.17]  \n",
       "2  [20.68, -35.49, 2.02, 0.57, 11.11]  \n",
       "3          [211.08, 7.15, 4.57, 0.87]  \n",
       "4                       [62.93, 7.57]  \n",
       "5                [50.19, 10.55, 0.93]  \n",
       "6                       [18.51, 5.75]  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_zippedList =  list(zip(models_name_list, coefficient_list, parameters_values))\n",
    "df_param = pd.DataFrame(param_zippedList, columns = ['model','coefficient','value'])\n",
    "df_param"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "b1e73a20",
   "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>model</th>\n",
       "      <th>coefficient</th>\n",
       "      <th>value</th>\n",
       "      <th>Explained Variance</th>\n",
       "      <th>MAE</th>\n",
       "      <th>RMSE</th>\n",
       "      <th>R2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Ostwald – de Waele power law</td>\n",
       "      <td>[K, n]</td>\n",
       "      <td>[41.55, 0.57]</td>\n",
       "      <td>0.97</td>\n",
       "      <td>1.10</td>\n",
       "      <td>1.27</td>\n",
       "      <td>0.9709</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Carreau</td>\n",
       "      <td>[eta0, eta_inf, L, n]</td>\n",
       "      <td>[72.59, 7.06, -1.44, 0.17]</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.9999</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Carreau-Yasuda</td>\n",
       "      <td>[eta0, eta_inf, L, n, a]</td>\n",
       "      <td>[20.68, -35.49, 2.02, 0.57, 11.11]</td>\n",
       "      <td>0.97</td>\n",
       "      <td>1.10</td>\n",
       "      <td>1.27</td>\n",
       "      <td>0.9709</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Cross</td>\n",
       "      <td>[eta0, eta_inf, k, n]</td>\n",
       "      <td>[211.08, 7.15, 4.57, 0.87]</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.07</td>\n",
       "      <td>0.9999</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Bingham</td>\n",
       "      <td>[yield_stress, plastic viscosity]</td>\n",
       "      <td>[62.93, 7.57]</td>\n",
       "      <td>0.95</td>\n",
       "      <td>0.92</td>\n",
       "      <td>1.82</td>\n",
       "      <td>0.9405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Herschel-Bulkley</td>\n",
       "      <td>[yield_stress, K, n]</td>\n",
       "      <td>[50.19, 10.55, 0.93]</td>\n",
       "      <td>0.99</td>\n",
       "      <td>0.28</td>\n",
       "      <td>0.58</td>\n",
       "      <td>0.9939</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Casson</td>\n",
       "      <td>[yield_stress, casson_viscosity]</td>\n",
       "      <td>[18.51, 5.75]</td>\n",
       "      <td>0.98</td>\n",
       "      <td>0.72</td>\n",
       "      <td>1.25</td>\n",
       "      <td>0.9721</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          model                        coefficient  \\\n",
       "0  Ostwald – de Waele power law                             [K, n]   \n",
       "1                       Carreau              [eta0, eta_inf, L, n]   \n",
       "2                Carreau-Yasuda           [eta0, eta_inf, L, n, a]   \n",
       "3                         Cross              [eta0, eta_inf, k, n]   \n",
       "4                       Bingham  [yield_stress, plastic viscosity]   \n",
       "5              Herschel-Bulkley               [yield_stress, K, n]   \n",
       "6                        Casson   [yield_stress, casson_viscosity]   \n",
       "\n",
       "                                value  Explained Variance   MAE  RMSE      R2  \n",
       "0                       [41.55, 0.57]                0.97  1.10  1.27  0.9709  \n",
       "1          [72.59, 7.06, -1.44, 0.17]                1.00  0.04  0.05  0.9999  \n",
       "2  [20.68, -35.49, 2.02, 0.57, 11.11]                0.97  1.10  1.27  0.9709  \n",
       "3          [211.08, 7.15, 4.57, 0.87]                1.00  0.06  0.07  0.9999  \n",
       "4                       [62.93, 7.57]                0.95  0.92  1.82  0.9405  \n",
       "5                [50.19, 10.55, 0.93]                0.99  0.28  0.58  0.9939  \n",
       "6                       [18.51, 5.75]                0.98  0.72  1.25  0.9721  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_results = pd.merge(df_param, df_results, left_on='model', right_on='Model Name')\n",
    "df_results = df_results.drop(['Model Name'], axis=1)\n",
    "df_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "5fb167aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Pandas Excel writer using XlsxWriter as the engine.\n",
    "writer = pd.ExcelWriter(\"%s_output.xlsx\" % output_name, engine='xlsxwriter')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f1e69f1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the dataframe to an XlsxWriter Excel object.\n",
    "df_final.to_excel(writer, sheet_name='data')\n",
    "df_results.to_excel(writer, sheet_name='results')\n",
    "\n",
    "# Close the Pandas Excel writer and output the Excel file.\n",
    "\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec135b82",
   "metadata": {},
   "source": [
    "Finally let's plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a87394a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Visualize the chosen model fit\n",
      "0 : Ostwald – de Waele power law\n",
      "1 : Carreau\n",
      "2 : Carreau-Yasuda\n",
      "3 : Cross\n",
      "4 : Bingham\n",
      "5 : Herschel-Bulkley\n",
      "6 : Casson\n"
     ]
    }
   ],
   "source": [
    "print('Visualize the chosen model fit')\n",
    "for i in range(len(models_name_list)):\n",
    "    print(i,':',models_name_list[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "09110fc9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter the number corresponding to the model: 1\n"
     ]
    }
   ],
   "source": [
    "model_number = input(\"Enter the number corresponding to the model: \")\n",
    "num = int(model_number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "7b74c695",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Shear Rate [1/s]')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAs+0lEQVR4nO3deXwV5dn/8c+VAIYAIiBSEJNQZXFhC7jiUrXWtS6ttk9NFYuaVmvdWrcntbX9NU/RWrdqbWOxokRxQ2vtIq5VcaGAiAoIUpOIRXYQCEtIrt8fMyeGkOUknJOTnPm+X695nTNz5sxckwPXzNz3Pfdt7o6IiERHRqoDEBGRtqXELyISMUr8IiIRo8QvIhIxSvwiIhGjxC8iEjFK/CIiEaPEL5FhZuea2Swz22hmy8zsH2Z2ZKrjEmlrSvwSCWZ2NXAH8H9APyAH+D1wRgu30ymeZSLtmRK/pD0z6wn8Evihu09z903uXuXuf3X3a8zsEDN708zWhXcCd5tZlzrfdzP7oZktBhab2VfMbKmZXWdmnwF/NrMMM7vezJaY2Woze8zMetfZxuNm9pmZrTezV83swDqfvWJmF9WZv8DMXm+bv45EkRK/RMHhQBbwVCOfVwNXAXuG6x4PXFpvnTOBQ4EDwvkvAb2BXKAQ+FG4zjHAAGAtcE+d7/8DGAzsBcwBSlt/OCK7RolfoqAPsMrdtzf0obvPdve33H27u5cBfyRI4HX92t3XuPvmcL4G+Lm7bw2X/QAocvel7r4VuAk4O1YM5O73u/uGOp+NDO9ERNqcyiYlClYDe5pZp4aSv5kNAW4DxgLZBP8vZtdb7ZN68yvdfUud+VzgKTOrqbOsGugXFgcVA+cAfQlOGhDcYaxv3SGJtJ6u+CUK3gS2EhTFNOReYCEw2N13B/4XsHrr1O/Gtv78J8DJ7r5HnSnL3T8FziWoRP4q0BPIC78T28cmghNOzJfiOSiR1lLil7Tn7uuBnwH3mNmZZpZtZp3N7GQzuwXoAXwObDSzYcAlrdjNH4BiM8sFMLO+ZhZrMdSD4MSzmiDB/1+9784FvhHGtR9wYSv2LxI3JX6JBHf/LXA18FNgJcEV+mXA08BPCK7KNwD3AY+2Yhd3As8A081sA/AWQWUwwINAOfApMD/8rK7bgW3AcmAyqviVJDMNxCIiEi264hcRiRglfhGRiFHiFxGJGCV+EZGI6RAPcO25556el5eX6jBERDqU2bNnr3L3vvWXd4jEn5eXx6xZs1IdhohIh2Jm5Q0tV1GPiEjEKPGLiESMEr+ISMR0iDJ+EUkfVVVVLF26lC1btjS/ssQlKyuLgQMH0rlz57jWV+IXkTa1dOlSevToQV5eHmb1O0GVlnJ3Vq9ezdKlSxk0aFBc30nbop7S0lLy8vLIyMggLy+P0lL1eyXSHmzZsoU+ffoo6SeImdGnT58W3UGl5RV/aWkphYWFVFZWAlBeXk5hYSEABQUFqQxNREBJP8Fa+vdMyyv+oqIiKisrORW4LlxWWVlJUVFRKsMSEWkX0jLxV1RUAHAicG0Dy0VEEuWVV17htNNOA+CZZ55h4sSJja67bt06fv/737d4HzfddBO33nprq2OsLy0Tf05ODhCMqtGjgeUiIs2prq5u8XdOP/10rr/++kY/b23iT7S0TPzFxcVkZ2ezEegMdAGys7MpLi5OcWQi0h6UlZUxbNgwCgoK2H///Tn77LOprKwkLy+P6667jvz8fB5//HGmT5/O4YcfTn5+Pueccw4bN24E4J///CfDhg0jPz+fadOm1W73gQce4LLLLgNg+fLlnHXWWYwcOZKRI0fyxhtvcP3117NkyRJGjRrFNddcA8BvfvMbDj74YEaMGMHPf/7z2m0VFxczZMgQjjzySD788MOEHn9aVu7GKnAX/+hHsHYt+w8cyDUTJ6piV6S9ufJKmDs3sdscNQruuKPZ1T788EMmTZrEuHHjmDBhQu2VeJ8+fZgzZw6rVq3iG9/4Bi+88ALdunXj5ptv5rbbbuPaa6/l4osv5qWXXmK//fbj29/+doPbv/zyyznmmGN46qmnqK6uZuPGjUycOJH333+fueExT58+ncWLFzNz5kzcndNPP51XX32Vbt26MXXqVObOncv27dvJz89nzJgxCfoDpWnihzD5b9sGEyYw97XXQL17ikgd++yzD+PGjQPgu9/9LnfddRdAbSJ/6623mD9/fu0627Zt4/DDD2fhwoUMGjSIwYMH1363pKRkp+2/9NJLPPjggwBkZmbSs2dP1q5du8M606dPZ/r06YwePRqAjRs3snjxYjZs2MBZZ51FdnY2EBQhJVLaJn4AuncPXsPbMxFpZ+K4Mk+W+k0gY/PdunUDggejTjjhBB555JEd1pubwDsUd+eGG27g+9///g7L70jy3yUty/hrKfGLSCMqKip48803AXj44Yc58sgjd/j8sMMOY8aMGXz00UcAbNq0iUWLFjFs2DDKyspYsmQJwE4nhpjjjz+ee++9FwgqitevX0+PHj3YsGFD7Tonnngi999/f23dwaeffsqKFSs4+uijefrpp9m8eTMbNmzgr3/9a0KPXYlfRCJp6NCh3HPPPey///6sXbuWSy65ZIfP+/btywMPPMB3vvMdRowYUVvMk5WVRUlJCaeeeir5+fnstddeDW7/zjvv5OWXX2b48OGMGTOG+fPn06dPH8aNG8dBBx3ENddcw9e+9jXOPfdcDj/8cIYPH87ZZ5/Nhg0byM/P59vf/jYjR47k5JNP5uCDD07osZu7J3SDyTB27Fhv1UAsc+fC6NEwbRqcdVbC4xKRlluwYAH7779/SmMoKyvjtNNO4/33309pHInU0N/VzGa7+9j66+qKX0QkYpT4RSRy8vLy0upqv6XSO/H3CJ/brVOZIiISdemd+Lt2BTNd8YuI1JHeiT8jA7p1U+IXEakjvRM/BMU9SvwiIrXSP/F3764yfhFJG91jjVZ2QTQSv674RaQd2759e5vuL/0Tv4p6RDq0ZIyfHeuW+YILLmDIkCEUFBTwwgsvMG7cOAYPHszMmTPZtGkTEyZM4JBDDmH06NH85S9/qf3uUUcdRX5+Pvn5+bzxxhsALFu2jKOPPppRo0Zx0EEH8dprrwE7XqE/8cQTXHDBBQBccMEF/OAHP+DQQw/l2muvZcmSJZx00kmMGTOGo446ioULFwLw8ccf1z7Z+9Of/nSXjx0IOglq79OYMWO81U45xX1Xvi8iCTV//vy4150yZYpnZ2c7UDtlZ2f7lClTdimGjz/+2DMzM33evHleXV3t+fn5/r3vfc9ramr86aef9jPOOMNvuOEGf+ihh9zdfe3atT548GDfuHGjb9q0yTdv3uzu7osWLfJYfrr11lv9V7/6lbu7b9++3T///HN3d+/WrVvtfh9//HEfP368u7uPHz/eTz31VN++fbu7ux933HG+aNEid3d/6623/Nhjj3V3969//es+efJkd3e/++67d9heXQ39XYFZ3kBOTe/eOSEo6gk7UxKRjiU2fnZdsfGzd3V8jUGDBjF8+HAADjzwQI4//njMjOHDh1NWVsbSpUt55plnaoc83LJlCxUVFQwYMIDLLruMuXPnkpmZyaJFiwA4+OCDmTBhAlVVVZx55pmMGjWq2RjOOeccMjMz2bhxI2+88QbnnHNO7Wdbt24FYMaMGTz55JMAnHfeeVx33XUNbqslopH4VdQj0iE1Nk52IsbP3m233WrfZ2Rk1M5nZGSwfft2MjMzefLJJxk6dOgO37vpppvo168f7777LjU1NWRlZQFw9NFH8+qrr/K3v/2NCy64gKuvvprzzz9/h+6ft2zZssO2Yl1A19TUsMceezTa5XP9LqR3VdLK+M0sy8xmmtm7ZvaBmf0iXD7IzN42s4/M7FEz65KsGACV8Yt0YI2Nk90W42efeOKJ/O53v8PDjizfeecdANavX0///v3JyMjgoYceqh2bt7y8nH79+nHxxRdz0UUXMWfOHAD69evHggULqKmp4amnnmpwX7vvvjuDBg3i8ccfB4Ii+HfffReAcePGMXXqVICE1G9Acit3twLHuftIYBRwkpkdBtwM3O7u+wFrgQuTGMMXzTk7QC+kIrKj2PjZdbXV+Nk33ngjVVVVjBgxggMPPJAbb7wRgEsvvZTJkyczcuRIFi5cWHvV/sorrzBy5EhGjx7No48+yhVXXAHAxIkTOe200zjiiCPo379/o/srLS1l0qRJjBw5kgMPPLC2MvnOO+/knnvuYfjw4Xz66aeJObiGCv4TPQHZwBzgUGAV0ClcfjjwXHPf36XK3V//2h3cKytbvw0RSZiWVO66BxW8ubm5bmaem5u7yxW76ardVO6aWSYwG9gPuAdYAqxz91ij1aXA3o18txAohF28rYt11LZxY9B3j4h0KAUFBbtckSs7Smo7fnevdvdRwEDgEGBYC75b4u5j3X1s3759Wx9ErA2tnt4VEQHa6AEud18HvExQtLOHmcXuNAYCCSq0aoT65Bdpd1x1bgnV0r9nMlv19DWzPcL3XYETgAUEJ4Czw9XGA39JVgzAjkU9IpJyWVlZrF69Wsk/Qdyd1atX1zYrjUcyy/j7A5PDcv4M4DF3f9bM5gNTzexXwDvApCTGoCt+kXZm4MCBLF26lJUrV6Y6lLSRlZXFwIED414/aYnf3ecBoxtY/h+C8v428bd//YtTgbNPPJFZubkUFxerokgkhTp37sygQYNSHUakpXUnbaWlpVzzy18C0J3gAYvCwsKEPQQhItIRpXXiLyoq4rPwEeme4bJYPx8iIlGV1n31VFRUYEAN0LvechGRqErrK/6cnBxqgPVAr3rLRUSiKq0Tf6yfj7V8ccXfVv18iIi0V2ld1BNrvbNxwgR6bdtGrlr1iIhgHeEhirFjx/qsWbNav4ETToBNmyAcIk1EJArMbLa7j62/PK2Lemr16gVr1qQ6ChGRdiEaib93b1i7NtVRiIi0C9FI/L16BYm/AxRriYgkW3QSf1VVUM4vIhJx0Uj8vcPGnCruERGJSOLvFT6+pcQvIhKxxK+WPSIiEUn8KuoREakVjcSvK34RkVrRSvy64hcRiUji79EDMjOV+EVEiEriN1O3DSIioWgkfvji6V0RkYiLTuJXfz0iIkCUEr+KekREgKglfl3xi4hEKPGrqEdEBIhI4i8tLeWuKVOoWbOGQbm5lJaWpjokEZGUSfvEX1paSmFhIWXr15MBrK2ooLCwUMlfRCIr7RN/UVERlZWVxKp1ewOVlZUUFRWlMiwRkZRJ+8RfUVEBUJv4+9RbLiISNWmf+HNycgBYFs73r7dcRCRq0j7xFxcXk52dzX/D+QFAdnY2xcXFqQxLRCRl0j7xFxQUUFJSQlZODjXA/j17UlJSQkFBQapDExFJiaQlfjPbx8xeNrP5ZvaBmV0RLr/JzD41s7nhdEqyYogpKChgSXk5GV/6ElecfbaSvohEWqckbns78GN3n2NmPYDZZvZ8+Nnt7n5rEvfdsAED4L//bX49EZE0lrTE7+7LCOtU3X2DmS0A9k7W/uIyYAB88klKQxARSbVGE7+ZXR3H9ze5+x+bW8nM8oDRwNvAOOAyMzsfmEVwV7BTXwpmVggUQgJb4AwYADNnJmZbIiIdVFNl/NcA3YEeTUw/bm4HZtYdeBK40t0/B+4F9gVGEdwR/Lah77l7ibuPdfexffv2jfd4mjZgAKxYAVVVidmeiEgH1FRRz0Pu/sumvmxm3Zr5vDNB0i9192kA7r68zuf3Ac/GH+4uGjAgeP3sM9hnnzbbrYhIe9LoFb+7X9vcl5tax8wMmAQscPfb6izvX2e1s4D34ws1AWKJXxW8IhJhzTbnNLMrzGx3C0wyszlm9rU4tj0OOA84rl7TzVvM7D0zmwccC1y1a4fQAkr8IiJxteqZ4O53mtmJQC+CZP4QML2pL7n764A18NHfWxxloijxi4jE9QBXLHmfQlDu/wENJ/T2r29fyMxU4heRSIsn8c82s+kEif+58GGsmuSGlSQZGdC/vxK/iERaPEU9FxI0vfyPu1eaWR/ge0mNKpn09K6IRFyzid/da4A5deZXA6uTGVRSDRgAH32U6ihERFKmVZ20mdmc5tdqp3TFLyIR16rE7+75iQ6kzQwYAGvWwJYtqY5ERCQl0r4//p3EmnQuW9b0eiIiaSqeB7gOM7N/m9lGM9tmZtVm9nlbBJcUassvIhEXzxX/3cB3gMVAV+Ai4J5kBpVUSvwiEnFxFfW4+0dAprtXu/ufgZOSG1YSKfGLSMTF046/0sy6AHPN7BaCrpQ7bt1A797QpYsSv4hEVjwJ/LvhepcBm4B9gG8mM6ikMlOTThGJtKZG4BoM3EowaMp7wE/c/RdtFVhSKfGLSIQ1dcV/P8EgKd8keHL3d20SUZKVlpbyt7lzWfDSS+Tl5VFaWprqkERE2lRTib+Hu9/n7h+6+2+AvDaKKWlKS0spLCzko8pKBgDl5eUUFhYq+YtIpDSV+LPMbLSZ5ZtZPtC13nyHU1RURGVlJf8FehIMKFxZWUlRUVGKIxMRaTtNter5DLitkXkHjktWUMlSUVEBwJJwfj9gbp3lIiJR0Gjid/evtGEcbSInJ4fy8nIWhvPDCBJ/Tk5O6oISEWljjRb1xFOc09GKfIqLi8nOzuYjoJog8WdnZ1NcXJziyERE2k5TRT1/NrOv0PQwi5OA0YkMKJkKCgqAoKz/4/Jy8rOzKSkpqV0uIhIF5u4Nf2BWRjDEYlOJf6W7H5KEuHYwduxYnzVrVmI3etppsHQpzJ2b2O2KiLQTZjbb3cfWX95UGX9eUiNKtaFD4cUXoaYmGItXRCQiopvxhg0LBmNRix4RiZhoJ36AhQubXk9EJM0o8Svxi0jExDMC1zQzO9XM0uskseeeQRfNH36Y6khERNpUPMn898C5wGIzm2hmQ5McU9swC676dcUvIhHTbOJ39xfcvQDIB8qAF8zsDTP7npl1TnaASaXELyIRFFfxjZn1AS4gGG/3HeBOghPB80mLrC0MHQqffQbr1qU6EhGRNhNPGf9TwGtANvB1dz/d3R919x8RdHDZccUqeFXOLyIREs8V/33ufoC7/9rdlwGY2W4ADT0RFmNm+5jZy2Y238w+MLMrwuW9zex5M1scvvZKyJG0hlr2iEgExZP4f9XAsjfj+N524MfufgBwGPBDMzsAuB540d0HAy+G86kxaBB07qzELyKR0tSYu18C9iYcgIUv+uzZnaDYp0nh3cGy8P0GM1sQbu8M4CvhapOBV4DrWhf+LurcGfbbT0U9IhIpTfXOeSJBhe5AdhyQZQPwvy3ZiZnlEfTi+TbQL1ZkRDC4S7+WbCuRSktL6VNWRs6CBZySl0dxcbF66hSRtNdUJ22Tgclm9k13f7K1OzCz7sCTwJXu/rnZF519urubWYPdg5pZIVAIyRkoJTb+7k83b+Y44NNw/F1AyV9E0lpT3TJ/192nmNmPCYZa3IG739bA1+pvozPwLPBcbH0z+xD4irsvM7P+wCvu3uRDYcnoljkvL4/y8nLOJyhvGgosAnJzcykrK0vovkREUqGxbpmbqtztFr52B3o0MDW3QyMYqGVBvZPEM8D48P144C/NRp8EsXF254Xz+fWWi4ikq6aKev4Yvv6ildseB5wHvGdmc8Nl/wtMBB4zswuBcuBbrdz+LomNv/seQaXFOGAqGn9XRNJfPA9w3WJmu5tZZzN70cxWmtl3m/ueu7/u7ubuI9x9VDj93d1Xu/vx7j7Y3b/q7msScygtExt/txp4CzgSjb8rItEQTzv+r7n758BpBH317Adck8yg2kJBQQElJSXk5uYyAxgO/PmOO1SxKyJpL57EHysOOhV43N3XJzGeNlVQUEBZWRk3Pf88mcC3VMwjIhEQT+J/1swWAmOAF82sL7AluWG1sUMPDcbdnTEj1ZGIiCRdPN0yXw8cAYx19ypgE8HTt+mjRw8YNQpefz3VkYiIJF08lbudge8Cj5rZE8CFwOpkB9bmxo2Dt9+GqqpURyIiklTxFPXcS1DM8/twyg+XpZcjj4TKSpg7N9WRiIgkVVN99cQc7O4j68y/ZGbvJiuglBk3LnidMQMOPji1sYiIJFE8V/zVZrZvbMbMvgxUJy+kFNl7b8jLUzm/iKS9eK74rwFeNrP/EHTNnAt8L6lRpcq4cfDii+AeDMYuIpKGmk387v6imQ0m6McM4EN335rcsFJk3DgoLYX//Af23bf59UVEOqB4WvX8EOjq7vPcfR6QbWaXJj+0FDjyyOBV7flFJI3FU8Z/sbuvi824+1rg4qRFlEoHHgg9e6qcX0TSWjyJP9PqjJ5iZplAl+SFlDqljzzCy1u38sF995GXl0dpaWmqQxIRSbh4Ev8/CR7eOt7MjgceCZelldiIXC9s2cKBQGU4IpeSv4ikm0ZH4KpdwSyDYAjEr4aLngf+5O5t1qQzGSNw1RcbkWsE8C5wCfAHNCKXiHRcjY3A1Wzir7eR3sDAsJK3zbRF4s/IyCD2t3gfWAMcDZgZNTU1Sd23iEgytGboxdgXXwkHYukNzAbuM7PbkxFkKtUdeeth4CggB43IJSLpJ54y/p7hQCzfAB5090OB45MbVtuLjcgFQeIHOL9zZ43IJSJpJ54ndzuZWX+CsXGLkhxPysRG3ioqKqK8ooLZXbpw9V570UsjcolImonniv+XwHPAR+7+77CvnsXJDSs1YiNy1dTUMObWW+n1ySfw3nupDktEJKHiGYjl8XDA9EvD+f+4+zeTH1qKnXMOZGbCI4+kOhIRkYRqtKjHzK5191vM7HfATk1/3P3ypEaWav36wVe/Cg8/DMXF6rRNRNJGU2X888PX5LajbM/OPRfGj4c334Qjjkh1NCIiCdFU4j/ZzNa6++Q2i6a9Oess+P73gx47lfhFJE00Vca/CLjVzMrM7BYzG91WQbUbPXrA6afDY49pLF4RSRuNJn53v9PdDweOIRhc/X4zW2hmPzezIW0WYaqdey6sWgUvvJDqSEREEiKeVj3l7n6zu48GvgOcCSxIdmDtxkknsTU7m2nnnENGRoZ67RSRDq/ZB7jMrBNwMvA/BE/svgLclNSo2pHSJ55g29atnFNdTRZQHvbaCV889CUi0pE02kmbmZ1AcIV/CjATmAr8xd03tV14gbbopK0xeXl5DCov52VgPPBguFy9dopIe9fi3jnN7CWCbmueDEfdSplUJv6MjAxwZw6wOzAMqEK9dopI+9fi3jnd/Th3/1Oqk36q5eTk4MD1wJeB79dZLiLSEcXTV0+rmNn9ZrbCzN6vs+wmM/vUzOaG0ynJ2n+ixHrtfA54CbgR2KtrV/XaKSIdVtISP/AAcFIDy29391Hh9Pck7j8hCgoKKCkpITc3lxuAvYDnTzpJFbsi0mElLfG7+6sEA1l1eLFeO992h7PPZsT06bB8earDEhFplWRe8TfmMjObFxYF9WpsJTMrNLNZZjZr5cqVbRlf04qLYcsW+NWvUh2JiEirtHXivxfYFxgFLAN+29iK7l7i7mPdfWzfvn3bKLw4DBkCF10Ef/gDLFmS6mhERFqsTRO/uy9392p3rwHuAw5py/0nzM9+Bp07w403pjoSEZEWa9PEHw7hGHMW8H5j67ZrAwbAVVcFg7TMmZPqaEREWiSZzTkfAd4EhprZUjO7ELjFzN4zs3nAscBVydp/0l17LVu7d+fVcePUh4+IdCjxDLbeKu7+nQYWT0rW/tpa6bPP8t7WrUysquJY4CX14SMiHUSjXTa0J6nssqExeXl5fFZezkKgmqCyYg3qw0dE2o8Wd9kgTauoqGArQZelA4HHCG6fKioqUhqXiEhzlPhbKdZXz9tAIUF/1b9FffiISPunxN9KsT58IOiq+TbgcmDqCSekMiwRkWYp8bdS3T58zIx7cnJYdtBBHDZ5MsyYkerwREQapcrdRFq7Fg45BD7/HGbNgn32SXVEIhJhqtxtC716wTPPwObNcOaZUFmZ6ohERHaixJ9o++8fPNH7zjtw4YXQAe6oRCRalPiT4dRT4de/hqlTubl3bz3ZKyLtStKe3I260r33plNmJtesW8frwLN6sldE2glV7iZJXl4ey8vLeQ04CLiEYEgyPdkrIm1FlbttrKKigi3AicDrwJ+Bu4D/lpenNC4RESX+JIk9wbuGYODh3wI/Al7bbTdYsSKFkYlI1CnxJ0ndJ3urgZ8AE7p0YUxNDYwdC7NnpzQ+EYkuJf4kqf9kb25uLsfffz+d3norWOHII+Ghh1IbpIhEkip3U2HFCvjWt+Bf/4Irr4Tf/AY6qYGViCSWKnfbk732guefh8svhzvu4LORIxm9zz5q7y8ibUKXmanSuTPceSdvbtnC6JISngK+Abyj9v4ikmS64k+x7zz3HEcRnIH/TTA2ZZ/KSoqKilIbmIikLSX+FKuoqGAWMBK4AzgXWAxcWV4OK1emMjQRSVNK/ClWt73/T4DBwEMEbf758pfhppuCbp5FRBJEiT/F6rb3B1gKXJGdzd9vuQVOOgl+8YvgBHD77bBlS+oCFZG0ocSfYg219y8pKeHr11wDjz8O//435OfD1VfDkCEwaRJs305paSl5eXlqCSQiLaZ2/B3FSy/BDTfAzJms79+fS1ev5uFt22o/zs7OpqSkRC2BRKSW2vF3dMcdB2+9BU89xfJVqyjdto2ZwDeBLKBSLYFEJE5K/B2JGZx5JgdUVTEe6As8ASwn6PL5gPJyqKpKZYQi0gEo8XdAA3NzeRDYD/gq8DhwBvB3gAED4JJL4NVXoaYmhVGKSHulxN8BxVoCVQMvAhcBg7p25V9XXQXHHw+TJ8Mxx0BuLvzkJzBr1k5j/6pyWCTC3L3dT2PGjHHZ0ZQpUzw3N9fNzHNzc33KlClffLhhg/vDD7t//evunTu7g/vgwe433ug+f75PmTLFs7OzHaidsrOzd9yGiHR4wCxvIKeqVU+6W7MGpk2DRx6Bl18Gd+Z37syDVVVMBeqOB6ZhIUXSi1r1RFXv3nDRRfDii/Dpp3DHHXxeVcVEoAyYB/wBOB/oXF6+U5GQiKSfpCV+M7vfzFaY2ft1lvU2s+fNbHH42itZ+5cG9O8PV1zB/+TmMgi4juBJ4W8Dkwn6CKJfPzjzzGCMgBkz9LSwSBpK5hX/AwTDzdZ1PfCiuw8mqJe8Pon7l0YUFxezIjubW4BTgN7AwVlZvD1hApxyCnzwAVx7bTBKWM+ecMQRQSXxtGnw2Wcpjl5EdlXSEr+7v0rQ91hdZxBcXBK+npms/Uvj6ncTkZOby5V/+hOHTpoEDzwAixfD8uXw9NNwxRWQkQF33w3f/GZw17DvvnDeeXDvvTBvHlRXq5WQSAeS1MpdM8sDnnX3g8L5de6+R/jegLWx+Qa+WwgUAuTk5IwpLy9vaDVpK1u3wpw58MYbwTRjRnByAKqysnh12zZeralhJvAhsLJrV/5w333qQkIkhRqr3E1Z4g/n17p7s+X8atXTDrnDxx/DG2/w0CWXMGLjRobzxS3kVqC8c2eGnHpq0Llc3WmvvYKnkEUkqRpL/G099OJyM+vv7svMrD+woo33L4liFnQX/eUvM/7883GgBzCKYEyBIcCQqiqGfPgh/O1vO3Yl0bPnzieDIUNg8GDo0SMVRyMSKW2d+J8BxgMTw9e/tPH+JQlycnIoLy9nA/BaOEHwXMBZ8+dDdTWUl8OiRTtOr78ODz+8YxPS/v1h6NCdTwqDBkGXLpSWllJUVERFRQU5OTkUFxerOEmkhZKW+M3sEeArwJ5mthT4OUHCf8zMLiR4duhbydq/tJ3i4mIKCwuprKysXZadnU1xcXEwk5lZe3fASfUaem3eDB99tPNJYdo0WLXqi/UyM/l8zz3ps3Ill9fUUA4sKy/ngQsvpPtnn3HGD34A3bol/2BF0oCe3JWESMqV+Jo1QQuj8GTw7G23MbCykiFAdkPr9+gR3DE0N+2xh+oYJBJSUrmbKEr8ApCRkUHs32tvoH84DQAmT5wIy5btPG3atPOGsrLgS19q+KQwYMAX7/fcM2jKWoeKmqQjaS+VuyKtFqtLgOABkTXABwR1CVx3XcNf2rCh4RPCf/8bvC5YEIxutm7dzt/t1Cl4kjk8ESzeuJGy117jxO3bWQ2sKy/n3gsvpPvy5ZwxfnxwJ5GZmZRjF0kkXfFLh1FaWtpgXUJChpzcvDl4Krmhk0Q4rXzvPfrU1DT91GOPHtCrV3AS2GOPL97Xf21oWbduu1wEpTsSqUtX/NLhxRJYUhJb165By6FBgxpdpV9GBhlAP6BXOO0Rvj54553BXcPatTu+fvzxF/MbNjQdQ6dOrTthhK+ljz22w4mxvLycwsJCACV/2YGu+EXilJeXR0NPkMfdnfX27bB+fXAiaOgk0dhr7P22bU1uvtKMNe6sAz4HNoUT2dmcee65wR1F3Sk7O75lKSq+0t3LrtMVv8guarbZanM6dYI+fYKppdyDnlKbOEn84Wc/oyfBHcjuQHdgL6BbZSU8+2xQ0b1pU8uH5OzSpXUnjXjW6dp1pwp02LlYT3cviaXELxKnpBY1NccsSJJduwYtjxpw16RJzd+RuAd3DrGTwKZNUFnZ9Hxj66xdC0uX7vh5ZWXLx3To2nWnE8S+8+bx2NatbAI2E3QBsrWykk2XXgrz58NuuwVTVtYX71syxb7XqVO7bdqbzDseFfWIpImkVn7Hyz2oKI/35NHI/L/+8Q+6ETyv0RXYDcgKX7tlZgZPgyeCWetPGrtywmls6tIFzBL2W6odv0gEpEu5eLP1KdXVQY+xTU1btjS/TqLW3749cQffpQsbqqrY4s5W4DzglfrHHyeV8YtEQEFBQYdM9PXF1Q1IdnYwtQc1NQk9uUy6+Wa6ENzlrKyzm4qKioSEq8QvIu1OSutTWiMj44s6mAS4Y+rUBu94cnJyErJ9DbYuIu1SQUEBZWVl1NTUUFZW1n6TfhIUFxeTXe9upkUtyJqhxC8i0s7UHx41Nzc3oZX0qtwVEUlTjVXu6opfRCRilPhFRCJGiV9EJGKU+EVEIkaJX0QkYjpEqx4zW0kwOHtT9gRWNbNOOtJxR4uOO3p25dhz3b1v/YUdIvHHw8xmNdRsKd3puKNFxx09yTh2FfWIiESMEr+ISMSkU+IvSXUAKaLjjhYdd/Qk/NjTpoxfRETik05X/CIiEgclfhGRiOnwid/MTjKzD83sIzO7PtXxJJOZ7WNmL5vZfDP7wMyuCJf3NrPnzWxx+Nor1bEmmpllmtk7ZvZsOD/IzN4Of/dHzaxLqmNMBjPbw8yeMLOFZrbAzA6PyO99Vfhv/H0ze8TMstLxNzez+81shZm9X2dZg7+vBe4Kj3+emeW3dr8dOvGbWSZwD3AycADwHTM7ILVRJdV24MfufgBwGPDD8HivB15098HAi+F8urkCWFBn/mbgdnffD1gLXJiSqJLvTuCf7j4MGEnwN0jr39vM9gYuB8a6+0FAJvA/pOdv/gBwUr1ljf2+JwODw6kQuLe1O+3QiR84BPjI3f/j7tuAqcAZKY4padx9mbvPCd9vIEgCexMc8+RwtcnAmSkJMEnMbCBwKvCncN6A44AnwlXS7pgBzKwncDQwCcDdt7n7OtL89w51ArqaWScgG1hGGv7m7v4qsKbe4sZ+3zOABz3wFrCHmfVvzX47euLfG/ikzvzScFnaM7M8YDTwNtDP3ZeFH30G9EtVXElyB3AtUBPO9wHWufv2cD5df/dBBGNt/zks5vqTmXUjzX9vd/8UuBWoIEj464HZROM3h8Z/34Tlu46e+CPJzLoDTwJXuvvndT/zoH1u2rTRNbPTgBXuPjvVsaRAJyAfuNfdRwObqFesk26/N0BYpn0GwYlvANCNnYtDIiFZv29HT/yfAvvUmR8YLktbZtaZIOmXuvu0cPHy2C1f+LoiVfElwTjgdDMrIyjKO46g3HuPsBgA0vd3Xwosdfe3w/knCE4E6fx7A3wV+NjdV7p7FTCN4N9BFH5zaPz3TVi+6+iJ/9/A4LC2vwtBBdAzKY4pacKy7UnAAne/rc5HzwDjw/fjgb+0dWzJ4u43uPtAd88j+H1fcvcC4GXg7HC1tDrmGHf/DPjEzIaGi44H5pPGv3eoAjjMzLLDf/Ox40773zzU2O/7DHB+2LrnMGB9nSKhlnH3Dj0BpwCLgCVAUarjSfKxHklw2zcPmBtOpxCUeb8ILAZeAHqnOtYkHf9XgGfD918GZgIfAY8Du6U6viQd8yhgVvibPw30isLvDfwCWAi8DzwE7JaOvznwCEE9RhXBHd6Fjf2+gBG0YlwCvEfQ6qlV+1WXDSIiEdPRi3pERKSFlPhFRCJGiV9EJGKU+EVEIkaJX0QkYpT4RUQiRolf2j0zKwq76J1nZnPN7NBweZmZ7dkG+y8zs/fC/f/LzHKbWT/PzM5t4T7yzGyzmc2ts2ynLnvrfHaYmd3XyLa6hn+nbW3x95GOR4lf2jUzOxw4Dch39xEEj/N/0vS3dml/nRr56Nhw/68AP21mM3lAixJ/aIm7j6oz/wCN91FzMvDPhj5w983hdv7bihgkApT4pb3rD6xy960A7r7K3esmtB+Z2ZzwinwYgJl1C6+WZ4a9Wp4RLs8zs9fC9eeY2RHh8q+Ey58h6BqgKW8S9ojY2PaAicBR4VX3VRYMIvMbM/t3eNfw/XgO3BvusjfmeOAFMzswPM654bYHx7NtiTYlfmnvpgP7mNkiM/u9mR1T7/NV7p5PMCjFT8JlRQR9+hwCHAv8JuzOeAVwQrj+t4G76mwnH7jC3Yc0E89JBF0n0MT2rgdec/dR7n47wWP46939YOBg4GIzG9SCv8EOwuKbKndfD/wAuDO8wh9L8Ni/SJMau60VaRfcfaOZjQGOIkjij5rZ9e7+QLhKrIfS2cA3wvdfI+jRM3YiyAJyCIo+7jazUUA1UDfJz3T3j5sI5WUz6w1sBG4Ml3VuYnt1fQ0YYWaxDsZ6Eoyi1NT+mvI1ghMiBHcgReFgNdPcfXErtykRosQv7Z67VxOUrb9iZu8R9Fj4QPjx1vC1mi/+PRvwTXf/sO52zOwmYDnBEIYZwJY6H29qJoxjgXVAKUEHYlcDVzWxvR12DfzI3Z9rZh/xOhm4DcDdHzaztwlGKPu7mX3f3V9K0H4kTamoR9o1Mxtar9x6FFDezNeeIyj7t3Abo8PlPYFl7l4DnEcwlmvcPBj96UqCrnF7N7G9DUCPevFcEo6lgJkNCYueWiw8phEEPbNiZl8G/uPudxF03zuiNduVaFHil/auOzDZzOab2TzgAOCmZr7z/wiKYeaZ2QfhPMDvgfFm9i4wjOav8nfiQf/njwA/bGJ784BqM3vXzK4iGCt4PjAnbJr5R+K42zazRwiKcoaa2VIzuxAYA7zjX3Sr+y3g/bAZ6EHAgy09Jokedcss0g5YMIbys+5+UDPr/RT4yN2nxrHNMoI+21clJEhJG0r8Iu2Ame0DvAGsrteWvzXb6kpwp9AXGO7ujTUJlYhS4hcRiRiV8YuIRIwSv4hIxCjxi4hEjBK/iEjE/H/Sr+J9rP1F/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('')\n",
    "plt.scatter(df_ares['Shear rate'],df_ares['Viscosity'], color = 'k', label = 'measured')\n",
    "plt.plot(x1, eta_list[num], color = 'r', label = 'predicted')\n",
    "\n",
    "plt.title(models_name_list[num])\n",
    "#plt.xscale('log')\n",
    "#plt.yscale('log')\n",
    "plt.legend()\n",
    "\n",
    "plt.ylabel('Viscosity [Pa.s]')\n",
    "plt.xlabel('Shear Rate [1/s]')\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9998c7c7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bac29cee",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}