waves-16-arma

Simulation of Standing and Propagating Sea Waves with Three-Dimensional ARMA Model
git clone https://git.igankevich.com/waves-16-arma.git
Log | Files | Refs

velocity-potentials.R (2622B)


      1 arma.middle_element <- function (x) {
      2   x[round(length(x) + 0.5)/2]
      3 }
      4 
      5 arma.plot_velocity_potential_field <- function (dir, ...) {
      6   args <- list(...)
      7   phi <- read.csv(file.path(dir, 'phi.csv'))
      8   left_top_x <- 0.1
      9   right_top_x <- max(phi$x)
     10   # slice time and Y ranges through the center
     11   slice_t <- arma.middle_element(unique(phi$t))
     12   slice_y <- arma.middle_element(unique(phi$y))
     13   print(paste('Middle elements (TY) = ', slice_t, slice_y))
     14   phi_slice <- phi[phi$t == slice_t & phi$y == slice_y & phi$x >= left_top_x & phi$z >= -8,]
     15   x <- unique(phi_slice$x)
     16   z <- unique(phi_slice$z)
     17   left_top_z <- max(phi_slice$z)
     18   right_top_z <- left_top_z
     19   print(paste('Velocity field size (XZ) = ', length(x), length(z)))
     20 
     21   # convert data frame to matrix
     22   seq_x <- seq_along(x)
     23   seq_z <- seq_along(z)
     24   indices <- data.matrix(expand.grid(seq_x, seq_z))
     25   u <- with(phi_slice, {
     26     out <- matrix(nrow=length(seq_x), ncol=length(seq_z))
     27     out[indices] <- phi
     28     out
     29   })
     30 
     31   # get wave profile
     32   zeta <- read.csv(file.path(dir, 'zeta.csv'))
     33   zeta_slice <- zeta[zeta$t == slice_t & zeta$y == slice_y & zeta$x >= left_top_x,]
     34 
     35   plot.new()
     36   plot.window(xlim=range(x),ylim=range(z),asp=1)
     37   axis(1); axis(2); box()
     38 
     39   .filled.contour(
     40     x, z, u,
     41     levels=args$levels,
     42     col=args$col
     43   )
     44 
     45 
     46   contour(
     47     x, z, u,
     48     levels=args$levels,
     49     asp=1,
     50     drawlabels=TRUE,
     51     add=TRUE
     52   )
     53 
     54   top_area_x <- c(left_top_x*0.90, left_top_x*0.90, zeta_slice$x, right_top_x*1.05)
     55   top_area_z <- c(left_top_z*1.20, zeta_slice$z[[1]], zeta_slice$z, right_top_z*1.20)
     56   polygon(top_area_x, top_area_z, lwd=4, border=NA, col='white')
     57   lines(zeta_slice$x, zeta_slice$z, lwd=4)
     58   box()
     59   title(xlab="x", ylab="z")
     60 }
     61 
     62 arma.plot_velocity_potential_field_legend <- function (...) {
     63   args <- list(...)
     64   levels <- args$levels
     65   col <- args$col
     66   plot.new()
     67   plot.window(
     68     xlim = c(0, 1),
     69     ylim = range(levels),
     70     xaxs = "i",
     71     yaxs = "i"
     72   )
     73   rect(0, levels[-length(levels)], 1, levels[-1L], col = col)
     74   axis(4)
     75   box()
     76 }
     77 
     78 arma.plot_velocity <- function(file1, file2, ...) {
     79   args <- list(...)
     80   data1 <- read.table(file1)
     81   data2 <- read.table(file2)
     82   ylim <- range(c(data1, data2))
     83   if ("ylim" %in% names(args)) {
     84     ylim <- args$ylim
     85   }
     86   plot.new()
     87   plot.window(
     88     xlim = c(0, nrow(data1)-1),
     89     ylim = ylim
     90   )
     91   lines(data1, lty=args$linetypes[[1]])
     92   lines(data2, lty=args$linetypes[[2]])
     93   axis(1)
     94   axis(2)
     95   box()
     96   title(xlab="x",ylab="u(x)")
     97   legend(
     98     "topleft",
     99     c(
    100       as.expression(bquote(u[1](x))),
    101       as.expression(bquote(u[2](x)))
    102     ),
    103     lty = paste(args$linetypes)
    104   )
    105 }