<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Galvanostatic Intermitten Titration Technique (GITT)
Step1:  Some functions
Custom step-table for GITT
Need to allow for several ocv steps with same step number in the same cycle. This is not calculated by default and we therefore have to (re-) run the make_step_table method selecting
all_steps=True.
Step2:  Checking
Step3:  Initial exploration
Overview of the full test
Curve with labels
Step4:  Holomap (with selector)
Step5:  Curve and HoloMap
Step6:  Processing a cycle
As we have seen, only some of the cycles contains the GITT curves. We need to select one of those. Lets start with cycle 5
Step7:  Some calculations are needed...
Step8:  The following equation to estimate the Diffusion Coefficients appears in a paper
Step11:  Comparing another reference (below) I assume that M is 1/number_of_moles ?
A function to automatically extract diffusion constant(s) from the steptable
This function is based on the Metrohm Autolab document. It will need some adjustments.
Step12:  Using the function to extract Diffusion constants vs capacity
Step13:  TODO
make a plot where we show each point that is used in the calculations
tweak the method so that we get realistic results
find a paper to compare with
Picking out titration curves
<ASSISTANT_TASK:>
Python Code:
%load_ext autoreload
%autoreload 2
import os
from pathlib import Path
from pprint import pprint
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cellpy
from cellpy import prms
from cellpy import prmreader
from cellpy import cellreader
from cellpy.utils import ocv_rlx
import holoviews as hv
%matplotlib inline
hv.extension("bokeh")
######################################################################
##                                                                  ##
##                       development                                ##
##                                                                  ##
######################################################################
if os.name == "nt":
    # Use these when working on my work PC:
    raw_data_path = r"C:\Scripting\MyFiles\development_cellpy\dev_data\gitt"
    out_data_path = r"C:\Scripting\MyFiles\development_cellpy\out"
else:
    # Use these when working on my MacBook:
    raw_data_path = "/Users/jepe/scripting/cellpy/dev_data/gitt"
    out_data_path = "/Users/jepe/scripting/cellpy/dev_data/out"
raw_data_path = Path(raw_data_path)
out_data_path = Path(out_data_path)
print(" SETTING SOME PRMS ".center(80, "="))
prms.Paths.db_filename = "cellpy_db.xlsx"
prms.Paths.cellpydatadir = out_data_path
prms.Paths.outdatadir = out_data_path
prms.Paths.rawdatadir = raw_data_path
prms.Paths.db_path = out_data_path
prms.Paths.filelogdir = out_data_path
pprint(prms.Paths)
pd.set_option("display.max_rows", 30)
pd.set_option("display.max_columns", 500)
pd.set_option("display.width", 1000)
fn = "20190403_cen59_04_rateGITT_01"
mass = 0.3  # mg
resn = fn + ".res"
cellpyn = fn + ".h5"
filename = prms.Paths.rawdatadir / resn
cellpyname = prms.Paths.cellpydatadir / cellpyn
cell = cellreader.get(filename, mass=mass, logging_mode="INFO")
cell.save(cellpyname)
cell = cellreader.get(cellpyname)
# print(cell)
Explanation: Galvanostatic Intermitten Titration Technique (GITT)
End of explanation
steptable = cell.make_step_table(all_steps=True).cell.steps
Explanation: Some functions
Custom step-table for GITT
Need to allow for several ocv steps with same step number in the same cycle. This is not calculated by default and we therefore have to (re-) run the make_step_table method selecting
all_steps=True.
End of explanation
steptable.loc[steptable.cycle == 5, ["index", "step", "ustep", "type"]].head(4)
Explanation: Checking
End of explanation
# retrieving the voltage vs time curve
t, v = cell.get_timestamp(), cell.get_voltage()
all_cycs = hv.Curve(
    (t, v), ("t", "time (sec)"), ("v", "voltage (v)"), label="voltage-time"
).opts()
# creating labels
cycle_label_df = steptable.drop_duplicates("cycle")
cycle_label_df["label"] = cycle_label_df.cycle.apply(str)
cycle_labels = hv.Labels(
    (cycle_label_df.test_time_first, cycle_label_df.voltage_first, cycle_label_df.label)
).opts(
    text_align="left",
)
Explanation: Initial exploration
Overview of the full test
Curve with labels
End of explanation
# creating a dictionary of curves (for each cycle)
cycs_dict = dict()
for c in cell.get_cycle_numbers():
    t = cell.get_timestamp(cycle=c)
    t = t - np.amin(t.values)  # setting first point to t=0
    curve = hv.Curve(
        (t, cell.get_voltage(cycle=c)),
        ("time", "time (seconds)"),
        ("voltage", "voltage (v vs. Li/Li+)"),
    )
    cycs_dict[c] = curve
