## read and pre-process data
d0 <- read.csv("bart16eAllData.csv", sep = ",")
d0 <- d0[d0$SourceE == "3D",]
d0 <- droplevels(d0)
d0$ht <- as.factor(d0$ht)
d0$Cyclone <- as.factor(d0$Cyclone)
boxplot(d0$E ~d0$ht + d0$Cyclone)

## means and deviations with outliers
aggregate(E~Cyclone+ht,data=d0,FUN=mean)
aggregate(E~Cyclone+ht,data=d0,FUN=sd)

## remove outliers:
d2 <- d0[d0$E < 100 & d0$E > 80,]
d2 <- droplevels(d2)

## plot
boxplot(d2$E ~d2$ht + d2$Cyclone)
boxplot(d2$E ~d2$Cyclone + d2$ht)
pdf(file="Box2.pdf")
par(mfrow=c(1,2))
boxplot(d2$E ~d2$Cyclone, horizontal = FALSE)
boxplot(d2$E ~d2$ht, horizontal = FALSE)
dev.off()

# values for column M-P in Table 5:
aggregate(E~Cyclone+ht,data=d2,FUN=mean)
aggregate(E~Cyclone+ht,data=d2,FUN=sd)