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 }