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/M36.png")
oc <- subset(clusters, name == "NGC 1960")
print(oc %>% select(-x,-y,-z))## name ra2000 de2000 l b classification distance diameter nc age
## 965 NGC 1960 84.075 34.14 174.534 1.072 1330 12 199 7.4Distance
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 NGC 1960
printf("%.0f pc / %.0f ly", oc$distance, as.ly(oc$distance))## [1] "1330 pc / 4336 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 NGC 1960
abv <- ab(distance, oc$distance)Number of open clusters nearer …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "543 (33.9%)"… and farther than NGC 1960
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "1058 (66.0%)"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 NGC 1960
printf("%.1f arcmin / %.0f pc / %.0f ly", oc$diameter, oc.diameter(oc), as.ly(oc.diameter(oc)))## [1] "12.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 NGC 1960
abv <- ab(diameter, oc.diameter(oc))Number of open clusters smaller …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "718 (44.8%)"… and larger than NGC 1960
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "884 (55.1%)"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 NGC 1960
printf("%.0f", oc$nc)## [1] "199"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 NGC 1960
abv <- ab(numberofstars, oc$nc)Number of open clusters with less …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "1014 (60.9%)"… and more stars than NGC 1960
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "651 (39.1%)"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 NGC 1960
printf("%.1f Myr", 10^(oc$age-6))## [1] "25.1 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 NGC 1960
abv <- ab(age, oc$age)Number of open clusters younger …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "273 (17.3%)"… and older than NGC 1960
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "1298 (82.2%)"Open Clusters Nearby
The 10 nearest open clusters to NGC 1960
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: FSR 0807 84.14166 31.85556 63 10.0 1300 60.38177
## 2: FSR 0775 81.39583 34.95750 63 8.3 1300 61.66832
## 3: FSR 0817 84.86250 30.89333 15 3.2 1359 82.96938
## 4: NGC 1912 82.16667 35.84833 694 22.0 1400 89.09545
## 5: NGC 2099 88.07500 32.55333 588 16.0 1383 102.30938
## 6: ASCC 15 80.65417 36.55000 649 26.0 1400 112.22239
## 7: FSR 0744 74.87500 38.01167 18 3.6 1250 204.17786
## 8: NGC 1778 77.01666 37.02333 116 10.0 1469 209.08338
## 9: NGC 1996 84.54166 25.81667 406 24.0 1400 210.28986
## 10: Dol Dzim 4 83.97500 25.95000 84 14.0 1220 212.60693Membership
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-NGC_1960.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: 621-025835 84.07970 34.13323 10 0.47 -0.9 -14.7 14.570 13.985 12.817
## 2: 621-025818 84.06582 34.14350 98 0.50 -2.1 -4.8 8.823 8.820 8.763
## 3: 621-025815 84.06445 34.13441 0 0.62 NA NA NA NA 14.976
## 4: 621-025848 84.08747 34.13584 37 0.67 4.1 -12.1 NA NA 11.395
## 5: 621-025809 84.06024 34.13997 99 0.73 -1.5 -6.3 NA NA 11.730
## ---
## 272: 621-025713 83.99681 34.06578 93 5.91 1.1 -1.4 14.222 13.044 10.954
## 273: 622-025006 84.06572 34.23833 99 5.92 0.8 -4.5 11.808 11.300 10.264
## 274: 621-025699 83.98109 34.07801 94 5.97 3.9 -5.7 NA NA 13.048
## 275: 622-024955 84.00912 34.22360 71 5.99 1.1 -16.1 NA NA 13.909
## 276: 621-025736 84.01273 34.05443 83 5.99 -4.4 5.4 NA NA 14.623nrow(stars)## [1] 276Membership probability of stars in the UCAC4 catalog near the center of NGC 1960
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] 152Calculate 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] 104p <- 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.