#!/usr/bin/Rscript
#library(minpack.lm)

# Usage:
#  data.get pek S11a,S12a,X1a,X2a 2017:1 2017:92 clean | da.export --cut-split --mode=R > test_input.dat
#  data.consolidate --source=edited pek 2017:1 2017:92 F1_S11 BsG_S11 U_S11 F1_S12 BsG_S12 Tu_S12 T_S12 U_S12 F1_X1 F1_X2 T_V13 U_V13|da.export --cut-split --mode=R >pekQ12017.dat
# to run: R --slave -f kappa.r > output_file

# Declare any inputs needed from the data
#inputNames <- c(
inputs <- c(
    "BsG0_S11",
    "BsG0_S12",
    "U0_S11",
    "Tu0_S12",
    "T0_S12",
    "T0_V13",
    "U0_V13",
    "U0_S12"
)

# Helper functions to calculate the logical values from the raw data.  These can also apply
# calibrations, etc.
SIMEXMeasurementError <- 0.3;
Tave<- function(data) {
 
   return(data[,T0_S12])
   }
TD <- function(data) {
   Tdewpt <-
  
243.04*(log(data[,U0_V13]/100)+((17.625*data[,T0_V13])/(243.04+data[,T0_V13])))/(17.625-log(data[,U0_V13]/100)-((17.625*data[,T0_V13])/(243.04+data[,T0_V13])))
   return( Tdewpt )
   }
   
wetRH <- function(data) {
   RHwet <- 100*(exp((17.625*TD(data))/(243.04+TD(data)))/exp((17.625*Tave(data))/(243.04+Tave(data))))
   return( RHwet )
}   

wetScattering <- function(data) {
    return( data[,BsG0_S12] )
}

dryRH <- function(data) {
    return( data[,U0_S11] )
}

dryScattering <- function(data) {
    return( data[,BsG0_S11] )
}


# Apply global filtering.  This is before scan selection, so this filtering applies globally
globalFilter <- function(data) {
    # Only allow non-MVC data
    validData <- rep(TRUE, length(data[,1]))
    for (variable in inputs) {
        validData <- validData & is.finite(data[,variable])
    }
    data <- data[validData,]

    # Exclude all points with the wet RH < 40
    data <- data[(wetRH(data) > 40),]   
    return( data )
    
}

# Apply scan constraint filtering.  This applies to each scan, returning true if the scan
# is acceptable for fitting.
scanFilter <- function(data) {
    # Abort if the scan has less than 14 points
    if (length(data[,1]) < 14)
        return (FALSE)
    # Abort if the scan doesn't get below 60% RH
    if (min(wetRH(data)) > 60)
        return (FALSE)
    # Abort if the scan doesn't get above 70% RH
    if (max(wetRH(data)) < 70)
        return (FALSE)
    # Abort if the scan doesn't span at least 20% RH
    if (max(wetRH(data))-min(wetRH(data)) < 20)
        return (FALSE)
    # Abort if the average dry scattering is less than 10 Mm-1
    if (mean(dryScattering(data)) < 10)
        return (FALSE)

    return (TRUE)
}

# Do the actual fit, returning the results list (including any confidence variables)
doFit <- function(data) {
    # Inputs to fitting
    fRH <- wetScattering(data) / dryScattering(data)
    RH <- wetRH(data)
    xdata <- wetRH(data)/(100-wetRH(data))
    result <- list()
   
   # fit <- nls(
       # The fit equation
       # fRH ~ a*(1-RH/100)^(-gamma),
       # The variables being fit
       # start=c(a=1,gamma=0.5),nls.control(maxiter=100, tol=1e-04)
 #  )
  #fitStat <- summary(fit)
 
    
   # fRH_fit <- function(fitstat) {
   #    fRH40_85 <- 4.0^(fitStat$coefficients[2,1])
    #   return( fRH40_85 )
 #} 
  	
  kappaFit <- lm(fRH ~ xdata, x=TRUE);  
    
    if (!is.null(SIMEXMeasurementError) && !is.na(SIMEXMeasurementError) &&
        SIMEXMeasurementError > 0) {
     library(simex);   
    
    sm <- simex(model=kappaFit, SIMEXvariable=("xdata"), 
        measurement.error=SIMEXMeasurementError);
    sdFitSlope <- sqrt(sm$variance.asymptotic[2,2]);
    sdFitIntercept <- sqrt(sm$variance.asymptotic[1,1]);
    fitSum <- summary(kappaFit);
    rSquared <- fitSum$r.squared;
   
    
} else {
    fitSum <- summary(kappaFit);
    sdFitSlope <- fitSum$coefficients[2,2];
    sdFitIntercept <- fitSum$coefficients[1,2];
   
} 
   
     kappaFRH_fit <- function(fitSum) {
        fRH_kappa <- (1+ fitSum$coefficients[2,1]*5.6667)/(1 +fitSum$coefficients[2,1]*0.6667)
        return( fRH_kappa )
      } 
    # Set the outputs columns
    result$BSG0 <- mean(dryScattering(data))
    result$minRH <- min(wetRH(data))
    result$maxRH <- max(wetRH(data))
    result$minT <- min(Tave(data))
    result$maxT <- max(Tave(data))
    result$TD <- mean(TD(data))
    result$U_S11 <- mean(dryRH(data))
    #result$gammaFRH <- fRH_fit(fitStat)
    #result$Gamma <- fitStat$coefficients[2,1]
    #result$A <- fitStat$coefficients[1,1]
    #result$Sigma_gamma <- fitStat$sigma
    result$kappaFrh <- kappaFRH_fit(fitSum)
    result$kappa <- kappaFit$coefficient[2]
    result$intercept <- kappaFit$coefficient[1]
    result$sdSlope<- sdFitSlope
    result$sdIntercept<- sdFitIntercept
    result$sigma_kappa<- fitSum$sigma
    result$rSquared <- rSquared
    return (result)
}