# creating a holomap object
hmap = hv.HoloMap(cycs_dict, "cycle")
Explanation: Holomap (with selector)
End of explanation
%%opts Curve [width=800, xformatter="%6.0f" tools=["hover"]]
(all_cycs * cycle_labels + hmap).cols(1)
Explanation: Curve and HoloMap
End of explanation
dfdata = cell.cell.raw
cyc5 = dfdata.loc[dfdata.Cycle_Index == 5, :]
gt5 = steptable.loc[steptable.cycle == 5, :]
cyc5curve = hv.Curve(cyc5, ("Test_Time", "time"), ("Voltage", "voltage")).opts(
    color="grey", alpha=0.5
)
cyc5points = hv.Scatter(cyc5, ("Test_Time", "time"), ("Voltage", "voltage")).opts(
    size=5, fill_alpha=0.2
)
slabels5 = hv.Labels(
    (gt5.test_time_first, gt5.voltage_first, gt5.ustep.apply(str))
).opts(text_align="left")
spikes5 = hv.Spikes(gt5.test_time_first).opts(
    line_alpha=0.3,
    spike_length=1,
)
%%opts Curve [width=1000, height=600, tools=['hover']]
print("the cycle starts where the label is")
cyc5curve * cyc5points * slabels5 * spikes5
Explanation: Processing a cycle
As we have seen, only some of the cycles contains the GITT curves. We need to select one of those. Lets start with cycle 5
End of explanation
gt5.head()
# approximating electrode-electrolyte contact area
diameter_si = 200.0 * 10e-7  # cm
rho_si = 2.32  # g/cm3
_mass = mass / 1000
area = (2 * 3 / diameter_si) * (_mass / rho_si)  # CHECK THIS
print(f"Diameter: {diameter_si:8.6f} cm")
print(f"Mass:     {mass:8.6f} mg")
print(f"Calculated contact area (ideal case): {area:8.6f} cm2")
# For reference (v, a, m is pr Si particle)
v = (4 / 3) * np.pi * ((diameter_si / 2) ** 3)
m = rho_si * v
n = mass / m
a = 4 * np.pi * (diameter_si / 2) ** 2
# estimate some values
# electrolyte-electrode area (cm2)
radius = 0.75  # cm
roughness = 1.0  # a dimensionless parameter ranging from 1 to a lot.
print("Should probably use BET to find out this?")
area = (1 * roughness) * np.pi * (radius ** 2)
print(f"S: {area:6.2f} cm2")
print(
    "You should compare this number to what was calculated in the cell above (the theoretical surface area for the set of same sized nanoparticles)"
)
# number of moles
# 1 mol Si weighs 28.0855 g
_mass = mass / 1000  # convert mass from mg to g
number_of_moles = _mass / 28.0855
print(f"n_m: {number_of_moles:8.6f} mol")
# molar volume
# The 2006 CODATA recommended value for the molar volume of silicon is 12.0588349(11)×10−6 m3/mol, with a relative standard uncertainty of 9.1×10−8
molar_volume = 12.06  # cm3/mol
print(f"V_m: {molar_volume:6.2f} cm3/mol")
Explanation: Some calculations are needed...
End of explanation
# Simple calculations by reading values from the graph
m = mass / 1000
molar_mass = 28.0855  # g/mol
density = 2.32  # g/cm3
Vm = molar_mass / density  # cm3/mol
S = 1.7  # cm2
t = 2880  # s
M = 1
Es1 = 0.378
Es2 = 0.344
DEs = Es2 - Es1
Et11 = 0.378  # starting the current pulse
Et12 = 0.309  # after (possible) IR drop
Et2 = 0.28
DEt1 = Et2 - Et11
DEt2 = Et2 - Et12
p1 = (4 / (np.pi * t)) * (m * Vm / (M * S)) ** 2
p21 = (DEs / DEt1) ** 2
p22 = (DEs / DEt2) ** 2
D1 = p1 * p21
D2 = p1 * p22
print(f"Assuming no IR drop: {D1}\nWith IR drop: {D2}")
print(f"A: {p1}")
print(f"DEt: {DEt1}")
print(f"DEs: {DEs}")
Explanation: The following equation to estimate the Diffusion Coefficients appears in a paper:
\begin{equation}
D = \frac{4}{\pi\cdot\tau} \cdot ( \frac {m \cdot V_m}{M \cdot S})^2 \cdot ( \frac {\Delta E_s}{ \Delta E_t})^2 \: for \: \tau << L^2 / D
\end{equation}
$ \Delta E_s $  = change of steady state voltage at the end of two sequential open-circuit relaxation periods
<br>
$ \Delta E_t $ = total change in the cell voltage $ E $ during the current pulse, neglecting the IR drop
<br>
$ V_s $ = mole volume of active material (cm^3/mol) (molar mass divided by the materials density)
<br>
$ M $ = ?
<br>
$ S $ = Surface area of the electrode (cm^2)
<br>
$ \tau $ = time of the current pulse
<br>
$ L $ = diffustion length  
Andrzej P. Nowak et al., Procedia Engineering 98 (2014) 8
End of explanation
def calc_A(n_m=0.0011, V_m=12.06, S=1.77):
    
    D = 4 /((pi)(tau)) * (n_m * V_m / S)^2 * (DEs / DEt)^2 = A /tau * (DEs / DEt)^2
    A = 4 /(pi) * (n_m * V_m / S)^2
    tau: duration of the current pulse (s)
    n_m: number of moles (mol)
    V_m: molar volume of electrode (cm3/mol)
    S: electrode-electrolyte contact area (cm2)
    DEs: steady state voltage change due to the current pulse (V)
    DEt: voltage change during the constant current pulse (eliminating the iR drop) (V)
    Ref.: application note from Metrohm Autolab b.v. pdf (BAT03)
    
    A = (4 / np.pi) * (n_m * V_m / S) ** 2
    return A
