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

finalTime           <- 1000 # s
convergenceVelocity <- 5e-5 # m/s

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']

  basalCoordinates <- h5file['/frictionalBoundary/coordinates'][]
  basalTrenchDistance <- basalCoordinates[,1]
  perm                <- order(basalTrenchDistance)
  sortedBasalTrenchDistance <- basalTrenchDistance[perm]

  {
    ## 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, velocityProxy,
                         indices = 1:dim(velocityProxy)[1], 1)
    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)
    plotMask   <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
    plotIndices<- which(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')))

    printlevels <- c('1000','100','10','1')
    levels      <- 1e-6 * as.numeric(printlevels) + convergenceVelocity
    ret <- contourLines(realTime[plotIndices],
                        sortedBasalTrenchDistance,
                        abs(velocityProxy[plotIndices,perm,1][,,1]),
                        levels = levels)

    for (i in seq(printlevels))
      writeContours(ret, levels[[i]],
                    file.path(directories[['output']],
                              paste.(pasteColon('2d-velocity-contours',
                                                'threequakes', basedir, 'level',
                                                printlevels[[i]]), 'tex')))
  }
  {
    ## We are interested in an enlarged time range around actual events here,
    ## (and no other quantities!) hence we pass a rather low velocity here.
    quakes <- findQuakes(300e-6 + convergenceVelocity, velocityProxy,
                         indices = 1:dim(velocityProxy)[1], 1)
    quakes$beginning <- realTime[quakes$beginningIndex]
    quakes$ending    <- realTime[quakes$endingIndex]
    quakes$duration  <- quakes$ending - quakes$beginning
    numQuakes        <- nrow(quakes)
    quake            <- quakes[numQuakes,]
    relaxedTime      <-
      c(quake[['beginning']] - 0.9*(quake[['ending']] - quake[['beginning']]),
        quake[['ending']]    + 0.1*(quake[['ending']] - quake[['beginning']]))
    plotMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
    plotIndices<- which(plotMask)

    printlevels <- c('3000','1000','300','100','30','10','3','1')
    levels      <- 1e-6 * as.numeric(printlevels) + convergenceVelocity
    ret <- contourLines(realTime[plotIndices],
                        sortedBasalTrenchDistance,
                        abs(velocityProxy[plotIndices,perm,1][,,1]),
                        levels = levels)

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