M7
Ralph Rogge (RRO)
02 June 2016
Load R libraries
library(data.table)
library(dplyr)
library(ggplot2) Parameters
Threshold for open cluster membership
oc.membership.threshold <- 68 # %Thickness of rings around center of open cluster for radial star density
oc.dr <- 1 # arcminSet Up
Load open clusters data DAML02 [1]
daml02.raw <- read.csv("../Data/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("../Data/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/M7.png")
oc <- subset(clusters, name == "NGC 6475")
print(oc %>% select(-x,-y,-z))## name ra2000 de2000 l b classification distance diameter nc age
## 1213 NGC 6475 268.4625 -34.79333 355.861 -4.501 301 82 43550 8.475Distance
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 6475
printf("%.0f pc / %.0f ly", oc$distance, as.ly(oc$distance))## [1] "301 pc / 981 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 6475
abv <- ab(distance, oc$distance)Number of open clusters nearer …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "26 (1.6%)"… and farther than NGC 6475
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "1576 (98.3%)"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 6475
printf("%.1f arcmin / %.0f pc / %.0f ly", oc$diameter, oc.diameter(oc), as.ly(oc.diameter(oc)))## [1] "82.0 arcmin / 7 pc / 23 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 6475
abv <- ab(diameter, oc.diameter(oc))Number of open clusters smaller …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "1108 (69.1%)"… and larger than NGC 6475
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "494 (30.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 NGC 6475
printf("%.0f", oc$nc)## [1] "43550"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 6475
abv <- ab(numberofstars, oc$nc)Number of open clusters with less …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "1658 (99.5%)"… and more stars than NGC 6475
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "7 (0.4%)"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 6475
printf("%.1f Myr", 10^(oc$age-6))## [1] "298.5 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 6475
abv <- ab(age, oc$age)Number of open clusters younger …
printf("%d (%.1f%%)", abv$n.lo, abv$p.lo)## [1] "789 (49.9%)"… and older than NGC 6475
printf("%d (%.1f%%)", abv$n.hi, abv$p.hi)## [1] "789 (49.9%)"Open Clusters Nearby
The 10 nearest open clusters to NGC 6475
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: ASCC 99 282.2708 -18.73000 2948 29.6 280 104.5280
## 2: Alessi 9 265.8750 -46.96667 54391 152.0 211 105.1324
## 3: Ruprecht 145 282.6500 -18.25000 4545 37.0 320 114.2804
## 4: BH 223 260.1708 -35.88334 58 8.0 414 120.5671
## 5: Mamajek 4 276.5000 -51.00000 8333 86.0 385 132.0970
## 6: Ruprecht 147 289.1750 -16.28333 933 27.0 295 136.4233
## 7: Mamajek 2 264.5000 -8.11167 78 12.0 161 173.4562
## 8: Collinder 350 267.0250 1.30000 1282 41.0 280 180.7383
## 9: NGC 6573 273.4208 -22.11833 58 3.8 460 181.2214
## 10: Collinder 359 270.2750 2.90000 70733 242.0 249 183.9903Membership
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_6475.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: 273-130019 268.0848 -35.40185 95 40.94 -4.2 -1.5 NA NA 14.716
## 2: 273-130057 268.0884 -35.40337 0 40.94 36.9 -19.4 NA NA 13.125
## 3: 273-130071 268.0897 -35.40185 0 40.83 44.3 52.3 NA NA 12.539
## 4: 273-130083 268.0907 -35.40508 0 40.98 -68.0 -41.8 NA NA NA
## 5: 273-130164 268.0978 -35.40246 82 40.69 14.1 -15.0 NA NA 13.796
## ---
## 73346: 280-141603 268.8619 -34.19564 0 40.95 -5.6 53.7 NA NA NA
## 73347: 280-141606 268.8623 -34.19686 0 40.90 NA NA NA NA 8.676
## 73348: 280-141611 268.8626 -34.19775 0 40.86 NA NA NA NA NA
## 73349: 280-141650 268.8666 -34.19703 0 40.99 -117.8 -37.9 NA NA 13.723
## 73350: 280-141678 268.8692 -34.19851 0 40.98 46.3 132.7 NA NA 14.154nrow(stars)## [1] 73350Membership probability of stars in the UCAC4 catalog near the center of NGC 6475
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] 39314Calculate 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] 3196p <- 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)
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.