3FLEX-FIT.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#############\n",
    "# \n",
    "# This Software was written by Brandon Bukowski at Northwestern University\n",
    "# If using this code, please cite the originial mathematical models from Brandani and Brandani et al.\n",
    "#\n",
    "# Brandani, S. Analysis of the Piezometric Method for the Study of Diffusion in Microporous Solids: Isothermal Case. Adsorption 1998, 4, 17–24.\n",
    "# Brandani, S.; Brandani, F.; Mangano, E.; Pullumbi, P. Using a Volumetric Apparatus to Identify and Measure the Mass Transfer Resistance in Commercial Adsorbents. Microporous and Mesoporous Materials 2020, 304, 109277.\n",
    "\n",
    "#\n",
    "# As well as our implimentation (currently in preparation)\n",
    "#\n",
    "# Son, F. A.; Bukowski, B. C.; Islamoglu, T.; Snurr, R. Q.; Farha, O. K. Rapid Quantification of Mass Transfer  Barriers in Metal-Organic Framework Crystals. in preparation.\n",
    "#\n",
    "#\n",
    "#########\n",
    "\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "from scipy.optimize import root_scalar\n",
    "import matplotlib.cm as cm\n",
    "from scipy.optimize import minimize_scalar\n",
    "import matplotlib.pylab as pylab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# User inputs are only required for the following 2 boxes. The remained can be run without modification. \n",
    "# if modifications are necessary, descriptions are provided to help troubleshooting\n",
    "\n",
    "# Common errors\n",
    "# RuntimeWarning: invalid value encountered in double_scalars app.launch_new_instance() \n",
    "#       - this can be safely ignored\n",
    "\n",
    "# UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data. fig.canvas.print_figure(bytes_io, **kw)\n",
    "#       - this can be safely ignored \n",
    "\n",
    "#UnboundLocalError: local variable 'j' referenced before assignment\n",
    "#       - some of the chambers in the csv are empty. scroll to the bottom of box 7 and choose the right size of \n",
    "#         valve_states to fix it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data Required\n",
    "\n",
    "###\n",
    "# name of the file you are reading, assumed to be a csv output file from 3Flex equipment\n",
    "###\n",
    "\n",
    "####\n",
    "# info to modify\n",
    "data_suffix = ' '  # Name of your file minus '.csv' goes here!\n",
    "####\n",
    "data = np.genfromtxt(f'{data_suffix}.csv', delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data Required\n",
    "\n",
    "\n",
    "# Chamber_select - choose which sample chamber to study, 0,1,2 corresponds to the first, second, third chamber\n",
    "# L_scale - particle size in meters, assuming 3 sample ports\n",
    "L_scale = [4e-6, 4e-6, 4e-6]\n",
    "# divide by 2 to make radius\n",
    "L_scale = [L/2 for L in L_scale]\n",
    "\n",
    "# max pressure - upper limit of pressure to consider - assuming mmHg\n",
    "max_pressure = [2., 2., 2.]\n",
    "\n",
    "# V_D - Manifold volume - cm^3\n",
    "V_D = 25.531\n",
    "\n",
    "# choose which sample chamber to solve for - remember Python starts from 0. \n",
    "# if 3 sample ports the choices are port 0, port 1, port 2\n",
    "Chamber_select = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "#\n",
    "# automatically reading the different chambers from the equipment, should not need modifying\n",
    "#\n",
    "\n",
    "man_pressure = data[:,0]\n",
    "time = data[:,3]\n",
    "time = time - time[0]\n",
    "time /= 1000\n",
    "U_pressures = [data[:,4], data[:,22], data[:,40]]\n",
    "\n",
    "valves = [data[:,14], data[:,32], data[:,50]]\n",
    "\n",
    "Q_sets = [data[:,11], data[:,29], data[:,47]]\n",
    "\n",
    "# V_Us = [data[0,10], data[0,28], data[0,46]]\n",
    "V_Us = [data[0,5], data[0,23], data[0, 41]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# simple plot showing the pressure profile in the manifold over the full run\n",
    "#\n",
    "plt.plot(time, U_pressures[0])\n",
    "plt.plot(time, U_pressures[1])\n",
    "plt.plot(time, U_pressures[2])\n",
    "plt.title(f'{data_suffix} : all chambers')\n",
    "plt.ylabel('Uptake vessel pressure (Torr)')\n",
    "plt.xlabel('time (s)')\n",
    "plt.legend(['Chamber 1', 'Chamber 2', 'Chamber 3'])\n",
    "#plt.ylim(0,400)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "# helper function to determine when the valve are open and closed\n",
    "##\n",
    "\n",
    "def valve_states(valve):\n",
    "    valve_indices = list(np.where(valve != 0))[0]\n",
    "    index_jump_list = []\n",
    "    valve_state = []\n",
    "    for i, index in enumerate(valve_indices[1:]):\n",
    "        if valve_indices[i] - valve_indices[i-1] != 1:\n",
    "            index_jump_list.append(i)\n",
    "    for j in range(len(index_jump_list) - 1):\n",
    "        valve_state.append(valve_indices[index_jump_list[j]:index_jump_list[j+1]])\n",
    "    valve_state.append(valve_indices[index_jump_list[j+1]:-1])\n",
    "    return valve_state\n",
    "\n",
    "valve_states = [valve_states(valves[0]), valve_states(valves[1]), valve_states(valves[2])]\n",
    "#valve_states = [valve_states(valves[0]), valve_states(valves[1])]\n",
    "#valve_states = [valve_states(valves[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "# function to determine the root-mean-square error between the experimental function (ex) and simulated function (h)\n",
    "## \n",
    "\n",
    "def RMS_calc(h, ex):\n",
    "    return np.sqrt(np.sum((h[2:] - ex[2:])**2) / h[2:].size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "####\n",
    "# Split the data into individual runs based on whether the valves are open or closed. \n",
    "# No user input necessary\n",
    "####\n",
    "\n",
    "\n",
    "P_U_0_overall = []\n",
    "P_D_0_overall = []\n",
    "P_inf_overall = []\n",
    "P_U_0_set = []\n",
    "P_D_0_set = []\n",
    "P_inf_set = []\n",
    "P_close = []\n",
    "scale_time_zip = []\n",
    "time_zip = []\n",
    "uptake_zip = []\n",
    "P_U_overall = []\n",
    "P_D_overall = []\n",
    "V_S_set = []\n",
    "Q_set = []\n",
    "P_Ds = []\n",
    "P_Us = []\n",
    "CJ = Chamber_select\n",
    "for j, vlv_group in enumerate(valve_states[CJ][:]):\n",
    "    dml = (Q_sets[CJ][vlv_group]-Q_sets[CJ][vlv_group[0]-1]) / (Q_sets[CJ][vlv_group[-1]] - Q_sets[CJ][vlv_group[0]-1])\n",
    "    if dml[0] < 0.1 and dml[0] > -0.1 and np.all(dml < 1.05) and np.all(dml > -0.05):\n",
    "        if U_pressures[CJ][vlv_group[0]-1] < max_pressure[CJ]:\n",
    "            P_D = man_pressure[vlv_group]\n",
    "            P_inf = U_pressures[CJ][vlv_group][-1]\n",
    "            if P_D[0]-P_inf > 0:\n",
    "                P_Ds.append(man_pressure[vlv_group])\n",
    "                P_Us.append(U_pressures[CJ][vlv_group])\n",
    "                P_D_0_set.append(man_pressure[vlv_group[0]-1])\n",
    "                P_U_0_set.append(U_pressures[CJ][vlv_group[0]-1])\n",
    "                P_inf_set.append(man_pressure[vlv_group][-1])\n",
    "                time_zip.append(time[vlv_group] - time[vlv_group][0])\n",
    "                scale_time_zip.append(time[vlv_group] - time[vlv_group][0])         \n",
    "                dml = (Q_sets[CJ][vlv_group]-Q_sets[CJ][vlv_group[0]-1]) / (Q_sets[CJ][vlv_group[-1]] - Q_sets[CJ][vlv_group[0]-1])\n",
    "                Q_set.append(Q_sets[CJ][vlv_group]- Q_sets[CJ][vlv_group[0]-1])\n",
    "                V_S_set.append(Q_sets[CJ][vlv_group[-1]]- Q_sets[CJ][vlv_group[0]-1])\n",
    "                uptake_zip.append(dml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "# plot the pressure in the uptake cell when the valves are opened\n",
    "##\n",
    "\n",
    "plt.figure()\n",
    "plt.figure(figsize=(12,12))\n",
    "import matplotlib.cm as cm\n",
    "colors = cm.rainbow(np.linspace(0,1,len(uptake_zip)))\n",
    "def sequential_color_plot(intensity_data, position_data, colors=colors, data_suffix=data_suffix):\n",
    "    fig, ax = plt.subplots(figsize=(8,8))\n",
    "    for i, c in zip(range(len(intensity_data)), colors):\n",
    "        ax.plot(np.sqrt(position_data[i]), intensity_data[i], color=c)\n",
    "        ax.scatter(np.sqrt(position_data[i]), intensity_data[i], color=c)\n",
    "    #ax.set_ylim(0,1)\n",
    "    ax.set_title(f'{data_suffix} : Chamber {Chamber_select+1}')\n",
    "    ax.set_xlim(0,5)\n",
    "    ax.set_xlabel('$t^{1/2}$ [$s^{1/2}$]')\n",
    "    ax.set_ylabel('fractional uptake')\n",
    "sequential_color_plot(uptake_zip[:], scale_time_zip[:])\n",
    "def sequential_color_plot(intensity_data, position_data, colors=colors, data_suffix=data_suffix):\n",
    "    fig, ax = plt.subplots(figsize=(8,8))\n",
    "    for i, c in zip(range(len(intensity_data)), colors):\n",
    "        ax.plot(position_data[i], intensity_data[i], color=c)\n",
    "        ax.scatter(position_data[i], intensity_data[i], color=c)\n",
    "    ax.set_title(f'{data_suffix} : Chamber {Chamber_select+1}')\n",
    "    ax.set_ylim(0,1)\n",
    "    ax.set_xlim(0,25)\n",
    "    ax.set_xlabel('$t$ [$s$]')\n",
    "    ax.set_ylabel('fractional uptake')\n",
    "sequential_color_plot(uptake_zip[:], time_zip[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "####\n",
    "# combined surface + intracrysalline objective function to minimze\n",
    "####\n",
    "\n",
    "def objective_function_comb(x0, uptake_bundle):\n",
    "    # requires - x0, initial guess for D_R and a_LDF : uptake bundle, parameters for solver\n",
    "    # returns - RMS, error between analytical solution and experimental data\n",
    "    \n",
    "    # D_R = dimensionless intracrystalline diffusion coefficient\n",
    "    D_R = x0[0] \n",
    "    # a_LDF = ratio between surface diffusion rate and intracrystalline diffusion rate - alpha_LDF\n",
    "    a_LDF = x0[1]\n",
    " \n",
    "    # pulling parameters from the bundle\n",
    "    # delta_gamma = delta and gamma determined from the fraction of moles in manifold / uptake cell\n",
    "    # pressure_bundle = set of initial and pressures- see below\n",
    "    # tau_times = experimental timesteps (s)\n",
    "    # scaled_time = timesteps for simulated values\n",
    "    # sig_D = sigma_D - experimental reduced pressure in the dosing volume\n",
    "    # omega_c = intrinsic valve constant that depends on the instument\n",
    "    # M = how many terms of the solution series we use\n",
    "    delta_gamma, pressure_bundle, tau_times, scaled_time, sig_D, omega_c, M = uptake_bundle\n",
    "    \n",
    "    # omega = the true valve constant\n",
    "    omega = omega_c / D_R\n",
    "\n",
    "    RMS = 0\n",
    "    # analytical solution \n",
    "    analytical = combined_diffusion(omega, D_R, delta_gamma, a_LDF, pressure_bundle, scaled_time, M)\n",
    "\n",
    "    # calculate RMS between analytical and experimental values\n",
    "    RMS = RMS_calc(analytical[:], sig_D[:])\n",
    "    \n",
    "    return(RMS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#######\n",
    "# The combined diffusion function from Brandani\n",
    "#######\n",
    "\n",
    "def combined_diffusion(omega, D_R, delta_gamma, a_LDF, pressure_bundle, times, M):\n",
    "    delta, gamma = delta_gamma\n",
    "    P_D_0, P_U_0, P_inf = pressure_bundle\n",
    "      \n",
    "    Beta_set = []\n",
    "    z_set = []\n",
    "    a_set = []\n",
    "\n",
    "    def f(x, valve_bundle):\n",
    "        gamma, delta, omega, a_LDF = valve_bundle\n",
    "        z = 1 + ((a_LDF*(gamma*x**2*(omega-x**2) + omega*delta*x**2)) / ((omega-x**2)*(a_LDF-gamma*x**2)-omega*delta*x**2))\n",
    "        return (x/np.tan(x))-z\n",
    "  \n",
    "    valve_bundle = [gamma, delta, omega, a_LDF] \n",
    "\n",
    "    a_q = gamma #a*c + b + d\n",
    "    b_q = -omega*gamma-a_LDF-omega*delta # a*c - b - d\n",
    "    c_q = omega*a_LDF    \n",
    "    \n",
    "    LB_wall = np.pi\n",
    "    window = 1\n",
    "    while len(Beta_set) < M:\n",
    "        LB = window * LB_wall + 1e-8\n",
    "        UB = (window + 1) * LB_wall - 1e-8\n",
    "        root_out = scipy.optimize.root_scalar(f, valve_bundle, bracket=[LB, UB], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "        Beta_set.append(root_out.root)\n",
    "        window += 1\n",
    "        \n",
    "    root1_old = np.sqrt((-b_q-np.sqrt((b_q**2)-4*a_q*c_q))/(2*a_q))\n",
    "    left_window = np.floor(root1_old / np.pi)\n",
    "    #print(f\"unrefined root: {root1_old}\")\n",
    "    trivial_root = scipy.optimize.root_scalar(f, valve_bundle, bracket=[left_window*np.pi + 1e-8, (left_window+1)*np.pi - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "    base_root = trivial_root.root\n",
    "    #print(f\"trivial root: {base_root}\")\n",
    "    root_error = np.inf\n",
    "    start = left_window*np.pi + 1e-8\n",
    "    #print(f\"guess root {root1_old}\")\n",
    "\n",
    "    if base_root < root1_old:\n",
    "        try: \n",
    "            br_1 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[left_window*np.pi+ 1e-8, base_root - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_1.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        try: \n",
    "            br_2 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[base_root + 1e-8, root1_old - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_2.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        try: \n",
    "            br_3 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[root1_old + 1e-8, (left_window+1)*np.pi - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_3.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "# -----------------------------------------------\n",
    "\n",
    "    elif root1_old < base_root:\n",
    "        try: \n",
    "            br_1 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[left_window*np.pi+ 1e-12, root1_old - 1e-12], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_1.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        try: \n",
    "            br_2 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[root1_old + 1e-12, base_root - 1e-12], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_2.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        try: \n",
    "            br_3 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[base_root + 1e-12, (left_window+1)*np.pi - 1e-12], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_3.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "                  \n",
    "    Beta_set.append(good_root)\n",
    "    \n",
    "    #print(\" ------------------- \")\n",
    "\n",
    "    root2_old = np.sqrt((-b_q+np.sqrt((b_q**2)-4*a_q*c_q))/(2*a_q))\n",
    "    left_window = np.floor(root2_old / np.pi)\n",
    "    trivial_root = scipy.optimize.root_scalar(f, valve_bundle, bracket=[left_window*np.pi + 1e-12, (left_window+1)*np.pi - 1e-12], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "    base_root = trivial_root.root\n",
    "    root_error = np.inf\n",
    "    start = left_window*np.pi + 1e-8\n",
    "\n",
    "    if base_root < root2_old:\n",
    "        try: \n",
    "            br_1 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[left_window*np.pi+ 1e-8, base_root - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_1.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        try: \n",
    "            br_2 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[base_root + 1e-8, root2_old - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_2.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        try: \n",
    "            br_3 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[root2_old + 1e-8, (left_window+1)*np.pi - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_3.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "# -----------------------------------------------\n",
    "\n",
    "    elif root2_old < base_root:\n",
    "        try: \n",
    "            br_1 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[left_window*np.pi+ 1e-8, root2_old - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_1.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        try: \n",
    "            br_2 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[root2_old + 1e-8, base_root - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_2.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "\n",
    "        try: \n",
    "            br_3 = scipy.optimize.root_scalar(f, valve_bundle, bracket=[base_root + 1e-8, (left_window+1)*np.pi - 1e-8], method='brentq', xtol=1e-12, rtol=1e-12)\n",
    "            good_root = br_3.root\n",
    "        except ValueError:\n",
    "            pass\n",
    "                \n",
    "                    \n",
    "    Beta_set.append(good_root)    \n",
    "    Beta_set = np.array(Beta_set)\n",
    " \n",
    "    def z_fun(Beta, gamma, delta, omega, a_LDF):\n",
    "        z = 1 + ((a_LDF*(gamma*Beta**2*(omega-Beta**2) + omega*delta*Beta**2)) / ((omega-Beta**2)*(a_LDF-gamma*Beta**2)-omega*delta*Beta**2))            \n",
    "        return z\n",
    "    z_set = z_fun(Beta_set, gamma, delta, omega, a_LDF)\n",
    "\n",
    "    def a_fun(Beta, gamma, delta, omega, z, a_LDF):\n",
    "        b2 = Beta**2\n",
    "        gb2 = gamma*b2\n",
    "        od = omega*delta\n",
    "        odb2 = od*b2\n",
    "        T1 = (z-1)*(od-gb2+a_LDF) + a_LDF*(od-gb2)\n",
    "        T2 = ((omega-b2)*(a_LDF-gb2)-odb2) / (2*b2)\n",
    "        T3 = (z**2 + b2 - z)\n",
    "        T4 = (z - 1) * (a_LDF+gamma*(omega-b2) - gb2 + od)\n",
    "        T5 =  a_LDF*(gamma*(omega-b2)-gb2+od)\n",
    "\n",
    "        return T1/(T2*T3+T4+T5)\n",
    "    a_set = a_fun(Beta_set, gamma, delta, omega, z_set, a_LDF)\n",
    "\n",
    "    #print(f\"total a_set sum {np.sum(a_set)}\")\n",
    "    \n",
    "    pre_factor = (3*delta) / (1 + 3*delta + 3*gamma)\n",
    "    y_intra = []\n",
    "    for new_time in times:\n",
    "        cumulant = pre_factor\n",
    "        for i in range(len(Beta_set)):\n",
    "            #print(a_set[i]+pre_factor)\n",
    "            cumulant += a_set[i] * np.exp(-(Beta_set[i]**2)*D_R*new_time)\n",
    "        y_intra.append(cumulant)\n",
    "    #print(y_intra[0])\n",
    "    rho_D_0 =(P_D_0 - P_U_0) / (P_inf - P_U_0) # 3.61 / 1.4\n",
    "    rho_D_intra = np.asarray(y_intra) * rho_D_0\n",
    "    sig_D_intra = (rho_D_intra *(P_inf-P_U_0)+P_U_0-P_inf) / (P_D_0-P_inf)\n",
    "    \n",
    "    return sig_D_intra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D_R_set = []\n",
    "aLDF_set = []\n",
    "omega_set = []\n",
    "omega_LDF = []\n",
    "rms_set = []\n",
    "rms2_set = []\n",
    "P_U_0_check = []\n",
    "P_inf_check = []\n",
    "P_D_0_check = []\n",
    "expts_data = []\n",
    "timeset_data = []\n",
    "fitted_data = []\n",
    "D_R_overall = []\n",
    "aLDF_overall = []\n",
    "rms_overall = []\n",
    "\n",
    "for exp in range(len(P_U_0_set[:])):\n",
    "    #it = 0\n",
    "    plt.figure(figsize=(5,5))\n",
    "\n",
    "\n",
    "    P_D_0 = P_D_0_set[exp] #[index_jump_list[0]]-1]\n",
    "    P_U_0 = P_U_0_set[exp] #[index_jump_list[0]]-1]\n",
    "    P_inf = P_inf_set[exp]\n",
    "    P_D = P_Ds[exp]\n",
    "    P_U = P_Us[exp]\n",
    "    V_U = V_Us[CJ]\n",
    "    #print(V_U)\n",
    "    V_S = V_S_set[exp]\n",
    "    time_cut = time_zip[exp]\n",
    "    sig_D = (P_D - P_inf) / (P_D_0 - P_inf)\n",
    "    #print(time_cut)\n",
    "    plt.scatter(time_cut, sig_D, c='k')\n",
    "    plt.plot(time_cut, sig_D, c='k')\n",
    "    expts_data.append(sig_D)\n",
    "    timeset_data.append(time_cut)\n",
    "    plt.yscale('log')\n",
    "    plt.ylabel('fractional uptake')\n",
    "    plt.xlabel('time (s)')\n",
    "    plt.ylim(0.01, 1)\n",
    "    plt.xlim(0,25)\n",
    "\n",
    "    #######################################################\n",
    "    # These are some initial conditions for the solver, can be modified if the fits are poor.\n",
    "    # M is the number of terms in the transcendental tan function. Can be increased if needed.\n",
    "\n",
    "    omega_c = 10\n",
    "    M = 100\n",
    "    D_R_list = [0.1]\n",
    "    a_LDF_list = [0.001]\n",
    "    #\n",
    "    #######################################################    \n",
    "\n",
    "\n",
    "    tau_times = np.linspace(0,25,100)\n",
    "    kai = (P_inf - P_U_0) / (P_D_0 - P_U_0)\n",
    "    eps = V_U/V_D\n",
    "    delt_temp = -kai / (3*kai+3*eps*kai-3)\n",
    "\n",
    "    gam_temp = delt_temp * eps\n",
    "    delta_gamma = [delt_temp, gam_temp]\n",
    "    new_preex = 3*delt_temp / (1+3*delt_temp+3*gam_temp)\n",
    "\n",
    "    pressure_bundle = [P_D_0, P_U_0, P_inf]\n",
    "    \n",
    "    ratio_set = np.array(D_R_list) * np.array(a_LDF_list)\n",
    "\n",
    "    uptake_bundle = [delta_gamma, pressure_bundle, tau_times, time_cut, sig_D, omega_c, M]\n",
    "    rms_test = np.inf\n",
    "  \n",
    "    for i in range(len(D_R_list)):\n",
    "        x0 = [D_R_list[i], a_LDF_list[i]]        \n",
    "        comb_solv = scipy.optimize.minimize(objective_function_comb, x0, args=uptake_bundle, method='Nelder-Mead', options={'xatol':1e-4, 'fatol':1e-4}) #})\n",
    "        rms = comb_solv['fun']        \n",
    "        if rms < rms_test:\n",
    "            rms_test = rms\n",
    "            D_R = comb_solv['x'][0]\n",
    "            a_LDF = comb_solv['x'][1]\n",
    "            k_r = a_LDF * D_R\n",
    "            omega = omega_c / D_R\n",
    "            \n",
    " \n",
    "    D_R_set.append(D_R)\n",
    "\n",
    "    aLDF_set.append(a_LDF)\n",
    "    P_U_0_check.append(P_U_0)\n",
    "    P_inf_check.append(P_inf) \n",
    "    P_D_0_check.append(P_D_0)\n",
    "\n",
    "    sig_D_comb = combined_diffusion(omega, D_R, delta_gamma, a_LDF, pressure_bundle, tau_times, M)\n",
    "    plt.plot(tau_times, sig_D_comb, c='purple')\n",
    "    fitted_data.append(sig_D_comb)\n",
    "    plt.show()\n",
    "\n",
    "    print(f\"Equilibrium pressure (mmHg) {P_inf}\")\n",
    "    print(f\"Intracrystalline diffusion coefficient (m2s-1) {D_R*L_scale[Chamber_select]**2}\")\n",
    "    print(f\"Surface barrier (ms-1) {a_LDF*D_R*L_scale[Chamber_select]}\")\n",
    "    print(f\"RMS error {rms}\")\n",
    "\n",
    "D_R_overall.append(D_R_set)\n",
    "aLDF_overall.append(aLDF_set)\n",
    "rms_overall.append(rms_set)\n",
    "P_U_0_overall.append(P_U_0_check)\n",
    "P_D_0_overall.append(P_D_0_check)\n",
    "P_inf_overall.append(P_inf_check)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D_R_overall_clean = []\n",
    "P_U_0_overall_clean = []\n",
    "P_inf_overall_clean = []\n",
    "P_D_0_overall_clean = []\n",
    "aLDF_overall_clean = []\n",
    "for i in range(1):\n",
    "    D_R_clean = []\n",
    "    P_U_0_clean = []\n",
    "    aLDF_clean = []\n",
    "    P_inf_clean = []\n",
    "    P_D_0_clean = []\n",
    "    for j in range(len(D_R_overall[i])):\n",
    "        D_R_clean.append(D_R_overall[i][j])\n",
    "        P_U_0_clean.append(P_U_0_overall[i][j])\n",
    "        P_D_0_clean.append(P_D_0_overall[i][j])\n",
    "        P_inf_clean.append(P_inf_overall[i][j])\n",
    "        aLDF_clean.append(aLDF_overall[i][j])\n",
    "    D_R_overall_clean.append(D_R_clean)\n",
    "    P_U_0_overall_clean.append(P_U_0_clean)\n",
    "    P_D_0_overall_clean.append(P_D_0_clean)\n",
    "    P_inf_overall_clean.append(P_inf_clean)\n",
    "    aLDF_overall_clean.append(aLDF_clean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitted_data\n",
    "expts_data\n",
    "fitted_data\n",
    "colors = cm.rainbow(np.linspace(0,1,len(uptake_zip)))\n",
    "def sequential_color_plot(intensity_data, position_data, fit_data, colors=colors, data_suffix=data_suffix):\n",
    "    plt.figure(figsize=(8,8))\n",
    "    for i, c in zip(range(len(intensity_data)), colors):\n",
    "        plt.scatter(position_data[i], intensity_data[i], color=c)\n",
    "        plt.plot(tau_times, fit_data[i], color=c)\n",
    "    plt.title(f'{data_suffix} : Chamber {Chamber_select+1}')\n",
    "    plt.yscale('log')\n",
    "    plt.ylim(1e-2, 1)\n",
    "    plt.xlim(0, 20)\n",
    "    plt.xlabel('$t$ [$s$]')\n",
    "    plt.ylabel('fractional uptake')\n",
    "sequential_color_plot(expts_data[:], timeset_data[:], fitted_data[:])\n",
    "\n",
    "\n",
    "print(\"Remember to cite Brandani, S.; Brandani, F.; Mangano, E.; Pullumbi, P. Using a Volumetric Apparatus to Identify and Measure the Mass Transfer Resistance in Commercial Adsorbents. Microporous and Mesoporous Materials 2020, 304, 109277.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.cm as cm\n",
    "colors = cm.Set1(np.linspace(0,1,len(P_U_0_overall_clean)))\n",
    "\n",
    "def sequential_color_plot(intensity_data, position_data, colors=colors, data_suffix=data_suffix):\n",
    "    fig, ax = plt.subplots()\n",
    "    for i, c in zip(range(len(intensity_data)), colors):\n",
    "        ax.scatter(position_data[i], intensity_data[i], color=c)\n",
    "    ax.set_title(f'{data_suffix} : Chamber {Chamber_select+1}')\n",
    "    ax.set_xlabel('Uptake cell pressure (Torr)')\n",
    "    ax.set_ylabel('Intrinsic intracrystalline diffusion coefficient')\n",
    "    plt.yscale(\"log\")\n",
    "    plt.ylim(0.001, 10)\n",
    "    plt.xlim(0, max_pressure[Chamber_select]+20)\n",
    "sequential_color_plot(D_R_overall_clean, P_inf_overall_clean)\n",
    "\n",
    "#######################\n",
    "#\n",
    "# plot of dimensionless intracrystalline diffusion coefficients against the chamber cell pressure\n",
    "#\n",
    "#######################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Intrinsic intracrystalline diffusivity\\n\")\n",
    "print(\"Pressure (Torr)  /  D/R^2 (s^-1)\\n\")\n",
    "for i in range(len(D_R_overall_clean[0])):\n",
    "    print(f\"{P_inf_overall_clean[0][i]} {D_R_overall_clean[0][i]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.cm as cm\n",
    "colors = cm.Set1(np.linspace(0,1,len(P_U_0_overall_clean)))\n",
    "#uptake_zip = uptake_zip[0]\n",
    "#time_zip = time_zip[0]\n",
    "def sequential_color_plot(intensity_data, position_data, colors=colors, data_suffix=data_suffix):\n",
    "    fig, ax = plt.subplots()\n",
    "    for i, c in zip(range(len(intensity_data)), colors):\n",
    "        ax.scatter(position_data[i], intensity_data[i], color=c)\n",
    "    ax.set_title(f'{data_suffix} : Chamber {Chamber_select+1}')\n",
    "    ax.set_xlabel('Uptake cell pressure (Torr)')\n",
    "    ax.set_ylabel('surface/intracrystalline ratio')\n",
    "    plt.yscale(\"log\")\n",
    "    plt.ylim(0.001, 100)\n",
    "    plt.xlim(0, max_pressure[Chamber_select])\n",
    "sequential_color_plot(aLDF_overall_clean, P_inf_overall_clean)\n",
    "\n",
    "#######################\n",
    "#\n",
    "# plot of ratio between intracrystalline diffusion coefficients and surface barriers against the chamber cell pressure\n",
    "#\n",
    "#######################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Surface / intracrystalline ratio\\n\")\n",
    "print(\"Pressure (Torr)  /  alpha_LDF\\n\")\n",
    "for i in range(len(aLDF_overall_clean[0])):\n",
    "    print(f\"{P_inf_overall_clean[0][i]} {aLDF_overall_clean[0][i]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.cm as cm\n",
    "colors = cm.Set1(np.linspace(0,1,len(P_U_0_overall_clean)))\n",
    "#uptake_zip = uptake_zip[0]\n",
    "#time_zip = time_zip[0]\n",
    "def sequential_color_plot(intensity_data, position_data, colors=colors, data_suffix=data_suffix):\n",
    "    fig, ax = plt.subplots()\n",
    "    for i, c in zip(range(len(intensity_data)), colors):\n",
    "        ax.scatter(position_data[i], np.array(intensity_data[i])*(L_scale[Chamber_select]**2), color=c)\n",
    "    ax.set_title(f'{data_suffix} : Chamber {Chamber_select+1}')\n",
    "    ax.set_xlabel('Uptake cell pressure (Torr)')\n",
    "    ax.set_ylabel('intracrystalline diffusion coefficient (m2/s)')\n",
    "    plt.yscale(\"log\")\n",
    "    plt.ylim(1e-14, 1e-8)\n",
    "    plt.xlim(0, max_pressure[Chamber_select])\n",
    "sequential_color_plot(D_R_overall_clean, P_inf_overall_clean)\n",
    "\n",
    "#######################\n",
    "#\n",
    "# plot of intracrystalline diffusion coefficients against the chamber cell pressure\n",
    "#\n",
    "#######################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Intracrystalline diffusivity\\n\")\n",
    "print(\"Pressure (Torr)  /  D (m^2 / s) \\n\")\n",
    "for i in range(len(D_R_overall_clean[0])):\n",
    "    print(f\"{P_inf_overall_clean[0][i]} {D_R_overall_clean[0][i]*L_scale[Chamber_select]**2}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.cm as cm\n",
    "colors = cm.Set1(np.linspace(0,1,len(P_U_0_overall_clean)))\n",
    "def sequential_color_plot(intensity_data, position_data, colors=colors, data_suffix=data_suffix):\n",
    "    fig, ax = plt.subplots()\n",
    "    for i, c in zip(range(len(intensity_data)), colors):\n",
    "        big_array = np.array(aLDF_overall_clean[i])*np.array(D_R_overall_clean[i])*L_scale[Chamber_select]\n",
    "        ax.scatter(position_data[i], np.array(big_array), color=c)\n",
    "    ax.set_title(f'{data_suffix} : Chamber {Chamber_select+1}')\n",
    "    ax.set_xlabel('Initial uptake cell pressure (Torr)')\n",
    "    ax.set_ylabel('surface permeability (m/s)')\n",
    "    plt.yscale(\"log\")\n",
    "    plt.ylim(1e-9, 1e-5)\n",
    "    plt.xlim(0, max_pressure[Chamber_select])\n",
    "sequential_color_plot(aLDF_overall_clean, P_inf_overall_clean)\n",
    "\n",
    "#######################\n",
    "#\n",
    "# plot of the surface permeabilities against the chamber cell pressure\n",
    "#\n",
    "#######################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Surface barrier\\n\")\n",
    "print(\"Pressure (Torr)  /  k (m / s) \\n\")\n",
    "for i in range(len(D_R_overall_clean[0])):\n",
    "    print(f\"{P_inf_overall_clean[0][i]} {aLDF_overall_clean[0][i]*D_R_overall_clean[0][i]*L_scale[Chamber_select]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Intrinsic surface barrier\\n\")\n",
    "print(\"Pressure (Torr)  /  k (1 / s) \\n\")\n",
    "for i in range(len(D_R_overall_clean[0])):\n",
    "    print(f\"{P_inf_overall_clean[0][i]} {aLDF_overall_clean[0][i]*D_R_overall_clean[0][i]}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}