Open Cluster Analysis
Ralph Rogge (RRO)
22 May 2016
Load R libraries
library(data.table)## Warning: package 'data.table' was built under R version 3.3.2library(dplyr)
library(ggplot2) ## Warning: package 'ggplot2' was built under R version 3.3.2Parameters
Threshold for open cluster membership
oc.membership.threshold <- 68 # %Thickness of rings around center of open cluster for radial star density
oc.dr <- 1 # arcminProject directories
data.dir <- "../Data/"
output.dir <- "../Output Data/"Create directories when needed
for (d in c(output.dir)) {
if (!dir.exists(d)) dir.create(d)
}
rm(d)Set Up
Load open clusters data DAML02 [1]
daml02.raw <- read.csv(paste0(data.dir, "/B_ocl.csv"),
comment.char="#",
sep=";",
as.is=TRUE,
strip.white=TRUE)
daml02 <- daml02.raw %>%
mutate(name = as.character(Cluster)) %>%
mutate(ra2000 = RAJ2000) %>%
mutate(de2000 = DEJ2000) %>%
mutate(l = as.numeric(X_Glon)) %>%
mutate(b = as.numeric(X_Glat)) %>%
mutate(classification = as.factor(Class)) %>%
mutate(nc = as.numeric(Nc)) %>%
mutate(diameter = as.numeric(Diam)) %>%
mutate(distance = as.numeric(Dist)) %>%
mutate(age = as.numeric(Age)) %>%
mutate(type = as.factor(TrType)) %>%
select(name, l, b, classification, distance, age) %>%
filter(!is.na(l))
print(data.table(daml02))## name l b classification distance age
## 1: Berkeley 58 116.753 -1.289 2700 8.470
## 2: NGC 7801 114.729 -11.315 r 1275 9.230
## 3: FSR 0459 116.468 -2.994 IR 3800 8.950
## 4: Stock 18 117.625 2.268 1242 8.100
## 5: Berkeley 59 118.221 4.996 1000 6.100
## ---
## 2163: Frolov 1 116.559 -0.568 nf 2000 7.588
## 2164: NGC 7790 116.589 -1.009 2944 7.749
## 2165: NGC 7795 116.376 -2.163 2105 8.650
## 2166: FSR 0460 116.575 -1.533 IR 5875 8.705
## 2167: NGC 7762 118.127 5.644 780 9.300nrow(daml02)## [1] 2167Load open clusters data DMCL14 [2]
dmcl14.raw <- read.csv(paste0(data.dir, "J_A+A_564_A79_table2.csv"),
comment.char="#",
sep=";",
as.is=TRUE,
strip.white=TRUE)
dmcl14 <- dmcl14.raw %>%
mutate(name = as.character(Cluster)) %>%
mutate(ra2000 = as.numeric(RAJ2000)) %>%
mutate(de2000 = as.numeric(DEJ2000)) %>%
mutate(l = as.numeric(X_Glon)) %>%
mutate(b = as.numeric(X_Glat)) %>%
mutate(diameter = 2*as.numeric(Rad)) %>%
mutate(n = as.numeric(N)) %>%
mutate(nc = as.numeric(Nc)) %>%
mutate(nf = as.numeric(Nf)) %>%
select(name, ra2000, de2000, diameter, nc) %>%
filter(!is.na(diameter))
print(data.table(dmcl14))## name ra2000 de2000 diameter nc
## 1: SAI 113 155.68333 -59.50555 2.0 18
## 2: Dutra Bica 10 267.57501 -28.89500 2.6 29
## 3: FSR 0357 337.72498 53.95639 2.6 23
## 4: FSR 0023 269.39584 -22.87556 2.7 18
## 5: DC 1 100.50417 -12.97417 2.7 11
## ---
## 1801: Alessi 22 357.02914 36.20500 84.0 1819
## 1802: Mamajek 4 276.50000 -51.00000 86.0 8333
## 1803: Blanco 1 1.02917 -29.83333 72.0 527
## 1804: NGC 6475 268.46249 -34.79333 82.0 43550
## 1805: Collinder 173 120.70417 -46.38334 372.0 157743nrow(dmcl14)## [1] 1805Merge open clusters catalog data of DAML02 and DMCL14
clusters <- daml02 %>%
merge(dmcl14, by="name") %>%
select(name, ra2000, de2000, l, b, classification, distance, diameter, nc, age)
print(data.table(clusters))## name ra2000 de2000 l b classification distance diameter nc age
## 1: AH03 J2011+26.7 302.98749 26.73556 65.737 -3.921 2567 6.5 32 7.35
## 2: Alessi 10 301.19168 -10.47833 31.654 -20.989 513 20.0 182 8.35
## 3: Alessi 12 310.95001 23.79167 67.465 -11.508 537 42.0 2313 8.10
## 4: Alessi 13 50.42500 -36.30000 238.549 -56.981 100 202.0 2752 NA
## 5: Alessi 19 274.61667 12.16667 40.216 12.758 550 90.0 9211 8.00
## ---
## 1796: Waterloo 2 82.00417 40.37167 168.410 3.094 550 6.0 33 8.33
## 1797: Waterloo 3 118.77499 -25.36500 242.562 1.443 5200 4.0 31 NA
## 1798: Waterloo 7 111.53333 -15.09333 230.279 0.622 2800 4.6 45 NA
## 1799: Waterloo 8 112.55000 -15.83333 231.397 1.130 NA 14.0 212 NA
## 1800: Westerlund 2 156.00833 -57.76667 284.276 -0.336 2850 4.0 84 6.30nrow(clusters)## [1] 1800Filter open cluster data by classification
clusters <- clusters %>%
filter(classification != "a") %>% # Possible asterism/dust hole/star cloud (no cluster)
filter(classification != "d") %>% # Objects considered doubtful by the DSS images inspection
filter(classification != "m") %>% # Possible moving group
filter(classification != "nf") %>% # Objects not found in the DSS images inspection
filter(nc >= 5) # Filter clusters w/ less than 5 stars
nrow(clusters)## [1] 1666Calculate open clusters galactic cartesian position
clusters$x <- clusters$distance * cos.deg(clusters$l) * cos.deg(clusters$b) / 1000 # kpc
clusters$y <- clusters$distance * sin.deg(clusters$l) / 1000 # kpc
clusters$z <- clusters$distance * cos.deg(clusters$l) * sin.deg(clusters$b) / 1000 # kpcOpen Cluster Analysis
knitr::include_graphics("DSS/M45.png")
oc <- subset(clusters, name == "Melotte 22")
print(oc %>% select(-x,-y,-z))## name ra2000 de2000 l b classification distance diameter nc age
## 908 Melotte 22 56.75 24.11667 166.571 -23.521 133 122 252 8.131Distance
Extract data for open clusters distance
distance <- clusters$distance
distance <- distance[!is.na(distance)]
length(distance)## [1] 1603Summary statistics of open cluster distance
summary(distance)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 45 1100 1760 2212 2700 14870Distance to Melotte 22
printf("%.0f pc / %.0f ly", oc$distance, as.ly(oc$distance))## [1] "133 pc / 434 ly"p <- histogram.plot(distance, c(0,10000), 100, oc$distance) +
labs(title="Distance of Open Clusters", x="Distance [pc]")
print(p)
Calculate the number of open clusters with distance above or below the value of Melotte 22
abv <- ab(distance, oc$distance)Number of open clusters nearer …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "6 (0.4%)"… and farther than Melotte 22
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "1596 (99.6%)"Diameter
Extract data for open clusters diameter
diameter <- oc.diameter(clusters)
diameter <- diameter[!is.na(diameter)]
length(diameter)## [1] 1603Summary statistics for open clusters diameter
summary(diameter)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.448 3.228 5.044 6.787 8.319 76.910Diameter of Melotte 22
printf("%.1f arcmin / %.0f pc / %.0f ly", oc$diameter, oc.diameter(oc), as.ly(oc.diameter(oc)))## [1] "122.0 arcmin / 5 pc / 15 ly"p <- histogram.plot(diameter, c(0,30), 0.25, oc.diameter(oc)) +
labs(title="Diameter of Open Clusters", x="Diameter [pc]")
print(p)
Calculate the number of open clusters with diameter above or below the value of Melotte 22
abv <- ab(diameter, oc.diameter(oc))Number of open clusters smaller …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "739 (46.1%)"… and larger than Melotte 22
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "863 (53.8%)"Number of Stars
Extract data for open clusters number of stars
numberofstars <- clusters$nc
numberofstars <- numberofstars[!is.na(numberofstars)]
length(numberofstars)## [1] 1666Summary statistics for open clusters number of stars
summary(numberofstars)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 9.0 63.0 133.5 1237.0 420.0 246200.0Number of stars of Melotte 22
printf("%.0f", oc$nc)## [1] "252"p <- histogram.plot(numberofstars, c(0,1000), 10, oc$nc) +
labs(title="Number of Stars of Open Clusters", x="Number of Stars")
print(p)
Calculate the number of open clusters with number of stars below or above the value of Melotte 22
abv <- ab(numberofstars, oc$nc)Number of open clusters with less …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "1100 (66.0%)"… and more stars than Melotte 22
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "563 (33.8%)"Age
Extract open clusters age
age <- clusters$age
age <- age[!is.na(age)]
length(age)## [1] 1580Summary statistics for open clusters age
summary(age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.000 7.757 8.475 8.290 8.900 10.100Age of Melotte 22
printf("%.1f Myr", 10^(oc$age-6))## [1] "135.2 Myr"p <- histogram.plot(age, c(5, 10), 0.1, oc$age) +
labs(title="Age of Open Clusters", x="Age [log yr]")
print(p)
Calculate the number of open clusters with age above or below the value of Melotte 22
abv <- ab(age, oc$age)Number of open clusters younger …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "568 (35.9%)"… and older than Melotte 22
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "1011 (64.0%)"Open Clusters Nearby
The 10 nearest open clusters to Melotte 22
nearest <- clusters %>%
filter(name != oc$name) %>%
mutate(r=sqrt((x-oc$x)^2 + (y-oc$y)^2 + (z-oc$z)^2) * 1000) %>%
arrange(r) %>%
head(10) %>%
select(name, ra2000, de2000, nc, diameter, distance, r)
print(data.table(nearest))## name ra2000 de2000 nc diameter distance r
## 1: Mamajek 3 81.79583 6.26667 22970 252 92 73.00665
## 2: Melotte 20 51.07917 49.86166 53143 302 185 86.91747
## 3: Melotte 25 66.72500 15.86667 18603 332 45 89.87130
## 4: Platais 2 18.45833 32.02833 18917 338 190 119.43271
## 5: Platais 3 70.40417 71.23666 20316 218 200 139.36315
## 6: Alessi 13 50.42500 -36.30000 2752 202 100 147.29803
## 7: Platais 4 76.84167 22.27833 21090 206 276 156.31590
## 8: Melotte 111 186.27501 26.10000 827 122 96 190.90327
## 9: Mamajek 1 130.52501 -79.02723 626 42 97 205.49086
## 10: Stock 2 33.67917 59.48500 3280 62 303 214.34022Membership
Compose file name for open cluster member stars
file.name <- paste0("../Data/J_A+A_564_A79_pm_ucac4-",sub(" ","_",oc$name),".csv")
print(file.name)## [1] "../Data/J_A+A_564_A79_pm_ucac4-Melotte_22.csv"stopifnot(file.exists(file.name))Load open cluster member data [2]
stars.raw <- read.csv(file.name, comment.char="#", sep=";", as.is=TRUE, strip.white=TRUE)
stars <- stars.raw %>%
mutate(name = as.character(UCAC4)) %>%
mutate(ra2000 = as.numeric(RAJ2000)) %>%
mutate(de2000 = as.numeric(DEJ2000)) %>%
mutate(p = as.numeric(P)) %>%
mutate(r = as.numeric(R)) %>%
mutate(pmRA = as.numeric(pmRA)) %>%
mutate(pmDE = as.numeric(pmDE)) %>%
mutate(B = as.numeric(Bmag)) %>%
mutate(V = as.numeric(Vmag)) %>%
mutate(J = as.numeric(Jmag)) %>%
select(name, ra2000, de2000, p, r, pmRA, pmDE, B, V, J) %>%
filter(!is.na(ra2000))
print(data.table(stars))## name ra2000 de2000 p r pmRA pmDE B V J
## 1: 571-008671 56.74844 24.10338 0 0.80 15.4 -7.8 16.888 16.119 14.496
## 2: 571-008667 56.72676 24.11186 0 1.31 35.7 -31.1 15.733 14.507 12.054
## 3: 571-008681 56.77490 24.11167 0 1.40 4.6 -5.9 17.149 16.185 14.198
## 4: 571-008682 56.77718 24.12113 0 1.51 2.0 -4.7 17.635 16.909 15.438
## 5: 571-008677 56.76706 24.13745 0 1.56 4.3 -3.9 NA NA 14.785
## ---
## 4200: 576-009117 57.14420 25.06685 0 60.93 0.6 -8.3 11.538 11.170 10.279
## 4201: 570-008730 57.80947 23.80996 0 60.93 9.7 1.4 17.046 16.153 15.109
## 4202: 574-008350 55.80445 24.65544 96 60.95 23.2 -44.7 NA NA 12.178
## 4203: 566-008101 56.92012 23.11284 0 60.95 -5.2 -3.2 13.132 12.661 11.665
## 4204: 575-008692 55.96730 24.84162 94 60.98 22.4 -43.4 15.785 14.420 11.470nrow(stars)## [1] 4204Membership probability of stars in the UCAC4 catalog near the center of Melotte 22
p <- histogram.plot(stars$p, c(0,100), 1, oc.membership.threshold) +
labs(x="Probability [%]", y="N", title=paste0("Membership (",oc$name,")"))
print(p)
Filter stars with a membership probability > 68%
oc.stars <- stars %>%
filter(p > oc.membership.threshold)
nrow(oc.stars)## [1] 258Calculate open cluster member offset from center
oc.stars <- oc.stars %>%
mutate(dRA = (oc$ra2000 - ra2000) * 60) %>% # arcmin
mutate(dDE = (oc$de2000 - de2000) * 60) # arcminColor Magnitude Diagram
Extract open cluster stars with B and V values defined
cmd <- oc.stars %>%
mutate(BV = B - V) %>%
filter(!is.na(BV)) %>%
select(B, V, BV)
nrow(cmd)## [1] 203p <- ggplot(cmd) +
geom_point(aes(BV,V), alpha=0.3) +
labs(x="B - V [mag]", y="V [mag]", title=paste0("Color Magnitude Diagram (",oc$name,")")) +
coord_cartesian(xlim=c(-0.5, 2.5)) +
scale_y_reverse() +
theme_bw()
print(p)
write.csv(cmd, paste0(output.dir,"Color Magnitude Diagram ",oc$name,".csv"), row.names = FALSE)Radial Star Density
dt <- table(cut_width(oc.stars$r, width=oc.dr, center=0))
n <- as.data.frame(dt)[,2]
r <- as.data.frame(dt, stringsAsFactors=FALSE)[,1]
r <- as.numeric(unlist(strsplit(substring(r,2,nchar(r)-1),",")))
index <- seq(length(r))
r1 <- r[index %% 2 == 1]
r2 <- r[index %% 2 == 0]
r <- (r1 + r2) / 2
d <- n / (pi * abs((r+oc.dr)^2 - r^2))p <- ggplot(data.frame(r=(r1+r2)/2, n=n)) +
geom_bar(aes(r, d), stat="identity") +
geom_hline(yintercept = mean(d)) +
labs(x="r [arcmin]", y="N [1/arcmin^2]", title=paste0("Radial Star Density (",oc$name,")")) +
theme_bw()
print(p)
Appendix
Function ab
Calculate the number of elements in a vector above and below a given value.
ab <- function(data, value, na.rm=TRUE) {
if (na.rm) data <- data[!is.na(data)]
n <- length(data)
n.lo <- sum(data < value)
n.hi <- sum(data > value)
p.lo <- 100 * n.lo / n
p.hi <- 100 * n.hi / n
list(n=n, n.lo=n.lo, n.hi=n.hi, p.lo=p.lo, p.hi=p.hi)
}Function as.ly
Convert parsec into lightyears.
as.ly <- function(x) 3.26*xFunction as.rad
Convert an angle in degrees into rad.
as.rad <- function(alpha) alpha * pi / 180.0Function cos.deg
Calculate the cosine of an angle given in degrees.
cos.deg <- function(alpha) cos(as.rad(alpha))Function histogram.plot
Create a histogram plot (ggplot2).
histogram.plot <- function(data, xlim=NULL, binwidth=1, xintercept=NULL) {
xlim.min <- ifelse(is.null(xlim), min(data, na.rm=TRUE), xlim[1])
xlim.max <- ifelse(is.null(xlim), max(data, na.rm=TRUE), xlim[2])
df <- data.frame(x=data[!is.na(data) & xlim.min <= data & data <= xlim.max])
p <- ggplot(df, aes(x)) +
geom_histogram(binwidth = binwidth) +
ylab("Count") +
theme_bw() +
coord_cartesian(xlim=c(xlim.min,xlim.max))
if (!is.null(xintercept)) p <- p + geom_vline(xintercept = xintercept)
p
}Function oc.diameter
Calculate the spatial diameter of an open cluster.
oc.diameter <- function(oc) oc$diameter * oc$distance * 60 / 206265Function printf
printf <- function(...) {
sprintf(...)
}Function sin.deg
Calculate the sine of an angle given in degrees.
sin.deg <- function(alpha) sin(as.rad(alpha))References
[1] Dias, W S; Alessi, B S; Moitinho, A; et al. (2002): „New catalogue of optically visible open clusters and candidates“. Astronomy and Astrophysics. EDP Sciences 389 (3), pp 871–873, DOI: 10.1051/0004-6361:20020668.
[2] Dias, W S; Monteiro, H; Caetano, T C; et al. (2014): „Proper motions of the optically visible open clusters based on the UCAC4 catalog“. Astronomy and Astrophysics. EDP Sciences 564, p A79, DOI: 10.1051/0004-6361/201323226.