M 13
Ralph Rogge (RRO)
- August 2015
library(astroxo)
library(data.table)
library(ggplot2)
library(dplyr)
library(mixtools)Position in the Milky May
Read globular cluster raw data from input file
raw <- read.csv("VII_202.tsv", sep=";", as.is=T, strip.white=T, na.strings="", comment.char="#",
col.names=c("name",
"ra2000", "de2000",
"l", "b", "rsun", "rgc",
"ebv",
"vt", "mvt",
"feh",
"vr", "e_vr",
"rc", "rh",
"log.tc", "log.th"))
raw$l <- as.numeric(raw$l)
raw$b <- as.numeric(raw$b)
raw$rsun <- as.numeric(raw$rsun)
raw$rgc <- as.numeric(raw$rgc)
raw$ebv <- as.numeric(raw$ebv)
raw$vt <- as.numeric(raw$vt)
raw$mvt <- as.numeric(raw$mvt)
raw$feh <- as.numeric(raw$feh)
raw$vr <- as.numeric(raw$vr)
raw$e_vr <- as.numeric(raw$e_vr)
raw$rc <- as.numeric(raw$rc)
raw$rh <- as.numeric(raw$rh)
raw$log.tc <- as.numeric(raw$log.tc)
raw$log.th <- as.numeric(raw$log.th)Use some more or less clever trick to remove textual annotations from raw data ;) and convert distance from kpc to kly
gc <- raw %>%
filter(!is.na(raw$l)) %>%
dplyr::select(name, ra2000, de2000, l, b, rsun, rgc, mvt, feh, rc) %>%
mutate(rsun=pc2ly(rsun), rgc=pc2ly(rgc))
print(data.table(gc))## name ra2000 de2000 l b rsun rgc mvt feh rc
## 1: NGC 104 00 24 05.2 -72 04 51 305.90 -44.89 14.02488 23.80968 -9.37 -0.76 0.37
## 2: NGC 288 00 52 47.5 -26 35 24 152.28 -89.38 26.41896 37.18224 -6.55 -1.24 1.42
## 3: NGC 362 01 03 14.3 -70 50 54 301.53 -46.25 27.07128 30.00672 -8.35 -1.16 0.17
## 4: NGC 1261 03 12 15.3 -55 13 01 270.54 -52.13 52.18560 58.38264 -7.76 -1.35 0.39
## 5: Pal 1 03 33 23.0 +79 34 50 130.07 19.03 31.63752 51.85944 -1.88 -0.80 0.32
## ---
## 143: NGC 7089 21 33 29.3 -00 49 23 53.38 -35.78 37.18224 33.59448 -8.97 -1.62 0.34
## 144: NGC 7099 21 40 22.0 -23 10 45 27.18 -46.83 25.76664 22.83120 -7.38 -2.12 0.06
## 145: Pal 12 21 46 38.8 -21 15 03 30.51 -47.68 60.99192 50.55480 -4.43 -0.93 0.20
## 146: Pal 13 23 06 44.4 +12 46 19 87.10 -42.70 85.78008 88.71552 -3.46 -1.65 0.48
## 147: NGC 7492 23 08 26.7 -15 36 41 53.39 -63.48 82.19232 79.25688 -5.72 -1.51 0.83rm(raw)Calculate x, y, z coordinates
gc <- cbind(gc, gc2xyz(gc$rsun, gc$l, gc$b))
print(data.table(gc %>% dplyr::select(name, x, y, z)))## name x y z
## 1: NGC 104 5.8262596 -8.0486620 -9.89803
## 2: NGC 288 -0.2530656 0.1329751 -26.41741
## 3: NGC 362 9.7896042 -15.9564226 -19.55532
## 4: NGC 1261 0.3019209 -32.0338534 -41.19561
## 5: Pal 1 -19.2527681 22.8877093 10.31583
## ---
## 143: NGC 7089 17.9934324 24.2105170 -21.73951
## 144: NGC 7099 15.6820143 8.0525420 -18.79231
## 145: Pal 12 35.3783619 20.8477658 -45.09719
## 146: Pal 13 3.1894312 62.9602997 -58.17259
## 147: NGC 7492 21.8864204 29.4593469 -73.54393Filter any globular cluster with distance <80 kly from the center of the Milky Way
gc <- gc %>% filter(rgc <= 80)
nrow(gc)## [1] 130Assume Milky May diameter 110 kly, distance from Sun to center 27 kly
milkyway <- circle(center=c(27, 0), diameter=110)
milkyway$z <- c(0)Galactical coordinates of Sgr A* (center of the Milky Way), and the Sun
sgra <- gc2xyz(r=27, l=0, b=0)
sun <- gc2xyz(r=0, l=0, b=0)Extract M13 from galactic clusters
m13.name <- "NGC 6205"
m13 <- gc %>% filter(name == m13.name)
gc <- gc %>% filter(name != m13.name)Plot clusters position in the Milky Way
ggplot() +
geom_point(data=gc, aes(x, y), color="blue", alpha=0.5, shape=20, size=2) +
geom_path(data=milkyway, aes(x, y), color="black") +
geom_point(data=sgra, aes(x, y), color="black", alpha=1.0, shape=20, size=2) +
annotate("text", x=sgra$x+1, y=sgra$y+1, hjust=0, vjust=0, label="Sgr A*") +
geom_point(data=sun, aes(x, y), color="yellow", alpha=1.0, shape=20, size=2) +
annotate("text", x=sun$x+1, y=sun$y+1, hjust=0, vjust=0, label="Sun") +
geom_point(data=m13, aes(x, y), color="red", alpha=1.0, shape=20, size=2) +
annotate("text", x=m13$x+1, y=m13$y+1, hjust=0, vjust=0, label="M 13") +
labs(title="Glalctic Clusters in the Milky Way (x, y)", x="x [kly]", y="y [kly]") +
coord_equal() +
theme_bw()
ggplot() +
geom_point(data=gc, aes(y, z), color="blue", alpha=0.5, shape=20, size=2) +
geom_path(data=milkyway, aes(y, z), color="black") +
geom_point(data=sun, aes(y, z), color="yellow", alpha=1.0, shape=20, size=2) +
annotate("text", x=sun$y+1, y=sun$z+1, hjust=0, vjust=0, label="Sun") +
geom_point(data=m13, aes(y, z), color="red", alpha=1.0, shape=20, size=2) +
annotate("text", x=m13$y+1, y=m13$z+1, hjust=0, vjust=0, label="M 13") +
labs(title="Galactic Clusters in the Milky Way (y, z)", x="y [kly]", y="z [kly]") +
coord_equal() +
theme_bw()
ggplot() +
geom_point(data=gc, aes(x, z), color="blue", alpha=0.5, shape=20, size=2) +
geom_path(data=milkyway, aes(x, z), color="black") +
geom_point(data=sgra, aes(x, z), color="black", alpha=1.0, shape=20, size=2) +
annotate("text", x=sgra$x+1, y=sgra$z+1, hjust=0, vjust=0, label="Sgr A*") +
geom_point(data=sun, aes(x, z), color="yellow", alpha=1.0, shape=20, size=2) +
annotate("text", x=sun$x+1, y=sun$z+1, hjust=0, vjust=0, label="Sun") +
geom_point(data=m13, aes(x, z), color="red", alpha=1.0, shape=20, size=2) +
annotate("text", x=m13$x+1, y=m13$z+1, hjust=0, vjust=0, label="M 13") +
labs(title="Galactic Clusters in the Milky Way (x, z)", x="x [kly]", y="z [kly]") +
coord_equal() +
theme_bw()
Absolute Visual Magnitude
ggplot() +
geom_histogram(data=gc, aes(mvt), binwidth=0.2, origin=-10.2, color="black", fill="white") +
geom_vline(xintercept=m13$mvt, color="darkgreen") +
annotate("text", x=m13$mvt-0.1, y=12, hjust=1, vjust=0, label="M13", color="darkgreen") +
labs(title="Globular Clusters - Distribution of Absolute Visual Magnitude",
x="Absolute Visual Magnitude [mag]",
y="Number of Stars per Bin") +
annotate("text", x=-4, y=10, hjust=0, vjust=0, label=paste0("N=",nrow(gc))) +
coord_cartesian(ylim=c(0,15)) +
theme_bw()
Number of globular clusters (absolute) brighter than M 13
n <- nrow(gc %>% filter(mvt < m13$mvt))
sprintf("%d %.1f%%", n, n / nrow(gc) * 100)## [1] "18 14.0%"rm(n)Metallicity
Calculate parameters for bimodal distribution
feh.data <- gc$feh[!is.na(gc$feh)]
feh.model <- normalmixEM(feh.data, lambda=0.5, mu=c(-1.5,-0.5), sigma=c(0.5,0.5))## number of iterations= 48feh.model[c("lambda", "mu", "sigma")]## $lambda
## [1] 0.6724072 0.3275928
##
## $mu
## [1] -1.5875239 -0.5785543
##
## $sigma
## [1] 0.3245163 0.2919215ggplot(data=gc, aes(x=feh)) +
geom_histogram(binwidth=0.1, origin=-2.35, color="black", fill="white") +
stat_function(fun=snorm, colour="red",
arg=list(scale=10*feh.model$lambda[1], mean=feh.model$mu[1], sd=feh.model$sigma[1])) +
stat_function(fun=snorm, colour="green",
arg=list(scale=10*feh.model$lambda[2], mean=feh.model$mu[2], sd=feh.model$sigma[2])) +
geom_vline(xintercept=m13$feh, color="darkgreen") +
annotate("text", x=m13$feh-0.1, y=12, hjust=1, vjust=0, label="M13", color="darkgreen") +
annotate("text", x=-0.5, y=10, hjust=0, vjust=0, label=paste0("N=",length(feh.data))) +
labs(title="Globular Clusters - Distribution of Metallicity",
x="Metallicity [Sun]",
y="Number of Stars per Bin") +
coord_cartesian(ylim=c(0,15)) +
theme_bw()
Number of globular clusters with an smaller metallicity than M 13
n <- nrow(gc %>% filter(feh < m13$feh))
sprintf("%d %.1f%%", n, n / nrow(gc) * 100)## [1] "43 33.3%"rm(n)Core Radius
ggplot(gc, aes(x=rc)) +
geom_histogram(binwidth=0.1, origin=-0.05, color="black", fill="white") +
geom_vline(xintercept=m13$rc, color="darkgreen") +
labs(title="Globular Clusters - Distribution of Core Radius",
x="Core Radius [arcmin]",
y="Number of Stars per Bin") +
annotate("text", x=m13$rc+0.1, y=10, hjust=0, vjust=0, label="M13", color="darkgreen") +
annotate("text", x=2.5, y=25, hjust=0, vjust=0, label=paste0("N=",nrow(gc))) +
coord_cartesian(ylim=c(0,30)) +
theme_bw()
Number of globular clusters with an smaller core radius than M 13
n <- nrow(gc %>% filter(rc < m13$rc))
sprintf("%d %.1f%%", n, n / nrow(gc) * 100)## [1] "97 75.2%"rm(n)Color Magnitude Diagram
Read color magnitude raw data from input file
raw <- read.csv("J_AJ_122_3219.tsv", sep=";", as.is=T, strip.white=T, na.strings="", comment.char="#",
col.names=c("ryl",
"v", "e_v",
"bv", "e_bv",
"ra", "de"))
raw$ryl <- as.numeric(raw$ryl)
raw$v <- as.numeric(raw$v)
raw$e_v <- as.numeric(raw$e_v)
raw$bv <- as.numeric(raw$bv)
raw$e_bv <- as.numeric(raw$e_bv)
raw$ra <- as.numeric(raw$ra)
raw$de <- as.numeric(raw$de)Use some more or less clever trick to remove textual annotations from raw data ;)
cm <- raw %>%
filter(!is.na(raw$ryl)) %>%
dplyr::select(ryl, v, bv) %>%
mutate(V=v-5*log10(ly2pc(1000*m13$rsun)))
print(data.table(cm))## ryl v bv V
## 1: 1 14.242 0.573 -4.9834902
## 2: 2 14.308 0.889 -4.9174902
## 3: 3 14.503 0.844 -4.7224902
## 4: 4 14.721 0.597 -4.5044902
## 5: 5 14.736 0.824 -4.4894902
## ---
## 2148: 5418 18.767 -0.324 -0.4584902
## 2149: 5419 18.774 -0.267 -0.4514902
## 2150: 5420 19.084 -0.285 -0.1414902
## 2151: 5422 20.255 0.138 1.0295098
## 2152: 5424 20.769 0.240 1.5435098rm(raw)Plot color magnitude diagram
ggplot(cm, aes(x=bv, y=V)) +
geom_point(aes(color=-bv), alpha=0.5) +
labs(title="M13 - Color Magnitude Diagram", x="Color Index B-V [mag]", y="Absolute Magnitude V [mag]") +
scale_color_gradientn(colours=rainbow(5, start=0, end=0.6), guide = F) +
scale_y_reverse() +
theme_bw()
Appendix
Scaled Normalized Gaussian
snorm <- function(x,scale,mean,sd) {
scale * dnorm(x,mean,sd)
}