# 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))

# vof <- function(t, base, A, T, tau) {
#   return(base + sum(fBasics::Heaviside(t-T)*A*(1-exp(-(t-T)/tau))))
# }
# oneexpR <- sapply(Time, vof, base=Vbase_one, A=A1_one, T=T1_one, tau=tau1_one)
# sum(OneExpModell-oneexpR)
# Ak <- c(TwoExpParams[[1]][[4]][1], TwoExpParams[[1]][[5]][1])
# tauk <- c(TwoExpParams[[1]][[6]][1], TwoExpParams[[1]][[7]][1])
# Tk <- c(TwoExpParams[[1]][[8]][1], TwoExpParams[[1]][[9]][1])
# twoexpR <- sapply(Time, vof, base=Vbase_two, A=Ak, T=Tk, tau=tauk)
# sum(TwoExpModell-twoexpR)

# 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"))