
args = commandArgs(trailingOnly=TRUE)
if (length(args) == 0) {
  cat("Generate a contour plot based on RMSD\n")
  cat("Usage:\n")
  cat("       Rscript contour_plot_afmfit.R rmsd_file cutoff_val \n")
  cat("       rmsd_file     file containing the RMSD\n") 
  cat("       cutoff_val    value above which the contour plot is saturated in color \n")
  cat("       A pdf file called contour_cutoff_val is generated \n")
  cat("       The rmsd_file is topPDB_StckHbRepul.txt generated by the script get_gbStats_scan2D.sh \n")
  quit()
}

library(dplyr)
# install.packages("readxl")
library("readxl")

file.res = args[1]
threshold = as.integer(args[2]) 
# file.res = "/Users/bhandariyr/Documents/projects/afm/scripts/lm/3Dscan/top_statsAfm_Pow_8_k_146_100FrameMod/topPDB_StckHbRepul.txt"
df = read.table(file.res, header = T, sep = "\t")
# dim(df); head(df)
## read excel files
# df = read_excel("/Users/bhandariyr/Documents/projects/afm/scripts/lm/3Dscan/top_statsAfm_Pow_8_k_146_100Frame/mod_100Frame.xlsx")
dim(df)

xy = sapply(strsplit(df$kappa, "_"), function(y){return(y[1])})
df$xy = as.numeric(substring(xy, 3)) ## remove first two characters

z = sapply(strsplit(df$kappa, "_"), function(y){return(y[2])})
df$z = as.numeric(substring(z, 2)) ## remove the first character
df$en = as.numeric(df$rmsd)

# df.sub = arrange(df.sub, x, y)
# library(reshape2)
# dcast(df.sub, xy ~ z, value.var="rmsd",fill=0)

# threshold = 7
df.sub = data.frame(x = df$xy, y = df$z, z=df$en)
plot(df.sub$z)
df.sub$z.mod = ifelse(df.sub$z > threshold, threshold, df.sub$z)

## Convert the dataframe into a 2D matrix for Contour Plot with EXCEL ####

# df.sub = arrange(df.sub, x, y)
# library(tidyr)
# df.fmt = pivot_wider(df.sub, names_from = y, values_from = z)
# write.table(df.fmt, file = "/Users/bhandariyr/Documents/projects/afm/scripts/lm/3Dscan/new/random1.tsv",
#             sep = "\t", row.names = F, col.names = T)


library(akima)
pdf(paste("contour_", threshold, ".pdf"))
fld <- with(df.sub, interp(x = x, y = y, z = z.mod, linear = T))
# fld <- with(df.sub, interp(x = x, y = y, z = z, linear = T))
filled.contour(x = fld$x,
               y = fld$y,
               z = fld$z,
               color.palette =
                 colorRampPalette(c("blue", "magenta")),
               xlab = "xy",
               ylab = "z",
               main = "En",
               key.title = title(main = "En Contour", cex.main = 1))
dev.off()

