4
$\begingroup$

I am performing a calculation in which I need to plot a generated response spectrum to a seismic tripartite graph. However, I am having a heck of a time finding a blank tripartite graph online. A low-resolution example is shown below.

tripartite graph

How might I go about generating a blank seismic tripartite graph? As far as I know, Excel does not have this functionality, since you need to be able to plot four separate axes in 2D with two of the axes at an angle. All axes are in logarithmic scale.

The unit system should be as such:

  • Spectral displacement: inches
  • Spectral velocity: in/sec
  • Spectral acceleration: % of $g$
  • Frequency: Hz

Otherwise, a link to a high-resolution graph would be acceptable.

$\endgroup$
  • $\begingroup$ Maybe you should construct it yourself in TikZ (LaTeX). $\endgroup$ – Karlo May 6 '16 at 15:12
  • $\begingroup$ @Karlo, if I knew how to, I would give it a go. My LaTeX knowledge is pretty sparse, and my TikZ is even less. $\endgroup$ – grfrazee May 6 '16 at 15:13
  • $\begingroup$ @Karlo, I opened a question on tex.SE for this as well. $\endgroup$ – grfrazee May 10 '16 at 15:57
2
$\begingroup$

You can try running the following R script which produces the figure below. tripartite graph paper

require(ggplot2)
require(data.table)

# The constant value grid lines
freqs = unique(c(seq(0.1, 0.2, by = 0.05), seq(0.2, 1.0, by = 0.1),
      seq(1.0, 2.0, by = 0.50), seq(2.0, 10.0, by = 1.0),
      seq(10.0, 20.0, by = 5.0), seq(20.0, 100.0, by = 10.0),
      seq(100.0, 200.0, by = 50.0), seq(200.0, 1000.0, by = 100.0)))
vels = unique(c(seq(1.0, 2.0, by = 0.50), seq(2.0, 10.0, by = 1.0),
      seq(10.0, 20.0, by = 5.0), seq(20.0, 100.0, by = 10.0),
      seq(100.0, 200.0, by = 50.0), seq(200.0, 1000.0, by = 100.0)))
accs = unique(c(seq(0.001, 0.002, by = 0.0005), seq(0.002, 0.01, by = 0.001),
      seq(0.01, 0.02, by = 0.005), seq(0.02, 0.1, by = 0.01),
      seq(0.1, 0.2, by = 0.05), seq(0.2, 1.0, by = 0.1),
      seq(1.0, 2.0, by = 0.50), seq(2.0, 10.0, by = 1.0),
      seq(10.0, 20.0, by = 5.0), seq(20.0, 100.0, by = 10.0),
      seq(100.0, 200.0, by = 50.0), seq(200.0, 1000.0, by = 100.0),
      seq(1000.0, 2000.0, by = 500.0), seq(2000.0, 10000.0, by = 1000.0)))
disps = unique(c(seq(0.0001, 0.0002, by = 0.00005), seq(0.0002, 0.001, by = 0.0001),
          seq(0.001, 0.002, by = 0.0005), seq(0.002, 0.01, by = 0.001),
          seq(0.01, 0.02, by = 0.005), seq(0.02, 0.1, by = 0.01),
          seq(0.1, 0.2, by = 0.05), seq(0.2, 1.0, by = 0.1),
          seq(1.0, 2.0, by = 0.50), seq(2.0, 10.0, by = 1.0),
          seq(10.0, 20.0, by = 5.0), seq(20.0, 100.0, by = 10.0),
          seq(100.0, 200.0, by = 50.0), seq(200.0, 1000.0, by = 100.0),
          seq(1000.0, 2000.0, by = 500.0), seq(2000.0, 10000.0, by = 1000.0)))

# Horizontal grid lines (const vel = in/s)
vel_grid = data.table()
freqmin = min(freqs)
freqmax = max(freqs)
for (vel in vels) {
  xx = c(freqmin, freqmax)
  yy = c(vel, vel)
  label = c(paste0("vel",vel), paste0("vel",vel)) 
  color = "vel"
  local_dt = data.table(x = xx, y = yy, label = label, color = color)
  vel_grid = rbind(vel_grid, local_dt)
}

