#Data from tests on 30/10/2019 at met city - including recent crack formation event. 

setwd("/Users/sarcher/Desktop/R data/1_MOSAiC/")
AX <- read.csv("/Volumes/Data_1/MOSAIC_2019/co2_ch4_chamber_portable/2019_12_02/2019_12_02_CO2.csv", header = TRUE)
library(tidyverse)
library(dplyr)
#simplyfying the df name
head(AX)
str(AX)
nrow(AX)


#narrow down the axes by altering XR1, XR2 and YR1, YR2
#X and Y -axis range lower, higher
XR1 <- 6000
XR2 <- 6500
YR1 <- 410
YR2 <- 418

#Plotting the truncated axes 

My_Theme = theme(
  axis.title.x = element_text(size = 20),
  axis.text.x = element_text(size = 16),
  axis.title.y = element_text(size = 20),
  axis.text.y = element_text(size = 16),
  legend.title = element_text(size =20),
  legend.text = element_text(size = 14),
  title = element_text(size = 20))

Plot_S1 <- ggplot(AX) + 
  geom_point(aes(x=sec, y = CO2)) +
  xlim(XR1, XR2) +
  ylim(YR1,YR2) +
  labs(title = "02-12_19_CO2")


Plot_S1 + My_Theme


#### Mode transition times: judged from change in voltage of switch

#location 1
# 80 s added to transitions - mismathc between data logger data and analysers??

Start <-4000
# M1 inlet
T1 <- 4350
# M2 chamber
T2 <- 4550
# M3 inlet
T3 <- 4750
# M4 chamber
T4 <- 4945
# M5 inlet
T5 <- 5145
# M6 chamber
T6 <- 5365
# M7 inlet
T7 <- 5565
# M8 chamber
T8 <- 5775
# M9 inlet
T9 <- 5980
# M10 chamber
T10 <- 6170
# M11 chamber
Finish <- 6430



#Selecting the sequences (modes) between switching valves, at the moment ~ 60s after each transition
# and ~ 20 secs before

M1 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec >Start + 0, sec < T1 - 25)#I
M2 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T1 + 25, sec < T2 - 25)#C 
M3 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec >T2 + 25, sec < T3 - 25)#I
M4 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T3 + 25, sec < T4 - 25)#C
M5 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T4 + 25, sec < T5 - 25)#I
M6 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T5 + 25, sec < T6 - 25)#C
M7 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T6 + 20, sec < T7 - 25)#I
M8 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T7 + 25, sec < T8 - 25)#C
M9 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T8 + 25, sec < T9 - 25)#I
M10 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T9 + 25, sec < T10 - 25)#C
M11 <- AX %>%  select(sec,CO2, H2O) %>%
  filter(sec > T10 + 25, sec < Finish - 0)#I


#Plotting each mode - if you want????


Plot_M <- ggplot(M1) + 
  geom_point(aes(x=sec, y = CO2)) +
  labs(title = "CO2")
Plot_M


#Trying to account for the inherent trends during  the measurements
#adding Chamber modes together:
# and Inlet modes together: 

Chamber <- rbind(M2,M4,M6,M8,M10)
head(Chamber)
str(Chamber)

Inlet <- rbind(M1,M3,M5,M7,M9,M11)
head(Inlet)
str(Inlet)

#Assigning mode to each of the sequences and combining them in one dataframe (CHIN)

A <- rep('red',760)
B <- rep('blue',1164)
mode <- c(A,B)

CHIN <- rbind(M2,M4,M6,M8,M10, M1,M3,M5,M7,M9,M11)
CHIN <- cbind(CHIN,mode)

head(CHIN)

#Plotting Chamber and Inlet on the same plot with separte regressions


Plot_In <- ggplot(data = CHIN) +
  geom_point(mapping = aes(x = sec, y = CO2), colour = mode) +
  geom_smooth(data = Inlet, method = loess, se = TRUE, aes(x = sec, y = CO2), colour = 'blue') +
  labs(title = '02122019_CO2 Location 2')
Plot_In + My_Theme

Plot_In <- ggplot(data = CHIN) +
  geom_point(mapping = aes(x = sec, y = CO2), colour = mode) +
  geom_smooth(data = Inlet, method = loess, se = TRUE, aes(x = sec, y = CO2), colour = 'blue') +
  geom_smooth(data = Chamber, method = loess, se = TRUE, aes(x = sec, y = CO2), colour = 'dark red') +
  ylim(414.2,414.45) +
  labs(title = '02122019_CO2 Location 2')
Plot_In + My_Theme

#Plotting theactual curve fits

i_loess <- loess(CO2 ~ sec, data = Inlet, degree = 2)
c_loess <- loess(CO2 ~ sec, data = Chamber, degree = 2) # not used
i_loess

i_pred <- predict(i_loess, data.frame(sec = seq(4000, 6430, 1)), se = TRUE) # generates data frame of prediction