def auto_calc_D(
    steptable, cycle_number, A=1.0, tau=None, ustep_first=None, ustep_last=None
):
    Function for extracting diffusion constant(s) and inserting into steptable.
    D = 4 /((pi)(tau)) * (n_m * V_m / S)^2 * (DEs / DEt)^2 = A /tau * (DEs / DEt)^2
    A = 4 /(pi) * (n_m * V_m / S)^2
    OBS! not corrected for IR drop yet.
    
    st = None
    if cycle_number is None:
        # This function is intended to only work on a pr. cycle basis
        # to prevent users from "poluting" the steptable for "non_GITT" experiments.
        print("no cycle number given")
        return
    st = steptable[steptable.cycle == cycle_number]
    st = st[st.type.isin(["charge", "discharge", "ocvrlx_up", "ocvrlx_down"])]
    if st.empty:
        print("the given cycle is not found")
        return
    if ustep_first is not None:
        st = st[st.ustep >= ustep_first]
    if ustep_last is not None:
        st = st[st.ustep <= ustep_last]
    # used for finding DE
    n3 = st["voltage_last"].shift(periods=-3)
    n2 = st["voltage_last"].shift(periods=-2)
    n1 = st["voltage_last"].shift(periods=-1)
    n0 = st["voltage_last"]
    #  st["n3"] = n3
    #  st["n2"] = st["voltage_last"].shift(periods=-2)
    #  st["n1"] = n1
    #  st["n0"] = st["voltage_last"].shift(periods=0)
    # used for finding tau
    if tau is None:
        tau = st["step_time_last"] - st["step_time_first"]
    # used for validating if proper GITT step
    t0 = st["type"]
    t1 = st["type"].shift(periods=-1)
    t2 = st["type"].shift(periods=-2)
    t3 = st["type"].shift(periods=-3)
    st["valid_D"] = (t0 == t2) & (t1.str.contains("ocv"))
    # calculating
    st["DEt"] = (
        st["voltage_last"] - st["voltage_first"]
    )  # should remove the IR drop here (maybe calculate it based on the measured IR time the current?)
    st["DEs"] = n3 - n1
    st["D"] = A / tau * (st["DEs"] / st["DEt"]) ** 2
    return st
