velocity-potentials.R (4450B)
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 if (!('contour_lwd' %in% names(args))) { 8 args$contour_lwd = 1 9 } 10 if (!('zeta_lwd' %in% names(args))) { 11 args$zeta_lwd = 4 12 } 13 if (!('sky_col' %in% names(args))) { 14 args$sky_col = 'white' 15 } 16 if (!('axis_args' %in% names(args))) { 17 args$axis_args = list() 18 } 19 if (!('box_args' %in% names(args))) { 20 args$box_args = list() 21 } 22 if (!('x_max' %in% names(args))) { 23 args$x_max = 0 24 } 25 if (!('z_min' %in% names(args))) { 26 args$z_min = -8 27 } 28 if (!('points_args' %in% names(args))) { 29 args$points_args = list(pch=21,col='black', bg='white') 30 } 31 if (!('title_args' %in% names(args))) { 32 args$title_args = list(xlab="x", ylab="z") 33 } 34 phi <- read.csv(file.path(dir, 'phi.csv')) 35 left_top_x <- 0.1 36 right_top_x <- max(phi$x) 37 # slice time and Y ranges through the center 38 slice_t <- arma.middle_element(unique(phi$t)) 39 slice_y <- arma.middle_element(unique(phi$y)) 40 print(paste('Middle elements (TY) = ', slice_t, slice_y)) 41 phi_slice <- phi[phi$t == slice_t & phi$y == slice_y & phi$x >= left_top_x & phi$z >= args$z_min,] 42 x <- unique(phi_slice$x) 43 z <- unique(phi_slice$z) 44 left_top_z <- max(phi_slice$z) 45 right_top_z <- left_top_z 46 print(paste('Velocity field size (XZ) = ', length(x), length(z))) 47 48 # convert data frame to matrix 49 seq_x <- seq_along(x) 50 seq_z <- seq_along(z) 51 indices <- data.matrix(expand.grid(seq_x, seq_z)) 52 u <- with(phi_slice, { 53 out <- matrix(nrow=length(seq_x), ncol=length(seq_z)) 54 out[indices] <- phi 55 out 56 }) 57 58 # get wave profile 59 zeta <- read.csv(file.path(dir, 'zeta.csv')) 60 zeta_slice <- zeta[zeta$t == slice_t & zeta$y == slice_y & zeta$x >= left_top_x,] 61 62 plot.new() 63 plot.window(xlim=range(c(x,args$x_max)),ylim=range(z),asp=1) 64 do.call(axis, c(list(side=1), args$axis_args)) 65 do.call(axis, c(list(side=2), args$axis_args)) 66 67 .filled.contour( 68 x, z, u, 69 levels=args$levels, 70 col=args$col 71 ) 72 73 74 contour( 75 x, z, u, 76 levels=args$levels, 77 asp=1, 78 drawlabels=TRUE, 79 add=TRUE, 80 lwd=args$contour_lwd 81 ) 82 83 top_area_x <- c(left_top_x*0.99, zeta_slice$x, right_top_x*1.01) 84 top_area_z <- c(left_top_z*1.10, zeta_slice$z, right_top_z*1.10) 85 polygon( 86 top_area_x, 87 top_area_z, 88 lwd=args$zeta_lwd, 89 border=args$sky_col, 90 col=args$sky_col 91 ) 92 lines(zeta_slice$x, zeta_slice$z, lwd=args$zeta_lwd) 93 94 # plot point having larger values than in the "compare_to" test 95 if ("compare_to" %in% names(args)) { 96 dir0 <- args$compare_to 97 phi0 <- read.csv(file.path(dir0, 'phi.csv')) 98 phi0_slice <- phi0[phi0$t == slice_t & phi0$y == slice_y & phi0$x >= left_top_x & phi0$z >= args$z_min,] 99 phi0_range <- range(phi0_slice[phi0_slice$z <= zeta_slice$z, "phi"]) 100 large_phi <- phi_slice[phi_slice$z <= zeta_slice$z & (phi_slice$phi < phi0_range[[1]] | phi_slice$phi > phi0_range[[2]]),] 101 do.call(points, c(list(x=large_phi$x, y=large_phi$z), args$points_args)) 102 } 103 104 do.call(title, args$title_args) 105 } 106 107 arma.plot_velocity_potential_field_legend <- function (...) { 108 args <- list(...) 109 levels <- args$levels 110 col <- args$col 111 plot.new() 112 plot.window( 113 xlim = c(0, 1), 114 ylim = range(levels), 115 xaxs = "i", 116 yaxs = "i" 117 ) 118 rect(0, levels[-length(levels)], 1, levels[-1L], col = col) 119 axis(4) 120 box() 121 } 122 123 arma.plot_velocity <- function(file1, file2, ...) { 124 args <- list(...) 125 data1 <- read.table(file1) 126 data2 <- read.table(file2) 127 ylim <- range(c(data1, data2)) 128 if ("ylim" %in% names(args)) { 129 ylim <- args$ylim 130 } 131 if (!("legend_x" %in% names(args))) { 132 args$legend_x <- "topleft" 133 } 134 if (!("legend_cex" %in% names(args))) { 135 args$legend_cex <- par("cex") 136 } 137 if (!('axis_args' %in% names(args))) { 138 args$axis_args = list() 139 } 140 if (!('title_args' %in% names(args))) { 141 args$title_args = list(xlab="x",ylab="u(x)") 142 } 143 plot.new() 144 plot.window( 145 xlim = c(0, nrow(data1)-1), 146 ylim = ylim 147 ) 148 lines(data1, lty=args$linetypes[[1]]) 149 lines(data2, lty=args$linetypes[[2]]) 150 do.call(axis, c(list(side=1), args$axis_args)) 151 do.call(axis, c(list(side=2), args$axis_args)) 152 box() 153 do.call(title, args$title_args) 154 legend( 155 args$legend_x, 156 c( 157 as.expression(bquote(u[1](x))), 158 as.expression(bquote(u[2](x))) 159 ), 160 lty = paste(args$linetypes), 161 cex = args$legend_cex 162 ) 163 }