Skip to content

Instantly share code, notes, and snippets.

@momeara
Created February 18, 2015 15:21
Show Gist options
  • Save momeara/c0ee5e75dcee6af15e73 to your computer and use it in GitHub Desktop.
Save momeara/c0ee5e75dcee6af15e73 to your computer and use it in GitHub Desktop.
sp3 BAH/CHI discontinuity at BAH=0
library(dplyr)
library(ggplot2)
bah_chi_compute_energy_sp3_talaris2013 <- function(
bah,
chi,
acc_don_scale
) {
acc_don_scale * 0.25 * (1 + cos(chi)) / 2
}
bah_chi_compute_energy_sp3_talaris2014 <- function(
bah,
chi,
acc_don_scale,
pinch
) {
acc_don_scale * (bah/(pinch + bah)) * 0.25 * (1 + cos(chi)) / 2
}
d <- expand.grid(
bah=seq(0, 2*pi, length.out=50), # 0 ~ linear
chi=seq(0,2*pi, length.out=50))
d <- rbind(
d %>% mutate(fn="talaris2013", energy = bah_chi_compute_energy_sp3_talaris2013(bah, chi, 1)),
d %>% mutate(fn="talaris2014 .1", energy = bah_chi_compute_energy_sp3_talaris2014(bah, chi, 1, .1)),
d %>% mutate(fn="talaris2014 .2", energy = bah_chi_compute_energy_sp3_talaris2014(bah, chi, 1, .2)),
d %>% mutate(fn="talaris2014 .25", energy = bah_chi_compute_energy_sp3_talaris2014(bah, chi, 1, .25)),
d %>% mutate(fn="talaris2014 .3", energy = bah_chi_compute_energy_sp3_talaris2014(bah, chi, 1, .3)),
d %>% mutate(fn="talaris2014 .4", energy = bah_chi_compute_energy_sp3_talaris2014(bah, chi, 1, .4)),
d %>% mutate(fn="talaris2014 .5", energy = bah_chi_compute_energy_sp3_talaris2014(bah, chi, 1, .5)),
d %>% mutate(fn="talaris2014 .8", energy = bah_chi_compute_energy_sp3_talaris2014(bah, chi, 1, .8)),
d %>% mutate(fn="talaris2014 1", energy = bah_chi_compute_energy_sp3_talaris2014(bah, chi, 1, 1)))
p <- ggplot( data=d) + theme_bw() +
geom_tile(aes(x=bah, y=chi, fill=energy)) +
scale_fill_gradient(low="blue", high="yellow") +
stat_contour(aes(x=bah, y=chi, z=energy), bins=15) +
scale_x_continuous("bah", breaks=c(0, pi/2, pi, 3*pi/2, 2*pi), labels=c("0", "pi/2", "pi", "3pi/2", "2pi")) +
scale_y_continuous("chi", breaks=c(0, pi/2, pi, 3*pi/2, 2*pi), labels=c("0", "pi/2", "pi", "3pi/2", "2pi")) +
facet_wrap(~fn) +
ggtitle("sp3 BAH/CHI potential discontinuity")
p
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment