###############################################################################
# Pathname to location of main source codes:
source("P:/Research_Tools/EMR-MRC/codes-EMR/EMR_calc_v39.1.r")
source("P:/current_development-EMR-MRC/codes-EMR-calc-plot/EMR_plot_v36.3b.r")

# Folders for input data and output:
inputlocation = "P:/current_hydrograph-sites/SH_CZO/inputs/"
outputlocation = "P:/current_hydrograph-sites/SH_CZO/outputs/"
###############################################################################

# File names for input data and for output (use "" for no output file):
inputfile = "SH-Q.2009.7mo.cd.csv" # units are d, m^3/d, mm
outputfile = "SH-Q.2009.7mo.out.80QC.csv"

# Type of hydrograph: 
# "WT" for water table, "QS" for streamflow with episode end set by the given 
# stabilization time, "QC" for streamflow with episode end set at the 
# intersection of the constant-slope hydrograph separation with the hydrograph,
# "QM" for the start set by a moving average criterion and end as for "QC"
hytype = "QC"

# Site and input information--input units, watershed area (for QS), labeling:
inpunits = c("d", "m^3/d", "mm") 
wsarea =  .079  # km^2 
probeloc = "Shale Hills"
timeper = "1/20/09-8/17/09" 

# MRC coefficient matrix: Each row represents a different part (portion of 
# range) of the MRC, with the intercept in the first column and the slope in  
# the second.
mrc_coef = matrix(nrow=2, ncol=2)
mrc_coef[1,] = c(0, -0.207971474833795)  # 2018-07-17 09:37:46   
mrc_coef[2,] = c(45.6112458589481, -0.810849050006656) # 2018-07-17 09:39:58  
  
# Specific yield (when data are with water-table levels) or other appropriate 
# measure of response sensitivity (unimportant for streamflow):
capacity = 1

# Lag time between input and response: 
lag_time = .7

# Fluctuation tolerance (the maximum dH/dt value allowable as system noise 
# rather than response to incoming water flux):
fluc_tol = 70

# Minimum amount of precip within precip window to allow inclusion as an 
# episode: Use a negative number to allow no exclusions on this basis.  
minprecip = .5 

# Parameters for moving minimum criterion: For use with QM type only. For 
# WT, QS, and QC, these parameters will be ignored and may be omitted. 
# The first parameter is the width of the window (units of time) preceding 
# an episode that will be examined for its minimum H. The second is the 
# threshold (units of H) by which H must exceed the moving minimum to indicate 
# that an episode has begun.
epst_par = c(.2, .5)

# Parameter used in determining the end of an episode:
# For QS, time allowed for stabilization of return to MRC behavior 
# For WT, as contingency, time allowed for stabilization  
# For QC or QM, slope [L^3 T^-2] of hydrograph separation line
epend_par = 1.5

# For smoothing of hydrograph-slope values, number of data points to left and 
# to right of the selected point that will be averaged (weighted in proportion 
# to proximity to that point): Use 0 for no smoothing. 
Nsmooth = 8

# Axis limits (set = for autogenerated limits):
limt = c(20, 230) 
limresp = c(0, 500)  
limprecip = c(0, 0) 
limsl = c(-100, 200) 

# End of frequently-changed parameters ++++++++++++++++++++++++++++++++++++++++

# Units for output, which can be different from input units:
outunits = vector("character")
outunits["precip"] = "mm"
outunits["prectime"] = "h" 
outunits["vol"] = "m^3" 

# Conversion factors, input to output units: 
unicon = vector("numeric") 
unicon["precip"] = 1
unicon["prectime"] = 1/24
unicon["epis_product"] = 1  # Product is typically recharge or stormflow 
unicon["volarea"] = (1/wsarea) * 1.e-3

###############################################################################

# Read data and convert labeling:
rawdata = read.csv(paste(inputlocation, inputfile, sep=""), header=TRUE)
axislabels = read.csv(paste(inputlocation, inputfile, sep=""), 
                      nrows = 1, header=FALSE)

# Calculate:
outfile_csv = paste0(outputlocation, outputfile)
if (outputfile == "") outfile_csv = NULL 

toplabel = paste(probeloc, timeper, hytype, fluc_tol, lag_time, epend_par, 
                 sep="  ") 
EMR_results = applyEMR(inputdata=rawdata, hytype=hytype, fluc_tol=fluc_tol, 
    lag_time=lag_time, epend_par=epend_par, epst_par=epst_par,  
    capacity=capacity, mrc_coef=mrc_coef, Nsmooth=Nsmooth, 
    precip_min=minprecip, outfile_csv=outfile_csv, toplabel=toplabel,  
    inputfile=inputfile, unicon=unicon, inpunits=inpunits, outunits=outunits) 

# Plot results:
plot_emr(EMR_results, plot_type="mrc", limt=limt, limresp=limresp, 
         limprecip=limprecip, limsl=limsl)
plot_emr(EMR_results, plot_type="hydro", limt=limt, limresp=limresp, 
         limprecip=limprecip, limsl=limsl)

###############################################################################