i_pred

i_pred.df <- data.frame(i_pred)
time <- seq(4000,6430,1) #combines time and predicitons
i_pred.df <- cbind(time, i_pred.df)
i_pred.df


Predicted_I <- ggplot(data = i_pred.df) +
  geom_point(mapping = aes(x = time, y = fit), size = 1, colour = 'dark blue')
Predicted_I


#Detrending the data?


D1 <- filter(i_pred.df, time > Start+0, time < T1- 25)
M1D1 <-  cbind(M1,D1)
M1DT <- mutate(M1D1, detrend = CO2 - fit)
head(M1DT)

D2 <- filter(i_pred.df, time > T1+25, time < T2- 25)
M2D2 <-  cbind(M2,D2)
M2DT <- mutate(M2D2, detrend = CO2 - fit)
head(M2DT)

D3 <- filter(i_pred.df, time > T2+25, time < T3- 25)
M3D3 <-  cbind(M3,D3)
M3DT <- mutate(M3D3, detrend = CO2 - fit)
head(M3DT)

D4 <- filter(i_pred.df, time > T3+25, time < T4- 25)
M4D4 <-  cbind(M4,D4)
M4DT <- mutate(M4D4, detrend = CO2 - fit)
head(M4DT)

D5 <- filter(i_pred.df, time > T4+25, time < T5- 25)
M5D5 <-  cbind(M5,D5)
M5DT <- mutate(M5D5, detrend = CO2 - fit)
head(M5DT)

D6 <- filter(i_pred.df, time > T5 + 25, time < T6 -25)
M6D6 <-  cbind(M6,D6)
M6DT <- mutate(M6D6, detrend = CO2 - fit)
head(M6DT)

D7 <- filter(i_pred.df, time >  T6 + 20, time < T7 -25)
M7D7 <-  cbind(M7,D7)
M7DT <- mutate(M7D7, detrend = CO2 - fit)
head(M7DT)

D8 <- filter(i_pred.df, time > T7 + 25, time < T8- 25)
M8D8 <-  cbind(M8,D8)
M8DT <- mutate(M8D8, detrend = CO2 - fit)
head(M8DT)

D9 <- filter(i_pred.df, time >  T8 + 25, time < T9- 25)
M9D9 <-  cbind(M9,D9)
M9DT <- mutate(M9D9, detrend = CO2 - fit)
head(M9DT)

D10 <- filter(i_pred.df, time > T9 + 25, time < T10- 25)
M10D10 <-  cbind(M10,D10)
M10DT <- mutate(M10D10, detrend = CO2 - fit)
head(M10DT)

D11 <- filter(i_pred.df, time > T10 + 25, time < Finish - 0)
M11D11 <-  cbind(M11,D11)
M11DT <- mutate(M11D11, detrend = CO2 - fit)
head(M11DT)


### putting the detrended data together...

Detrend_C <- rbind(M2DT,M4DT, M6DT, M8DT, M10DT)
Detrend_Cm <- add_column(Detrend_C, valve = 'chamber')

Detrend_I <- rbind(M1DT,M3DT, M5DT, M7DT, M9DT, M11DT)
Detrend_Im <- add_column(Detrend_I, valve = 'inlet')
head(Detrend_Im)

Detrend_B <- rbind(Detrend_Cm,Detrend_Im)

Plot_DT <- ggplot(Detrend_B) +
  geom_point(mapping=(aes(x = sec, y = detrend, colour = valve)))+
  labs(title = "231119 CO2 detrended_Site 2")
Plot_DT

Plot_DT + My_Theme
Plot_DT + My_Theme + scale_colour_distiller(palette = 'Spectral')

head(Detrend_B)

# trying to show as a box chart - differences between sequences:

??brewer

