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

finalTime           <- 1000 # s
convergenceVelocity <- 5e-5 # m/s
printlevels <- c('1','2','3','5',
                 '10','20','30','50',
                 '100','200','300','500',
                 '1000')
criticalVelocities  <- 1e-6*as.numeric(printlevels) + 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("rtol=1e-5_diam=1e-2")) {
  dir          <- file.path(directories[['simulation']],
                            '3d-lab', basedir)
  h5file       <- h5::h5file(file.path(dir, 'output.h5'), 'r')
  relativeTime <- h5file['relativeTime'][]
  realTime     <- finalTime * relativeTime

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

  ## 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(200e-6 + convergenceVelocity, velocityProxy,
                       ## Note: We only look at the last 1000 timesteps here because
                       ## we're only interested in the last few events. This
                       ## dramatically reduces RAM usage and runtime.
                       indices = seq(dim(velocityProxy)[1]-1000,
                                     dim(velocityProxy)[1]), c(1,3))
  quake <- quakes[nrow(quakes)-1,] # Q: why did we not need the -1 in julia?
  
  weakPatchGridVelocityProxy <- h5file["/weakPatchGrid/velocity"]

  stepSize <- 30 ## note: should/might differ by time/spatial resolution
  ts <- seq(quake$beginningIndex, quake$endingIndex, by=stepSize)
  dd <- data.frame(timeSteps = ts,
                   times = realTime[ts],
                   timeOffsets = realTime[ts] - realTime[quake$beginningIndex])

  fname = pasteColon('3d-velocity-contours', basedir)
  write.csv(dd, row.names = FALSE, quote = FALSE,
            file = file.path(directories[['output']],
                             paste.(pasteColon(fname, 'times'), 'csv')))

  weakPatchGridXCoordinates <- h5file["/weakPatchGrid/xCoordinates"][]
  weakPatchGridZCoordinates <- h5file["/weakPatchGrid/zCoordinates"][]
  
  for (t in ts) {
    m <- weakPatchGridVelocityProxy[t,,,]
    s <- sqrt(m[,,,1]^2 + m[,,,3]^2)
    ret <- contourLines(weakPatchGridXCoordinates, weakPatchGridZCoordinates, s,
                        level=criticalVelocities)

    for (i in seq(printlevels))
      writeContours(ret, criticalVelocities[[i]],
                    file.path(directories[['output']],
                              paste.(pasteColon(fname,
                                                'step', t,
                                                'level', printlevels[[i]]),
                                     'tex')))
  }
  h5::h5close(h5file)
}