source('tools/support/findQuakes.R')
source('tools/support/writeContours.R')
Rcpp::sourceCpp('tools/support/trapezoidal.cpp')

finalTime              <- 1000 # s
specialTrenchDistances <- c(0.15,0.30,0.45) # m
convergenceVelocity    <- 5e-5 # m/s

last      <- function(x) x[[length(x)]]
paste.    <- function(...) paste(..., sep='.')
pasteColon<- function(...) paste(..., sep=':')

directories <- ini::read.ini('config.ini')$directories
dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
for (basedir in c("rfpitol=100e-7")) {
  dir          <- file.path(directories[['simulation']],
                            '2d-lab-fpi-tolerance', basedir)
  h5file       <- h5::h5file(file.path(dir, 'output.h5'), 'r')
  relativeTime <- h5file['relativeTime'][]
  realTime     <- finalTime * relativeTime

  calcMask     <- relativeTime >= 0.7
  calcTime     <- realTime[calcMask]
  calcDuration <- last(calcTime) - calcTime[[1]]
  calcRange    <- range(which(calcMask))
  calcLength   <- sum(calcMask)

  ## We are interested in an enlarged time range around actual events here,
  ## (and no other quantities!) hence we pass a very low velocity here.
  quakes           <- findQuakes(1e-6 + convergenceVelocity,
                                 h5file['/frictionalBoundary/velocity'][,,],
                                 indices = calcRange[1]:calcRange[2])
  quakes$beginning <- realTime[quakes$beginningIndex]
  quakes$ending    <- realTime[quakes$endingIndex]
  quakes$duration  <- quakes$ending - quakes$beginning
  numQuakes        <- nrow(quakes)

  relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
                               quakes[[numQuakes,  'ending']]), f=0.02)
  threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])

  plotMask   <- threeQuakeTimeMask
  plotTime   <- realTime[plotMask]
  plotRange  <- range(which(plotMask))
  plotLength <- sum(plotMask)

  write(relaxedTime[[1]],
        file.path(directories[['output']],
                  paste.(pasteColon('timeframe', 'min', 'threequakes', basedir),
                         'tex')))
  write(relaxedTime[[2]],
        file.path(directories[['output']],
                  paste.(pasteColon('timeframe', 'max', 'threequakes', basedir),
                         'tex')))

  surfaceCoordinates    <- h5file['/surface/coordinates'][]
  surfaceTrenchDistance <- surfaceCoordinates[,1] / cos(atan(.27))
  perm                  <- order(surfaceTrenchDistance)

  displacement <- h5file['/surface/displacement']
  # This is the displacement imposed through the Dirichlet condition at the last node
  displacementOffset <- displacement[calcRange[1]:calcRange[2],last(perm),1]
  balancedSurfaceDisplacement <- matrix(nrow = plotLength, ncol = displacement@dim[2])
  for (k in 1:displacement@dim[2]) {
    d <- displacement[calcRange[1]:calcRange[2],k,1:2]
    d[,,1] <- d[,,1] - displacementOffset
    # We're in a tilted coordinate system
    dv <- -sin(atan(.27))*d[,,1] + cos(atan(.27))*d[,,2]
    meanVertSurfaceDisplacement  <- trapezoidal(calcTime, dv) / calcDuration
    balancedSurfaceDisplacement[,k] <- (dv - meanVertSurfaceDisplacement)[
      (plotRange[1] - calcRange[1] + 1):(plotRange[2] - calcRange[1] + 1)
    ]
  }

  ## Interpolate velocity to special points from surrounding velocities
  {
    data <- matrix(nrow=plotLength, ncol=1+length(specialTrenchDistances))
    data[,1] <- plotTime
    for (k in seq(plotLength)) {
      interpolant <- approxfun(surfaceTrenchDistance[perm],
                               balancedSurfaceDisplacement[k,perm])
      for (i in seq(specialTrenchDistances)) {
        specialTrenchDistance <- specialTrenchDistances[[i]]
        data[k,i+1] <- interpolant(specialTrenchDistance)
      }
    }
    singleDataName <- file.path(directories[['output']],
                                paste.(pasteColon('dip-single-points', basedir),
                                       'csv'))
    write.csv(data, singleDataName, row.names = FALSE)
  }
  h5::h5close(h5file)

  printlevels <- c('-40','-20','-10','-5','-2.5','0','2.5','5','10','20','40')
  levels      <- as.numeric(printlevels) * 1e-6

  ret <- contourLines(plotTime,
                      surfaceTrenchDistance[perm],
                      balancedSurfaceDisplacement[,perm],
                      levels=levels)

  for (i in seq(printlevels))
    writeContours(ret, levels[[i]],
                  file.path(directories[['output']],
                            paste.(pasteColon('2d-dip-contours', basedir,
                                              'level', printlevels[[i]]),
                                   'tex')))
}