source('tools/support/findQuakes.R')

finalTime              <- 1000 # s
convergenceVelocity    <- 5e-5 # m/s
criticalVelocity       <- 1000e-6 + convergenceVelocity

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

  velocityProxy<- h5file['/frictionalBoundary/velocity']

  quakes <- findQuakes(criticalVelocity, velocityProxy,
                       indices = 1:dim(velocityProxy)[1], 1)

  basalCoordinates <- h5file['/frictionalBoundary/coordinates'][]

  maximumVelocities<- maxVelocity(velocityProxy[,,1])
  displacement     <- h5file['/frictionalBoundary/displacement']
  numVertices      <- dim(displacement)[[2]]
  for (quakeID in seq(nrow(quakes))) {
    qb <- quakes[[quakeID,'beginningIndex']]
    qe <- quakes[[quakeID,'endingIndex']]
    localSlippingTimes <- abs(velocityProxy[qb:qe,,1][,,1]) > criticalVelocity

    qd   <- displacement[qb:qe,,1]
    slip <- vector(mode='numeric', length=numVertices)
    for (i in seq(numVertices)) {
      if (any(localSlippingTimes[,i])) {
        begs <- positiveStarts(localSlippingTimes[,i])
        ends <- negativeStarts(localSlippingTimes[,i])
        slip[i]<- sum(qd[ends,i,] - qd[begs,i,])
      }
    }
    quakes[quakeID,'peakSlip']     <- max(abs(slip))
    quakes[quakeID,'peakSlipRate'] <- max(maximumVelocities[qb:qe])

    maxRuptureWidth <- 0
    for (ts in seq(dim(localSlippingTimes)[[1]])) {
      st <- localSlippingTimes[ts,]
      if (!any(st))
        next;

      slippingCoordinates <- basalCoordinates[st,1] # x-coordinate only
      maxRuptureWidth <- max(maxRuptureWidth, diff(range(slippingCoordinates)))
    }
    quakes[quakeID,'ruptureWidth'] <- maxRuptureWidth
  }
  h5::h5close(h5file)

  quakes$beginning <- realTime[quakes$beginningIndex]
  quakes$ending    <- realTime[quakes$endingIndex]

  write.csv(quakes[c('beginning','ending',
                     'peakSlip','peakSlipRate',
                     'ruptureWidth')],
            row.names = FALSE, quote = FALSE,
            file = file.path(directories[['output']],
                             paste.(pasteColon('events', basedir), 'csv')))
}