# Grimm spectrometer ####
# Load packages
library('parallel')
library('ggplot2')
library('openair')
library('lubridate')
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library('plot3D')
# Load data ####
data_files <- c('20180129_parsed',
                '20180130_parsed',
                '20180131_parsed',
                '20180201_parsed',
                '20180202_parsed')
for (data_file in data_files){
  data_path <- '~/data/grimm/nelson_dust/WORKING/'
  grimm <- read.csv(paste0(data_path,data_file,'.txt'), sep=";", quote="", na.strings="NaN")
  grimm$date<-as.POSIXct(grimm$date,tz='Etc/GMT-12')
  nfield=length(grimm[1,])
  # Convert to particles/cc (originally as particles/litre)
  grimm[,3:nfield] <- grimm[,3:nfield]/1000
  grimm_dndlogdp = grimm
  Dp = as.numeric(sub("n","",names(grimm)[3:length(names(grimm))]))
  if (length(Dp)<10){
    # 107 3 channels
    dlogDp = c(1,1,1)
  } else if (length(Dp)<20){
    # Silver 108 units
    dlogDp=c(0.124939,0.096910,0.113943,0.090177,0.096910,0.204120,0.096910,0.176091,0.124939,0.096910,0.176091,0.124939,0.176091,0.124939,1)
  } else {
      # Green 109 units
      dlogDp=c(0.0492180227,0.0299632234,0.0669467896,0.057991947,0.0511525224,0.0457574906,0.0644579892,0.0494853631,0.0321846834,0.057991947,0.096910013,0.1139433523,0.0901766303,0.096910013,0.096910013,0.079181246,0.0669467896,0.057991947,0.096910013,0.1139433523,0.0621479067,0.0543576623,0.0705810743,0.096910013,0.079181246,0.0669467896,0.057991947,0.096910013,0.079181246,0.0280287236,1)
  }
  
  for (i in 1:length(Dp)-1){
    # Calculate dN/dLogDp
    grimm_dndlogdp[,i+2] = (grimm[,i+2] - grimm[,i+3])/dlogDp[i]
  }

  grimm_dndlogdp$grimmN<-grimm[,3]
  grimm_dndlogdp[grimm_dndlogdp==0] <- NA
  # Plotting summaries
  
  timePlot(mydata = grimm,
           pollutant = names(grimm)[c(3)],
           group = TRUE,
           y.relation = 'free',
           main = paste(data_file,'1 minute'),
           frame.plot=TRUE)
  
#  png(paste0(data_path,data_file,'_time_series.png'),width = 1200,height = 600)
  timePlot(mydata = grimm_dndlogdp,
           pollutant = names(grimm_dndlogdp)[c(3,5,7,10,12)],
           avg.time = '15 min',
           group = TRUE,
           y.relation = 'free',
           main = paste(data_file,'15 minute dN/dLogDp'),
           frame.plot=TRUE)
#  dev.off()
#  png(paste0(data_path,data_file,'_shades.png'),width = 1200,height = 600)
  image2D(data.matrix(log10(grimm_dndlogdp[,3:length(Dp)-1])),
          x=grimm_dndlogdp$date,
          y=Dp[1:length(Dp)-1],
          xlab='Date',
          ylab='Dp',
          log="y",
          clab = paste(data_file,"dN/dLogDp (log10 scale) 1 minute data"),
          colkey = list(dist = 0, shift = 0,
                        side = 3, length = 0.5, width = 0.5,
                        cex.clab = 1.2, col.clab = "black", line.clab = 2,
                        col.axis = "black", col.ticks = "black", cex.axis = 0.8),
          frame.plot=TRUE)
#  dev.off()
#  png(paste0(data_path,data_file,'_time_variation.png'),width = 1200,height = 1200)
  timeVariation(mydata = grimm_dndlogdp,
                pollutant = names(grimm_dndlogdp)[c(3,5,7,10,12)],
                normalise = TRUE,
                main = data_file,
                frame.plot=TRUE)
#  dev.off()
  rm(grimm)
  rm(grimm_dndlogdp)
  rm(Dp)
  rm(dlogDp)
}