# Simulating EIS in R

I've been asked on a few occasions about how I go about simulating impedance spectra in the R language. It's essentially based on straightforward maths with some functions defined to make the more complicated calculations simple. This page here gives an example of the approach I normally use, and the approach that underpins the other simulations shown on this website.

First, we can create a function to create a vector of values for the angular frequency, $$\omega$$, based on minimum and maximum values of log frequency and a number of points per decade:

omega_range <- function(low.logf, high.logf, p.per.dec) {
2 * pi * 10^seq(low.logf, high.logf, length.out = (p.per.dec * (high.logf - low.logf)) + 1)
}


### Equivalent circuit elements

Now we can create functions which calculate the complex impedance for each of the individual circuit elements. If you're unclear on how these are calculated, then I recommend looking at the other pages in this section for the mathematical definitions.

A function for the impedance of a resistor is not strictly necessary but I find it helps for readability later:

Z_R <- function(R) {
complex(real = R, imaginary = 0)
}


All the following elements are functions of angular frequency. First the capacitor:

Z_C <- function(C, omega) {
complex(real = 0, imaginary = -1/(omega * C))
}


Constant phase element:

Z_Q <- function(Q, n, omega) {
1/(Q * (complex(real = 0, imaginary = 1) * omega)^n)
}


Semi-infinite Warburg:

Z_W <- function(sigma, omega) {
j = complex(real = 0, imaginary = 1)
(1 - j) * sigma * omega^-0.5
}


Finite length Warburg:

Z_FLW <- function(Z0, tau, omega) {
j = complex(real = 0, imaginary = 1)
Z0 * (j * omega * tau)^-0.5 * tanh((j * omega * tau)^0.5)
}


Finite space Warburg (note: you must have the pracma package installed, which contains the required coth() function):

Z_FSW <- function(Z0, tau, omega) {
require(pracma)
j = complex(real = 0, imaginary = 1)
Z0 * (j * omega * tau)^-0.5 * coth((j * omega * tau)^0.5)
}


### Simple simulation

First we need to define which frequencies we will calculate, using the omega_range() function defined earlier. For example:

omega <- omega_range(low.logf = -0.7, high.logf = 5.2, p.per.dec = 8)


And to make life easier, we can create a function so we can add the impedances of circuit elements in parallel:

par <- function(a, b) {
1/((1/a) + (1/b))
}


Then simulating a typical equivalent circuit might look something like this, using the Randles circuit as an example:

z_rand <- Z_R(5) + par((Z_R(45) + Z_W(40, omega)), Z_C(2E-5, omega))


The z_rand variable is a vector of complex numbers which we can't plot very easily:

str(z_rand)

##  cplx [1:49] 85.6-35.9i 80.9-31.2i 76.7-27.1i ...


So we can create a data frame with the numbers we want using the Re() and Im() functions in R for extracting the real and imaginary parts from the complex numbers.

df <- data.frame(Re = Re(z_rand), Im = -Im(z_rand), f = omega/(2 * pi), omega = omega)

str(df)

##  'data.frame':   49 obs. of  4 variables:
##   $Re : num 85.6 80.9 76.7 73.2 70.1 ... ##$ Im   : num  35.9 31.2 27.1 23.6 20.6 ...
##   $f : num 0.2 0.265 0.351 0.466 0.619 ... ##$ omega: num  1.25 1.66 2.21 2.93 3.89 ...


And now it's just a question of creating the plot:

library(ggplot2)
ggplot(df, aes(x = Re, y = Im)) +
geom_point() + coord_fixed() The coord_fixed() function for ggplot is extremely useful here in ensuring the proportionality of the axes in the Nyquist plot. 