# Vertical grid lines (freq = cycles/sec)
freq_grid = data.table()
velmin = min(vels)
velmax = max(vels)
for (freq in freqs) {
  xx = c(freq, freq)
  yy = c(velmin, velmax)
  label = c(paste0("freq",freq), paste0("freq",freq)) 
  color = "freq"
  local_dt = data.table(x = xx, y = yy, label = label, color = color)
  freq_grid = rbind(freq_grid, local_dt)
}

# -45 grid lines (constant acc = omega * vel = (2*pi*freq) radians/sec * vel in/sec  
#                              = (2*pi*freq)*vel in/sec^2 = (2*pi*freq)*vel/g  Gs)
gravity = 32.2*12.0  # in/s^2
acc_grid = data.table()
for (acc in accs) {
  for (freq in freqs) {
    f = freq
    vel = acc*gravity/(freq*2*pi)
    if (vel < velmin) {
      vel = velmin
      f = acc*gravity/(vel*2*pi)
    }
    if (vel > velmax) {
      vel = velmax
      f = acc*gravity/(vel*2*pi)
    }
    xx = c(f)
    yy = c(vel)
    label = c(paste0("acc",acc)) 
    color = "acc"
    local_dt = data.table(x = xx, y = yy, label = label, color = color)
    acc_grid = rbind(acc_grid, local_dt)
  }
}

## The const acc=1 line is the disp axis
acc_1 = acc_grid[acc_grid$label == "acc1"]
disp_axis = data.table()
for (disp in disps) {
  acc = 1.2
  omega = sqrt(acc*gravity/disp)
  freq = omega/(2*pi)
  vel = omega*disp
  xx = c(freq)
  yy = c(vel)
  label = c(paste0("disp_axis",acc)) 
  color = "disp_axis"
  local_dt = data.table(x = xx, y = yy, z = disp, label = label, color = color)
  disp_axis = rbind(disp_axis, local_dt)
}

# +45 grid lines (const disp = vel/omega = vel/(2*pi*freq) in)
disp_grid = data.table()
for (disp in disps) {
  for (freq in freqs) {
    f = freq
    vel = disp*freq*(2*pi)
    if (vel < velmin) {
      vel = velmin
      f = vel/(disp*2*pi)
    }
    if (vel > velmax) {
      vel = velmax
      f = vel/(disp*2*pi)
    }
    xx = c(f)
    yy = c(vel)
    label = c(paste0("disp",disp)) 
    color = "disp"
    local_dt = data.table(x = xx, y = yy, label = label, color = color)
    disp_grid = rbind(disp_grid, local_dt)
  }
}

## The const disp=0.1 line is the acc axis
disp_1 = disp_grid[disp_grid$label == "disp0.1"]
acc_axis = data.table()
for (acc in accs) {
  disp = 0.12
  omega = sqrt(acc*gravity/disp)
  freq = omega/(2*pi)
  vel = omega*disp
  xx = c(freq)
  yy = c(vel)
  label = c(paste0("acc_axis",acc)) 
  color = "acc_axis"
  local_dt = data.table(x = xx, y = yy, z = acc, label = label, color = color)
  acc_axis = rbind(acc_axis, local_dt)
}

# Join vertical and horizontal
grid = rbind(vel_grid, freq_grid, acc_grid, disp_grid)

### Axis label points
disp = 0.08
acc = 1000
omega = sqrt(acc*gravity/disp)
freq = omega/(2*pi)
vel = omega*disp
acc_axis_label = data.table(x = c(freq), y = c(vel), label="Spectral acceleration (g)")

disp = 200
acc = 0.8
omega = sqrt(acc*gravity/disp)
freq = omega/(2*pi)
vel = omega*disp
disp_axis_label = data.table(x = c(freq), y = c(vel), label="Spectral displacement (in)")

