require(graphics)
setwd('C://Users//dweise//Documents//SERDP//RC2640//Wind tunnel')
rm(list=ls(all=TRUE))
library(psych);library(compositions);library(MASS);library(lme4);library(car);library(DescTools)
library(zCompositions);library(openair);library(dplyr)
library(ggplot2);library(emmeans);library(robCompositions)
library(gplots);library(gdata);library(lattice);library(gtools);library(Formula);library(pROC);
library(magrittr);library(lmerTest);library(EnvStats);library(gginnards);library(grid);library(chron)
hm <- function(x) 1/mean(1/x,na.rm=TRUE)
sdhm <- function(x) sqrt((mean(1/x,na.rm=TRUE))^(-4)*var(1/x,na.rm=TRUE)/length(x[!is.na(x)]))

wt_fuel_temps<-read.csv(file="WindTunnelFLIRTemperatureTime.csv")
#
# convert to absolute temperature
#
wt_fuel_temps$TempK<-wt_fuel_temps$TempC+273.15
wt_temps<-subset(subset(subset(wt_fuel_temps,Fire>7),Fire!=17),Fire!=84)
DF<-wt_temps
DF$ctime<-chron(dates.=DF$Date,times.=DF$Time)
DF2<-subset(DF,select=c(Fire,ctime))
elaps <- by(DF2,DF2$Fire,function(x) apply(x,2,function(x) x-min(x)))
elaps<- do.call(rbind,elaps)
DF$elapsed<-elaps[,2]*24*3600  # convert elapsed time from days to seconds
temp_rate<- DF %>%
  group_by(Fire) %>%
  summarise(trate=(TempC-lag(TempC,1))/(elapsed-lag(elapsed,1)),.groups="keep") %>%
  filter(trate > 0 & !is.na(trate)) %>% 
  summarise(mxtrate=max(trate))
windows(9,9)
FireLab<-function(string){
  string<-paste("WT",string)
  string
}
ggplot(DF[DF$Fire<34,],aes(x=elapsed,y=TempK))+
  geom_point(color="red")+
  geom_line()+
  facet_wrap(~Fire,labeller=labeller(Fire=FireLab))+
  theme_bw()+
  labs(x='Elapsed time (s)',y="Fuel temperature (K)")+
  theme(strip.text.x=element_text(size=10,face="bold.italic"),
        axis.text.x=element_text(size=12,angle=90,hjust=1,vjust=0.5,face="bold",color="black"),
        axis.text.y=element_text(size=12,vjust=0.5,face="bold",color="black"),
        axis.title=element_text(size=15))

savePlot(filename = "Time_IRTemp_8_33",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

ggplot(DF[DF$Fire>33 & DF$Fire<59,],aes(x=elapsed,y=TempK))+
  geom_point(color="red")+
  geom_line()+
  facet_wrap(~Fire,labeller=labeller(Fire=FireLab))+
  theme_bw()+
  labs(x='Elapsed time (s)',y="Fuel temperature (K)")+
  theme(strip.text.x=element_text(size=10,face="bold.italic"),
        axis.text.x=element_text(size=12,angle=90,hjust=1,vjust=0.5,face="bold",color="black"),
        axis.text.y=element_text(size=12,vjust=0.5,face="bold",color="black"),
        axis.title=element_text(size=15))

savePlot(filename = "Time_IRTemp_34_58",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

ggplot(DF[DF$Fire>58 & DF$Fire<85,],aes(x=elapsed,y=TempK))+
  geom_point(color="red")+
  geom_line()+
  facet_wrap(~Fire,labeller=labeller(Fire=FireLab))+
  theme_bw()+
  labs(x='Elapsed time (s)',y="Fuel temperature (K)")+
  theme(strip.text.x=element_text(size=10,face="bold.italic"),
        axis.text.x=element_text(size=12,angle=90,hjust=1,vjust=0.5,face="bold",color="black"),
        axis.text.y=element_text(size=12,vjust=0.5,face="bold",color="black"),
        axis.title=element_text(size=15))

savePlot(filename = "Time_IRTemp_59_83",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)

ggplot(DF[DF$Fire>83,],aes(x=elapsed,y=TempK))+
  geom_point(color="red")+
  geom_line()+
  facet_wrap(~Fire,labeller=labeller(Fire=FireLab))+
  theme_bw()+
  labs(x='Elapsed time (s)',y="Fuel temperature (K)")+
  theme(strip.text.x=element_text(size=10,face="bold.italic"),
        axis.text.x=element_text(size=12,angle=90,hjust=1,vjust=0.5,face="bold",color="black"),
        axis.text.y=element_text(size=12,vjust=0.5,face="bold",color="black"),
        axis.title=element_text(size=15))

savePlot(filename = "Time_IRTemp_85_97",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)#
#
# Joe O'Brien IR temperatures
#
load("~/SERDP/RC2640/FJ Leaf Temperatures/.RData")

FJ<-function(string){
  string<-paste("Ft. Jackson",string)
  string
}

ggplot(leaftemp[leaftemp$reltime<=2000,],aes(reltime,TempC+273.15,col=leafplot,label=fire))+
#  ggplot(FJ_temps,aes(x=elapsed,y=TempK,color=Locat))+
#  geom_point()+
  geom_line()+
#  geom_smooth(method="gam")+
  facet_wrap(~fire,labeller=labeller(fire=FJ))+
  scale_color_brewer(
    name="Leaf Plot",
    aesthetics="col",
    type="qual",
    guide="legend",
    palette="Dark2") +
  theme_bw()+
  labs(x='Relative time (s)',y="Fuel temperature (K)")+
  theme(strip.text.x=element_text(size=15,face="bold.italic"),
        axis.text.x=element_text(size=15,angle=90,hjust=1,vjust=0.5,face="bold",color="black"),
        axis.text.y=element_text(size=15,vjust=0.5,face="bold",color="black"),
        axis.title=element_text(size=20),
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=20,face="bold"),
        legend.position=c(0.8,0.2),
        legend.key = element_rect(fill = "white", colour = "black"),
        legend.key.size=unit(20,"bigpts"))


savePlot(filename = "Time_IRTemp_FJ",
         type = "jpg",
         device = dev.cur(),
         restoreConsole = TRUE)#

maxTs<-aggregate(TempK~Fire+Locat,data=FJ_temps,max)
maxTs$TempC<-maxTs$TempK-273.15

