MQ gas sensor datasheet correlation function coefficients estimation and constants estimations

# An R script to estimate MQ gas sensors correlation curve and compute Ro, min and max Rs/Ro
#
# Copyright (c) Davide Gironi, 2016
#
# Released under GPLv3.
# Please refer to LICENSE file for licensing information.

# How to use this script:
# 1) set limits as datasheet curve ("xlim" and "ylim")
#    ex.
#      xlim = c(10, 1000)
#      ylim = c(0.1, 10)
# 2) find out datasheet curve points, and write it out (to "pointsdata")
#    each line it's a point on cartesian coordinate system
#    the useful WebPlotDigitizer app can help you extract points from the graph
#    ex.
#      pointsdata = "
#        10.052112405371744, 2.283698378106183
#        20.171602728600178, 1.8052797165878915
#        30.099224396434586, 1.5715748803154423
#        50.09267987761949, 1.3195287228519417
#        80.38812026903305, 1.1281218760133969
#        90.12665922665023, 1.0815121769656304
#        100.52112405371739, 1.0430967861855598
#        199.62996638292853, 0.8000946404902397
#      "
# 3) optional for Ro estimation: measure the sensor resistance (set it to "mres" ohm value) at a know amount of gas
#    set it to 0 if you do not need the Ro estimation
#    ex.
#      mres = 26954
# 4) optional for Ro estimation: set the know amount of gas for the resistance measure of the previous step (to "mppm")
#    set it to 0 if you do not need the Ro estimation
#    ex.
#      mppm = 392
# 5) optional for min-max Rs/Ro estimation: set the minand max amount of gas the sensor will react to (as "minppm" and "maxppm")
#    set it to 0 if you do not need the min-max Rs/Ro estimation
#    ex.
#      minppm = 10
#      maxppm = 200
library(data.table)

#remove old variables
rm(list=ls())

#set input values
xlim = c(0.1, 10)
ylim = c(0.1, 10)
minppm = 0
maxppm = 0
mres = 0
mppm = 0
pointsdata = "
YOUR_POINTS_HERE
"

#load points using fread
setnames(points <- fread(pointsdata, sep=",", sep2="\n"), c("x","y"))

#set named list of points, and swapped list of points
#points will be used to plot and compute values as datasheet figure
#pointsrev will be used to plot and compute values for the correlation function, it's the datasheet figure with swapped axis
x <- as.vector(points[,x])
y <- as.vector(points[,y])
points = list(x=x, y=y)
pointsrev = list(x=y, y=x)

#the nls (Nonlinear Least Squares) it's used to perform the power regression on points
#in order to work, nls needs an estimation of staring values
#we use log-log slope estimation to find intitial values

#estimate fit curve initial values
xfirst = head(points$x, n=1)
xlast = tail(points$x, n=1)
yfirst = head(points$y, n=1)
ylast = tail(points$y, n=1)
bstart= log(ylast/yfirst)/log(xlast/xfirst)
astart = yfirst/(xfirst^bstart)
#perform the fit
fit <- nls("y~a*x^b", start=list(a=astart,b=bstart), data=points)

#estimate fitref curve initial values
xfirstrev = head(pointsrev$x, n=1)
xlastrev = tail(pointsrev$x, n=1)
yfirstrev = head(pointsrev$y, n=1)
ylastrev = tail(pointsrev$y, n=1)
bstartrev = log(ylastrev/yfirstrev)/log(xlastrev/xfirstrev)
astartrev = yfirstrev/(xfirstrev^bstartrev)
fitrev <- nls("y~a*x^b", start=list(a=astartrev,b=bstartrev), data=pointsrev)

#plot fit curve (log-log scale)
fiteq = function(x){coef(fit)["a"]*x^(coef(fit)["b"])}
plot(points, log="xy", col="blue", xlab="ppm", ylab="Rs/Ro", xlim=xlim, ylim=ylim, panel.first=grid(equilogs=FALSE))
curve(fiteq, col="red", add=TRUE)

#plot fitrev curve (log-log scale)
fiteqrev = function(x){coef(fitrev)["a"]*x^(coef(fitrev)["b"])}
plot(pointsrev, log="xy", col="blue", xlab="Rs/Ro", ylab="ppm", xlim=ylim, ylim=xlim, panel.first=grid(equilogs=FALSE))
curve(fiteqrev, col="red", add=TRUE)

#plot fit curve (linear scale)
fiteq = function(x){coef(fit)["a"]*x^(coef(fit)["b"])}
plot(points, col="blue", xlab="ppm", ylab="Rs/Ro", panel.first=grid(equilogs=FALSE))
curve(fiteq, col="red", add=TRUE)

#plot fitrev curve (linear scale)
fiteqrev = function(x){coef(fitrev)["a"]*x^(coef(fitrev)["b"])}
plot(pointsrev, col="blue", xlab="Rs/Ro", ylab="ppm", panel.first=grid(equilogs=FALSE))
curve(fiteqrev, col="red", add=TRUE)

#estimate min Rs/Ro
cat("\nCorrelation function coefficients")
cat("\nEstimated a\n")
cat("  ")
cat(coef(fitrev)["a"])
cat("\nEstimated b\n")
cat("  ")
cat(coef(fitrev)["b"])
cat("\n")

#estimate min Rs/Ro
if (minppm != 0) {
    minRsRo = (maxppm/coef(fitrev)["a"])^(1/coef(fitrev)["b"])
    cat("\nEstimated min Rs/Ro\n")
    cat("  ")
    cat(minRsRo)
    cat("\n")
}

#estimate max Rs/Ro
if (maxppm != 0) {
    maxRsRo = (minppm/coef(fitrev)["a"])^(1/coef(fitrev)["b"])
    cat("\nEstimated max Rs/Ro\n")
    cat("  ")
    cat(maxRsRo)
    cat("\n")
}

#estimate Ro
if (mppm != 0 && mres != 0) {
    Ro = mres*(coef(fitrev)["a"]/mppm)^(1/coef(fitrev)["b"])
    cat("\nEstimated Ro\n")
    cat("  ")
    cat(Ro)
    cat("\n")
}