fiteo.R

In mendivilg/eisr: eisR - Fitting EIS

# FUNCIONES CON NÚMEROS COMPLEJOS The functions Re, Im, Mod, Arg and Conj have their usual interpretation
# as returning the real part, imaginary part, modulus, argument and complex conjugate for complex values.
# The modulus and argument are also called the polar coordinates. If z = x + i y with real x and y, for r
# = Mod(z) = √(x^2 + y^2), and φ = Arg(z), x = r*cos(φ) and y = r*sin(φ).



j <- 1i  # unidad compleja imaginaria



# Impedancia de un capacitor
zc <- function(f, c) {
  #-j/(2 * pi * f * c)  # o
  1/(j*2*pi*f*c)
}

# Impedancia de un inductor
zl <- function(f, l) {
  j * 2 * pi * f * l
}

# Impedancia Warburg
zw <- function(f, Yo) {
  1/(sqrt(j * 2 * pi * f) * Yo)
}

# Impedancia CPE -Elemento de fase constante Q
zq <- function(f, Yo, n) {
  #(j*2*pi*f)^(-n)/Yo
  1/(((j*2*pi*f)^n)*Yo)
}

# Impedancia FSW - Finite length diffusion elements.  The FSW (Finite Space Warburg) T
zt <- function(f, Yo, B) {
  omega <- 2 * pi * f
  argu <- B * sqrt(j * omega)
  (cosh(argu))/sinh(argu)/Yo * sqrt(j * omega)
}

# Impedancia FLW - Finite length diffusion elements.  The FLW (Finite Length Warburg) O
zo <- function(f, Yo, B) {
  omega <- 2 * pi * f
  argu <- B * sqrt(j * omega)
  tanh(argu)/Yo * sqrt(j * omega)
}

# Impedancia Gerisher G.
zg <- function(f, Yo, k) {
  omega <- 2 * pi * f
  (k + j * omega)^(-0.5)/Yo
}

# Impedancia Gerisher Fractal F.
zf <- function(f, Yo, k, alfa) {
  omega <- 2 * pi * f
  (k + j * omega)^(alfa)/Yo
}

# Calculo de dos impedancias en paralelo
paral <- function(imp1, imp2) (imp1*imp2)/(imp1+imp2)

## Genera barrido de frecuencias
sweep <- function(fi, ff, step){
  ## fi: frecuencia inicial
  ## ff: frecuencia final
  ## step: muestras por decada
  ## TO DO: falta modo lineal y paso/octava
  return(c(10^(seq(log10(fi), log10(ff), by = -(step/100))),ff))
}

# Minimal code for impedance fitting using the R
# Kiyoshi Kobayashi (National Institute for Materials Science)
# datos en DOI: 10.13140/RG.2.1.3540.2085, test imp data, Exmpl impedance fitting
################################################

RC <- function(frec, r,c, n)
  return(paral(r,  zq(frec, c, n)))

frec <- sweep(1000, 0.05,10)
x <- frec
impedancia <- 200 + RC(frec, 3150, 0.0005, 0.1)


# Import minpac.lm for Levenberg-Marquardt algorithm
library(minpack.lm)
library(nlstools)
library(eisr)
# Open dialog to select text data file
data <- openEIS()  # Archivo data en directorio de trabajo

# Separate the data into f, zr, and zi
f <- data[,2,1]
zr <- data[,3,1]
zi <- data[,4,1]

# Complex impedance data list
complejo <- 14 * zr + 1i*zi

# Declare impedance model function
# PARÁMETROS MODELO 1:
# p[1] : R1
# p[2] : CPE1 Yo
# p[3] : CPE1 n
# p[4] : R2
# p[5] : C1
# p[6] : R3
# p[7] : C2
# p[8] : R4

imp <- function(frec, p){
   z <- p[1]
   z <- z + paral(zq(frec,p[2], p[3]), p[4])
   z <- z + paral(zc(frec, p[5]), p[6])
   z <- z + paral(zc(frec, p[7]), p[8])
   return(z)
 }


# Declare the residual function
rsd <- function(p,z,x) {
  errz = imp(x,p)-z
  err = c(Re(errz), Im(errz))
  return (err)
}

# Input an initial guess of the fitting parameters
# p[1] : R1
# p[2] : CPE1 Yo
# p[3] : CPE1 n
# p[4] : R2
# p[5] : C1
# p[6] : R3
# p[7] : C2
# p[8] : R4

p0 <- c( 77, 4.7e-4, 0.78, 26, 3.9e-2, -9, 5.1e-1, 90)
linf <- c( 40, 4.7e-4, 0.78e-1, 0, 3.9e-6, 0, 5.1e-3, 0)
lsup <- c( 7700, 4.7, 8.78, 526, 3.9,123, 5.1, 3490)

# start fitting
nls.out <- nls.lm(p = p0, lower = linf, upper = lsup, fn = rsd, z = complejo, x = f, control =
                    nls.lm.control(nprint = 1, maxfev = 10000, maxiter = 1000, ftol = 0.000001))


summary(nls.out)
aa <- summary(nls.out)

pfit <- coef(nls.out)
perr <- aa[[8]][11:20]
fcov <- vcov(nls.out)
RSS <- sum(abs(imp(pfit,f)-zc)*abs(imp(pfit,f)-zc))

while (aa[[5]] != 1){
  p0 <- pfit
  nls.out <- nls.lm(p = p0, lower = pl, upper = pu, fn = rsd, z = complejo, x = f, control =
                      nls.lm.control(nprint = 1, maxfev = 1000, maxiter = 100, ftol = 0.000001))
  summary(nls.out)
  aa <- summary(nls.out)
}

summary(nls.out)
aa <- summary(nls.out)
pfit <- coef(nls.out)
perr <- aa[[8]][11:20]
RSS <- sum(abs(imp(pfit,f)-complejo)*abs(imp(pfit,f)-complejo))
fcov <- vcov(nls.out)


zc_cal <- imp(pfit,f)
zr_cal <- Re(zc_cal)
zi_cal <- Im(zc_cal)

plot(f,zr,log="x",ann=F)
#par(new=T)
plot(f,zr_cal,log="x",  type="l")

#dev.new()
plot(f,zi,log="x", xlim=c(10^-2, 10^6),ylim=c(-2.5, 0.5),ann=F)
#par(new=T)
lines(f,zi_cal,log="x", xlim=c(10^-2, 10^6),ylim=c(-2.5, 0.5),type="l", col="red")

#dev.new()
plot(zr,-zi,xlim=c(10,25),ylim=c(-2, 5),ann=F, asp = 1)
#par(new=T)
lines(zr_cal,-zi_cal, xlim=c(10,25),ylim=c(-2, 5), type="l", asp = 1, col ="red")

pfit
perr
RSS
fcov