# Plot horizontal grid
xticks = freqs[c(TRUE, FALSE)]
yticks = vels[c(TRUE, FALSE)]

plt = ggplot() +
      geom_line(data = grid, aes(x = x, y = y, group = label, color = color),
                size=0.25) +
      geom_line(data = acc_1, aes(x = x, y = y), color = 1, size=0.5) + 
      geom_line(data = disp_1, aes(x = x, y = y), color = 1, size=0.5) + 
      xlab("Frequency (Hz)") + 
      ylab("Spectral velocity (in/sec)") + 
      scale_color_discrete(guide = FALSE) +
      scale_x_log10(limits = c(freqmin, freqmax),
                    expand = c(0,0),
                    breaks=xticks,
                    labels= c("0.1","0.2","0.4","0.6","0.8", xticks[xticks >= 1])) +
      scale_y_log10(limits = c(velmin, velmax),
                    expand = c(0,0),
                    breaks=yticks) +
      coord_fixed() +
      theme_bw() +
      theme(panel.grid.major = element_blank(),
            panel.grid.minor = element_blank()) +
      annotate("text", x = acc_axis$x, y = acc_axis$y, label = as.character(acc_axis$z), 
               size = 3, angle = -45) + 
      annotate("text", x = disp_axis$x, y = disp_axis$y, label = as.character(disp_axis$z), 
               size = 3, angle = 45)  +
      annotate("text", x = acc_axis_label$x, y = acc_axis_label$y, label = acc_axis_label$label,
               size = 4, angle = 45)  +
      annotate("text", x = disp_axis_label$x, y = disp_axis_label$y, label = disp_axis_label$label,
               size = 4, angle = -45) 
print(plt)
ggsave("seismic_tripartite_paper1.pdf")
| improve this answer | |
$\endgroup$
  • $\begingroup$ Perfect. This is about exactly what I was looking for. $\endgroup$ – grfrazee May 8 '16 at 17:22
  • $\begingroup$ @grfrazee You can then accept the answer so that the question gets removed from the unanswered list. $\endgroup$ – Biswajit Banerjee May 8 '16 at 19:59
  • 1
    $\begingroup$ I'm well aware of how the Question process works here. You can rest assured that I'll accept an answer when I'm good and ready to do so. $\endgroup$ – grfrazee May 9 '16 at 1:07
  • $\begingroup$ Then I presume my answer is not really what you're looking for. Could you explain what you are looking for, exactly? $\endgroup$ – Biswajit Banerjee May 9 '16 at 2:27
  • $\begingroup$ There is more than one way to solve a problem. Someone else might have a solution that's different than yours, and it could be beneficial to me or to someone else. Accepting yours now will detract others who might have another answer. $\endgroup$ – grfrazee May 9 '16 at 3:10
2
$\begingroup$

Tripartite Paper - Period on X-axis This is the only place on the internet that I could find that actually gave decent resolution images for tripartite plotting paper. However, it had some limits that affected use in our area: 1) the design earthquakes top out in our area at about 1g and 2) we generally are only concerned with periods between 0.01s and 10s, or 0.1 Hz to 100 Hz. This resulted in plots of the earthquake crammed down in one corner, and the 144-yr operating basis earthquake barely showed up at all. Also, we generally plot with the period, rather than the frequency, on the X-axis.

I took a weekend to learn just enough R to modify Biswajit Banerjee's excellent script to come up with a plot that would serve my purposes. I'm going to throw it up in case it's useful for anybody else. If you've not worked with R, you can install the basic R interpreter, cut and paste either Mr. Banerjee's or my code into a script, and run it to get a higher-quality PDF output.

It didn't take long to learn enough to actually plot the earthquakes using relatively easy code, but the way I did it is extremely user-unfriendly, if effective.

require(ggplot2)
require(data.table)