Box_plot_Detrend <- ggplot(data = Detrend_B, x = sec, y = detrend)+
  geom_point(mapping=aes(x = sec, y = detrend, colour = valve))+
  geom_boxplot(data = M1DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M2DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M3DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M4DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M5DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M6DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M7DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M8DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M9DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M10DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  geom_boxplot(data = M11DT,mapping = aes(x = sec, y = detrend), notch = TRUE)+
  labs(title = "02122019 Detrended CO2 S2")

Box_plot_Detrend + My_Theme + scale_colour_manual(values = c("red", "blue"))
Box_plot_Detrend + My_Theme + scale_colour_manual(values = c("black", "blue"))

?scale_color_manual

Cavg <- pull(Detrend_C, detrend)
Cx = mean(Cavg)
Csd <- sd(Cavg)
Iavg <- pull(Detrend_I, detrend)
Ix = mean(Iavg)
Isd <- sd(Iavg)
Cx
Csd
Ix
Isd

Cx-Ix

#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#Ambient concentration

CO2ambient <- pull(Detrend_I, CO2)
CO2air = mean(CO2ambient)
CO2sd <- sd(CO2ambient)
CO2air
CO2sd

H2Oambient <- pull(Detrend_I, H2O)
H2Oair = mean(H2Oambient)
H2Osd <- sd(H2Oambient)
H2Oair
H2Osd






#ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
#trying to determine whether the Chamber data lies above or below the trend in inlet CO2?
#t.test for all chamber versus inlet data:

TtestG <- t.test(Cavg, Iavg,alternative = c("greater"), mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95)
t.test(Cavg, Iavg,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95)
t.test(Cavg, Iavg,alternative = c("less"), mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95)


#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

#converting ppb to ng/l
#V = volume per mole V = nRT/P, K is Kelvin, 1 atm
K <-  273
R <- 0.082057
n <- 1
P <- 1
V <-  n*R*K/P
V
W <- 44
#CO2mass is ug/m3
#for CO2 Cx and Ix uits ppm = umol/mol
Cx_mass <- W * Cx     #ug/mol
Cx_vol <- Cx_mass/V    # ug/liter
Cx_m3 <- Cx_vol*1000    #ug/m3
Cx_m3

Ix_mass <- W * Ix #ug/mol
Ix_vol <- Ix_mass/V # ug/liter
Ix_m3 <- Ix_vol*1000 #ug/m3
Ix_m3







#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX


$CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

#Flow rate




#Add flow rates here: 
Q1 <- FR2_Avg  #SLPM
Qsd <- FR2_sd
Qcv <- 100*Qsd/Q1

#converting to (m3 min-1)
Qm3 <- Q1/1000
Qm3

#Add chamber diameter: inches
D1 <- 26+(5/16)
D2 <- 26+(0/16)
D3 <- 26+(0/16)
D4 <- 25+(14/16)
Diam <- c(D1,D2,D3,D4)
D_x <- mean(Diam)
D_sd <- sd(Diam)
D_cv <-  100*D_sd/D_x
D_x
D_sd
D_cv

#Chamber area (Converts D_x to meters) A is the area in m2.
D_xm <- (D_x * 2.54)/100
D_xm
A <- pi*(D_xm/2)*(D_xm/2)
A
pi

#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#Calculating flux:
#F = Q * (C-I) / A

#ng/m2/min

Flux <- Qm3 *(Cx_m3 - Ix_m3)/A
Flux                      #ug/m2/min

Flux_day <- Flux*24*60
Flux_day                  #ug/m2/d
Flux_umol_day <- Flux_day/44
Flux_umol_day             #umol/m2/d


#Listing the results

Date <- '02Dec19 < 3160 s'
Location <- 'Site 1'
Gas <- "CO2"
Vs <- c(Date, Location, Gas, Flux_day, Flux_umol_day, Q1, Qsd, CO2air, CO2sd, H2Oair, H2Osd, Cx, Csd, Ix, Isd,  T5, FR2_Avg, FR2_sd, Start, Finish)
list(Vs)

Parameters <- c('Date', 'Location', 'Gas', 'Flux/day (ug/m2)', 'Flux/day (umol/m2)', 'Q1 (SLPM)', 'Qsd', 'CO2air (ppm)', 'CO2sd', 'H2Oair (ppm)', 'H2Osd', 'Cx (ppm)', 'Csd', 'Ix (ppm)', 'Isd',  'T5', 'Flow', 'Flow_sd','Start', 'Finish')


Results <- data.frame('variables'= Parameters, 'Obs' = Vs)
head(Results)
View(Results)


# extracting the coefficients and calcualting growth rates etc
intercept <- round(summary(hyperfit)$coefficients)
summary(hyperfit)$coefficients
Intercept <- round(summary(hyperfit)$coefficients [1], 6)
Intercept
Intercept_SE <- round(summary(hyperfit)$coefficients [1,2], 6)
Intercept_SE
Max_U <- 1/Intercept
Max_U

Slope <- summary(hyperfit)$coefficient [2,1]
Slope
Net_U <- 1/(Intercept+Slope)
Net_U

summary()


#ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
#T-tests between detrended CO2 data for the different modes: should maybe do a matrix comparison:
  
#test 1
t.test(MD1, MD2,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 2
t.test(M4D, M6D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 3
t.test(M4D, M7D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 4
t.test(M4D, M8D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 5
t.test(M4D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 6
t.test(M4D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 7
t.test(M5D, M6D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 8
t.test(M5D, M7D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 9
t.test(M5D, M8D,alternative =c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 10
t.test(M5D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 11
t.test(M5D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 12
t.test(M6D, M7D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 13
t.test(M6D, M8D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 14
t.test(M6D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 15
t.test(M6D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 16
t.test(M7D, M8D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 17
t.test(M7D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 18
t.test(M7D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 19
t.test(M8D, M9D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 20
t.test(M8D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
#test 21
t.test(M9D, M10D,alternative = c("two.sided"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)