comparison:lab-sim.R

paste.    <- function(...) paste(..., sep='.')
pasteColon<- function(...) paste(..., sep=':')

basedir <- 'rfpitol=100e-7'

directories <- ini::read.ini('config.ini')$directories
labdata <- within(
  read.table(unz(file.path(directories[['experiment']],
                           'B_Scale-model-earthquake-data.zip'),
                 'scaleEQs_12-31_catalogue.txt'),
             sep='\t', header=FALSE, comment.char='%',
             col.names=c('frame',
                         'meanSlipNature', 'meanSlipLab',
                         'peakSlipNature', 'peakSlipLab',
                         'ruptureWidthNature', 'ruptureWidthLab',
                         'ignored', 'ignored', 'ignored', 'ignored', 'ignored')),
  {
    recurrenceLab <- c(NA, diff(frame))/10 # (10Hz camera: 10 frames ~ 1s)
    meanSlipLab   <- meanSlipLab / 1e6     # micro m -> m
    peakSlipLab   <- peakSlipLab / 1e6})   # micro m -> m
## NB: We skip the first two quakes here, for compatibility with my old code
##     Maybe skipping the first was intentional so that each quake would have
##     a recurrence time. But skipping the second was certainly an accident
labdata <- labdata[3:nrow(labdata),]

simdata <- read.csv(file.path(directories[['output']],
                              paste.(pasteColon('events', basedir), 'csv')))

report <- function(lab, sim, quantity) {
  write.table(lab, file.path(directories[['output']],
                             paste.(pasteColon('boxplot-data', 'lab',
                                               quantity), 'tex')),
              row.names=FALSE, col.names=FALSE)
  write.table(sim, file.path(directories[['output']],
                             paste.(pasteColon('boxplot-data', 'simulation',
                                               basedir, quantity), 'tex')),
              row.names=FALSE, col.names=FALSE)
}
report(labdata$recurrenceLab, diff(simdata$beginning), 'recurrence')
report(labdata$ruptureWidthLab, simdata$ruptureWidth, 'ruptureWidth')
# meters -> millimeters
report(labdata$peakSlipLab * 1000, simdata$peakSlip*1000, 'peakSlip')