# The constant value grid lines
periods = unique(c(seq(0.01, 0.02, by = 0.005), seq(0.02, 0.1, by = 0.01),
      seq(0.1, 0.2, by = 0.05), seq(0.2, 1.0, by = 0.1),
      seq(1.0, 2.0, by = 0.50), seq(2.0, 10.0, by = 1.0)))
vels = unique(c(seq(0.01, 0.02, by = 0.005), seq(0.02, 0.1, by = 0.01),
    seq(0.1, 0.2, by = 0.05), seq(0.2, 1.0, by = 0.1),
    seq(1.0, 2.0, by = 0.50), seq(2.0, 10.0, by = 1.0),
      seq(10.0, 20.0, by = 5.0)))
accs = unique(c(seq(0.0001, 0.0002, by = 0.00005),
    seq(0.0002, 0.001, by = 0.0001),
    seq(0.001, 0.002, by = 0.0005), seq(0.002, 0.01, by = 0.001),
      seq(0.01, 0.02, by = 0.005), seq(0.02, 0.1, by = 0.01),
      seq(0.1, 0.2, by = 0.05), seq(0.2, 1.0, by = 0.1),
      seq(1.0, 2.0, by = 0.50), seq(2.0, 10.0, by = 1.0),
      seq(10.0, 20.0, by = 5.0), seq(20.0, 100.0, by = 10.0)))
disps = unique(c(seq(0.00001, 0.00002, by = 0.000005), seq(0.00002, 0.0001, by = 0.00001),
          seq(0.0001, 0.0002, by = 0.00005), seq(0.0002, 0.001, by = 0.0001),
          seq(0.001, 0.002, by = 0.0005), seq(0.002, 0.01, by = 0.001),
          seq(0.01, 0.02, by = 0.005), seq(0.02, 0.1, by = 0.01),
          seq(0.1, 0.2, by = 0.05), seq(0.2, 1.0, by = 0.1),
          seq(1.0, 2.0, by = 0.50), seq(2.0, 10.0, by = 1.0),
          seq(10.0, 20.0, by = 5.0)))

# Horizontal grid lines (const vel = in/s)
vel_grid = data.table()
periodmin = min(periods)
periodmax = max(periods)
for (vel in vels) {
  xx = c(periodmin, periodmax)
  yy = c(vel, vel)
  label = c(paste0("vel",vel), paste0("vel",vel)) 
  color = "vel"
  local_dt = data.table(x = xx, y = yy, label = label, color = color)
  vel_grid = rbind(vel_grid, local_dt)
}

# Vertical grid lines (freq = cycles/sec)
period_grid = data.table()
velmin = min(vels)
velmax = max(vels)
for (period in periods) {
  xx = c(period, period)
  yy = c(velmin, velmax)
  label = c(paste0("period",period), paste0("period",period)) 
  color = "period"
  local_dt = data.table(x = xx, y = yy, label = label, color = color)
  period_grid = rbind(period_grid, local_dt)
}

# 45 grid lines (constant acc = omega * vel = (2*pi/period) radians/sec * vel in/sec  
#                              = (2*pi/period)*vel in/sec^2 = (2*pi/period)*vel/g  Gs)
gravity = 32.2*12.0  # in/s^2
acc_grid = data.table()
for (acc in accs) {
  for (period in periods) {
    Ts = period
    vel = acc*gravity*period/(2*pi)
    if (vel < velmin) {
      vel = velmin
      Ts = vel*2*pi/(acc*gravity)
    }
    if (vel > velmax) {
      vel = velmax
      Ts = vel*2*pi/(acc*gravity)
    }
    xx = c(Ts)
    yy = c(vel)
    label = c(paste0("acc",acc)) 
    color = "acc"
    local_dt = data.table(x = xx, y = yy, label = label, color = color)
    acc_grid = rbind(acc_grid, local_dt)
  }
}