Explanation: Comparing another reference (below) I assume that M is 1/number_of_moles ?
A function to automatically extract diffusion constant(s) from the steptable
This function is based on the Metrohm Autolab document. It will need some adjustments.
End of explanation
# Lets just run it and see what happens
A = calc_A(n_m=0.0011, V_m=12.06, S=1.77)
print(f"A: {A}")
x = auto_calc_D(steptable, 5, A=A, ustep_first=27)
discharge = x.loc[(x.type == "discharge") & (x.valid_D)]
charge = x.loc[(x.type == "charge") & (x.valid_D)]
# discharge.head(25)
%%opts Curve [width=600 xformatter="%6.0f" tools=["hover"]]
X = discharge["discharge_avr"]*1_000_000/mass
Y = discharge["D"]
discharge_diffcurve = hv.Scatter((X,Y), ("capacity"), ("diffusion coef. (cm2 s-1)"), label="discharge").opts(size=12) * hv.Curve((X,Y)).opts(alpha=0.4)
X = charge["charge_avr"]*1_000_000/mass
Y = charge["D"]
charge_diffcurve = hv.Scatter((X,Y), ("capacity"), ("diffusion coef. (cm2 s-1)"), label="charge").opts(size=12) * hv.Curve((X,Y)).opts(alpha=0.4)
discharge_diffcurve * charge_diffcurve
Explanation: Using the function to extract Diffusion constants vs capacity
End of explanation
cycle = 5
titcurve_ustep = 28
# filtering wrt cycle number to get the Data_Points for the step
datapoints = x.loc[
    (x.type == "discharge") & (x.ustep == titcurve_ustep), ["point_first", "point_last"]
]
datapoints
# Using the data points to select the data from dfdata
first = datapoints.iloc[0, 0]
last = datapoints.iloc[0, 1]
dftit = dfdata.loc[first:last, :]
dftit.head()
dfcurve = (
    hv.Curve(dftit, ("Step_Time", "time"), ("Voltage", "voltage"))
    .opts(color="grey", alpha=0.5)
    .opts(width=1000, xformatter="%6.0f")
)
min_time = hv.VLine(100.0)
max_time = hv.VLine(1000.0)
a = hv.Arrow(x=500, y=0.31, text="fit", direction="v")
dfcurve * min_time * max_time * a
# plotting dE vs sqrt(t) and doing a linear regression to find the slope
dE = dftit["Voltage"]
t = np.sqrt(dftit["Step_Time"])
decurve = hv.Scatter((t, dE), "sqrt-time", "voltage").opts(
    width=1000, xformatter="%6.0f"
)
decurve
Explanation: TODO
make a plot where we show each point that is used in the calculations
tweak the method so that we get realistic results
find a paper to compare with
Picking out titration curves
End of explanation
<END_TASK>