# VO2 Modelling based on "Kinetic analysis of oxygen dynamics under variable work rate" (https://www.uni-konstanz.de/mmsp/pubsys/publishedFiles/ArBeBr17.pdf)

# open .mat file (taken from PB19_T0202_cwr300_cut_151123.mat)
data <- R.matlab::readMat('VO2Data.mat')
Power <- data[[1]][[1]][1,]
Time <- data[[1]][[2]][1,]
VO2 <- data[[1]][[3]][1,]

OneExpParams <- R.matlab::readMat('OneExpModel.mat')
Vbase_one <- OneExpParams[[1]][[3]][1] # or unlist() first
A1_one <- OneExpParams[[1]][[4]][1]
tau1_one <- OneExpParams[[1]][[5]][1]
T1_one <- OneExpParams[[1]][[6]][1]

OneExpModell <- Vbase_one + fBasics::Heaviside(Time-T1_one)*A1_one*(1-exp(-(Time-T1_one)/tau1_one))

TwoExpParams <- R.matlab::readMat('TwoExpModel.mat')
Vbase_two <- TwoExpParams[[1]][[3]][1]
A1_two <- TwoExpParams[[1]][[4]][1]
A2_two <- TwoExpParams[[1]][[5]][1]
tau1_two <- TwoExpParams[[1]][[6]][1]
tau2_two <- TwoExpParams[[1]][[7]][1]
T1_two <- TwoExpParams[[1]][[8]][1]
T2_two <- TwoExpParams[[1]][[9]][1]

TwoExpModell <- Vbase_two + fBasics::Heaviside(Time-T1_two)*A1_two*(1-exp(-(Time-T1_two)/tau1_two)) + fBasics::Heaviside(Time-T2_two)*A2_two*(1-exp(-(Time-T2_two)/tau2_two))
TwoExpFirst <- Vbase_two + fBasics::Heaviside(Time-T1_two)*A1_two*(1-exp(-(Time-T1_two)/tau1_two))

# plot data
library(ggplot2)
source("defines.R")
ggplot() + plot_style +
  geom_line(aes(x = Time, y = Power, color="Power"), size = 5) + 
  geom_line(aes(x = Time, y = VO2*100, color="VO2"), size = 5) +
  geom_line(aes(x = Time, y = OneExpModell*100, color="OneExpModell"), size = 3) +
  geom_line(aes(x = Time, y = TwoExpModell*100, color="TwoExpModell"), size = 3) +
  geom_line(aes(x = Time, y = TwoExpFirst*100, color="TwoExpFirst"), size = 2, linetype = "longdash") +
  scale_color_manual(name = "Values", values = c("Power"="blue", "VO2"="red", "OneExpModell"="grey", "TwoExpModell"="yellow", "TwoExpFirst"="black")) + 
  scale_y_continuous(breaks=seq(0, 300, by=100), limits=c(50,300), sec.axis = sec_axis(~./100, name = "VO2"))

# Euler Verfahren
# Sekundensampling ==> h = Time[2]-Time[1] = 1
x <- numeric(length(Time))
x[1] <- 0
for (i in 2:length(Time))
{
  x[i] = x[i-1] + 1/tau1_one * (A1_one - x[i-1])
}
OneEuler <- Vbase_one + approx(x = Time, y = x, xout = (Time-T1_one), yleft = 0)[[2]]

# Euler Verfahren
# Sekundensampling ==> h = Time[2]-Time[1] = 1
x1 <- numeric(length(Time))
x1[1] <- 0
for (i in 2:length(Time))
{
  x1[i] = x1[i-1] + 1/tau1_two * (A1_two - x1[i-1])
}
x2 <- numeric(length(Time))
x2[1] <- 0
for (i in 2:length(Time))
{
  x2[i] = x2[i-1] + 1/tau2_two * (A2_two - x2[i-1])
}
TwoEuler <- Vbase_two + approx(x = Time, y = x1, xout = (Time-T1_two), yleft = 0)[[2]] + approx(x = Time, y = x2, xout = (Time-T2_two), yleft = 0)[[2]]

# plot data
ggplot() + plot_style +
  geom_line(aes(x = Time, y = Power, color="Power"), size = 5) + 
  geom_line(aes(x = Time, y = VO2*100, color="VO2"), size = 5) +
  geom_line(aes(x = Time, y = OneExpModell*100, color="OneExpModell"), size = 3) +
  geom_line(aes(x = Time, y = OneEuler*100, color="OneEuler"), size = 3, linetype = "dashed") +
  geom_line(aes(x = Time, y = TwoExpModell*100, color="TwoExpModell"), size = 3) +
  geom_line(aes(x = Time, y = TwoExpFirst*100, color="TwoExpFirst"), size = 2, linetype = "longdash") +
  geom_line(aes(x = Time, y = TwoEuler*100, color="TwoEuler"), size = 3, linetype = "dashed") +
  scale_color_manual(name = "Values", values = c("Power"="blue", "VO2"="red", "OneExpModell"="grey", "OneEuler"="green", "TwoExpModell"="yellow", "TwoExpFirst"="black", "TwoEuler"="brown")) + 
  scale_y_continuous(breaks=seq(0, 300, by=100), limits=c(50,300), sec.axis = sec_axis(~./100, name = "VO2"))