# Declare scan indexing, all points with the same index are part of the same scan
declareScans <- function(data) {
    return(floor(data$EPOCH / 3600) * 3600)
}
# Calculate the fractional year of times
calculateFractionalYear <- function(time) {
    yearStart <- time
    yearStart$sec <- 0;
    yearStart$min <- 0;
    yearStart$hour <- 0;
    yearStart$mday <- 1;
    yearStart$mon <- 0;
    yearEnd <- yearStart;
    yearEnd$year <- yearEnd$year + 1;
    yearLength <- as.numeric(yearEnd - yearStart, units="secs");
    return(yearStart$year + as.numeric(time - yearStart, units="secs") / yearLength + 1900);
}
# Generate the output labels for the scans that have been declared
outputLabels <- function(scans) {
    startTimes <- as.POSIXlt(as.integer(scans),
        tz="GMT", origin=ISOdatetime(1970,1,1,0,0,0,tz="GMT"))
    result <- list()
    #result$Epoch <- formatC(as.integer(scans), format="d")
    #result$FractionalYear <- formatC(calculateFractionalYear(startTimes), format="f", digits=9)
    result$DateTime <- strftime(startTimes, format="%F %T")
    #result$Year <- formatC(as.integer(startTimes$year + 1900), format="d")
    #result$DOY <- formatC((startTimes$yday + 1 + (startTimes$hour/24) +(startTimes$min/(60*24)) + (startTimes$sec/(60*60*24))), format="f", digits=5, width=9)
    return(result)
}

# Read and reformat data into the matrix.  Then declare variables with the associated names
# for ease of use.
#logicalInput <- read.csv("stdin")

#********************************
logicalInput <- read.csv(file='pek2017Q1.dat', header=TRUE, sep=",")
dataMatrix <- matrix(ncol=(length(inputs)+1), nrow=length(logicalInput$EPOCH))
dataMatrix[,1] <- declareScans(logicalInput)
inputNames <- inputs

for (variableIndex in seq(1,length(inputNames))) {
    variableName <- inputNames[variableIndex]
    dataMatrix[,(variableIndex+1)] <- logicalInput[[variableName]]
    assign(variableName, variableIndex+1)
}
inputs <- seq(2,length(inputNames)+1)
# Apply filtering
dataMatrix <- globalFilter(dataMatrix)

# Now split the matrix so we have a set of matrices with the row index being the scan ID
scanData <- lapply(split(dataMatrix[,seq(2,length(inputs)+1)], dataMatrix[,1]), matrix, ncol=length(inputs))
inputs <- seq(1,length(inputNames))
for (variableIndex in seq(1,length(inputNames))) {
    variableName <- inputNames[variableIndex]
    assign(variableName, variableIndex)
}

# Select valid scans
scanData <- scanData[unlist( lapply(scanData, scanFilter) )]

# Get the fitting results

results <- lapply(scanData, doFit)
# Construct output
labels <- outputLabels(names(results))
columns <- unique(unlist(lapply(results, names)))
outputData <- matrix(ncol=(length(labels)+length(columns)), nrow=length(results))
outputData[,seq(1,length(labels))] <- matrix(unlist(labels), ncol=length(labels))
outputData[,seq(length(labels)+1,length(labels)+length(columns))] <- matrix(unlist(lapply(results,
function(scan) {
    return( format(round(unlist(scan[columns]),3)) )
})), ncol=length(columns), byrow=TRUE)
columns <- c(names(labels), columns)
# Write it
write.table(outputData, file=stdout(), col.names=columns, row.names=FALSE, sep=", ", quote=FALSE)