## The const acc=1 line is the disp axis
acc_1 = acc_grid[acc_grid$label == "acc0.01"]
disp_axis = data.table()
for (disp in disps) {
  acc = 0.012
  omega = sqrt(acc*gravity/disp)
  period = (2*pi)/omega
  vel = omega*disp
  xx = c(period)
  yy = c(vel)
  label = c(paste0("disp_axis",acc)) 
  color = "disp_axis"
  local_dt = data.table(x = xx, y = yy, z = disp, label = label, color = color)
  disp_axis = rbind(disp_axis, local_dt)
}

# -45 grid lines (const disp = vel/omega = vel/(2*pi/period) in)
disp_grid = data.table()
for (disp in disps) {
  for (period in periods) {
    Ts = period
    vel = disp*(2*pi)/period
    if (vel < velmin) {
      vel = velmin
      Ts = disp*2*pi/vel
    }
    if (vel > velmax) {
      vel = velmax
      Ts = disp*2*pi/vel
    }
    xx = c(Ts)
    yy = c(vel)
    label = c(paste0("disp",disp)) 
    color = "disp"
    local_dt = data.table(x = xx, y = yy, label = label, color = color)
    disp_grid = rbind(disp_grid, local_dt)
  }
}

## The const disp=0.1 line is the acc axis
disp_1 = disp_grid[disp_grid$label == "disp0.1"]
acc_axis = data.table()
for (acc in accs) {
  disp = 0.12
  omega = sqrt(acc*gravity/disp)
  period = 2*pi/omega
  vel = omega*disp
  xx = c(period)
  yy = c(vel)
  label = c(paste0("acc_axis",acc)) 
  color = "acc_axis"
  local_dt = data.table(x = xx, y = yy, z = acc, label = label, color = color)
  acc_axis = rbind(acc_axis, local_dt)
}

# Join vertical and horizontal
grid = rbind(vel_grid, period_grid, acc_grid, disp_grid)

### Axis label points
disp = 0.08
acc = 0.002
omega = sqrt(acc*gravity/disp)
period = 2*pi/omega
vel = omega*disp
acc_axis_label = data.table(x = c(period), y = c(vel), label="Spectral acceleration (g)")

disp = 0.006
acc = 0.008
omega = sqrt(acc*gravity/disp)
period = 2*pi/omega
vel = omega*disp
disp_axis_label = data.table(x = c(period), y = c(vel), label="Spectral displacement (in)")

# Plot horizontal grid
xticks = periods[c(TRUE, FALSE)]
yticks = vels[c(TRUE, FALSE)]


plt = ggplot() +
      geom_line(data = grid, aes(x = x, y = y, group = label, color = color),
                size=0.25) +
      geom_line(data = acc_1, aes(x = x, y = y), color = 1, size=0.5) + 
      geom_line(data = disp_1, aes(x = x, y = y), color = 1, size=0.5) + 
      xlab("Period (s)") + 
      ylab("Spectral velocity (in/sec)") + 
      scale_color_discrete(guide = FALSE) +
      scale_x_log10(limits = c(periodmin, periodmax),
                    expand = c(0,0),
                    breaks=xticks,
                    labels= c("0.01","0.02","0.04","0.06","0.08", xticks[xticks >= 0.1])) +
      scale_y_log10(limits = c(velmin, velmax),
                    expand = c(0,0),
                    breaks=yticks) +
      coord_fixed() +
      theme_bw() +
      theme(panel.grid.major = element_blank(),
            panel.grid.minor = element_blank()) +
      annotate("text", x = acc_axis$x, y = acc_axis$y, label = as.character(acc_axis$z), 
               size = 3, angle = 45) + 
      annotate("text", x = disp_axis$x, y = disp_axis$y, label = as.character(disp_axis$z), 
               size = 3, angle = -45)  +
      annotate("text", x = acc_axis_label$x, y = acc_axis_label$y, label = acc_axis_label$label,
               size = 4, angle = -45)  +
      annotate("text", x = disp_axis_label$x, y = disp_axis_label$y, label = disp_axis_label$label,
               size = 4, angle = 45)

print(plt)
ggsave("TripartitePaper.pdf",width = 10.0, height = 10.0, units = "in")
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.