# 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 limists as datasheet curve ("xlim" and "ylim")
# 2) find out datasheet curve points, and write it out (to "pointsdata")
# 3) measure the sensor resistance (set it to "mres")
# 4) set (to "mres") the know amount of gas for the resistance measure of the previous step

library(data.table)

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

#set limists as datasheet curve graph here
#ex.
#  xlim = c(10, 1000)
#  ylim = c(0.1, 10)
xlim = c(0.1, 1000)
ylim = c(0.1, 1000)

#set measured resistance (ohm) at a know gas amount (ppm)
#ex.
#  mres = 26954
#  mppm = 392
mres = 20000
mppm = 100

#extraced datasheet curve points
#each line it's a point on cartesian coordinate system
#ex.
#4.2323, 19.2312
#5.2313, 28.9231
#the useful WebPlotDigitizer app can help you extract points from the graph
pointsdata = "
YOUR_VALUES_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)

#display fit curve coefficients (log-log scale)
cat("\nFit coefficients\n")
coef(fit)
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)

#display fitrev curve coefficients (log-log scale)
cat("\nFit coefficients for flipped curve\n")
coef(fitrev)
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)

#display fit curve coefficients (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)

#display fitrev curve coefficients (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 Ro
Ro = mres*(coef(fitrev)["a"]/mppm)^(1/coef(fitrev)["b"])
cat("\nEstimated Ro\n")
cat(Ro)
cat("\n")

#estimate min Rs/Ro
minppm = xlim[2]
minRsRo = (minppm/coef(fitrev)["a"])^(1/coef(fitrev)["b"])
cat("\nEstimated min Rs/Ro\n")
cat(minRsRo)
cat("\n")

#estimate max Rs/Ro
maxppm = xlim[1]
maxRsRo = (maxppm/coef(fitrev)["a"])^(1/coef(fitrev)["b"])
cat("\nEstimated max Rs/Ro\n")
cat(maxRsRo)
cat("\n")