Load external libraries
library(data.table)
library(dplyr)
library(ggplot2)
library(knitr)
library(png)
source("R/vizor.R")
Read telescope data from file and preprocess
raw <- read.csv("telescopes.csv", na.strings="NA", as.is=TRUE, strip.white=TRUE)
T <- raw %>%
mutate(name = F(name)) %>%
mutate(N = round(focal.length / aperture, 1)) %>%
mutate(reflectivity = N(reflectivity)) %>%
mutate(central.obstruction = N(central.obstruction)) %>%
mutate(transmission = N(transmission))
Print telescope data
## Telescope Aperture [mm] Focal Length [mm] f/ Price [€]
## N 254 254 1250 4.9 3750
## N 200 200 1000 5.0 2000
## MC 809 200 2000 10.0 4150
## MN 152 152 1200 7.9 2400
## R 152 152 1200 7.9 3500
## R 140 140 890 6.4 3100
## R 130 130 910 7.0 1900
Save telescope data to file
write.csv(DF, "output/telescopes.csv", row.names=FALSE)
Wavelength of maximum sensitivity of photopic vision
lambda <- photopic.vision$wavelength[which.max(photopic.vision$v)]
Telescope Dawes limit for wavelength at maximum sensitivity of photopic vision at 555 nm
T$limiting.r <- dawes.limit(T$aperture, lambda)
Telescope limiting magnitude with effective apperture
obstruction <- ifelse (T$nmirrors>0, T$aperture*T$central.obstruction/100, 0)
aperture.eff <- sqrt(T$aperture^2 - obstruction^2)
T$limiting.V <- telescope.limiting.magnitude(aperture.eff, "north")
Print telescope properties
## Telescope Dawes Limit [arcsec] Limiting Magnitude [mag]
## N 254 0.46 14.99
## N 200 0.58 14.53
## MC 809 0.58 14.52
## MN 152 0.77 14.08
## R 152 0.77 14.10
## R 140 0.83 13.94
## R 130 0.90 13.80
Save telescope properties to file
write.csv(DF, "output/properties.csv", row.names=FALSE)
Points to assign in telescope ranking
POINTS <- c(10, 7, 5, 4, 3, 2, rep(1, ifelse(nrow(T)>6, nrow(T)-6, 0)))
Read WDS catalog data from file and preprocess
raw <- read.csv("data/WDS.csv", sep=";", comment="#", as.is=TRUE, strip.white=TRUE)
WDS <- raw %>%
mutate(observed.at = N(Obs2)) %>%
mutate(PA = N(pa2)) %>%
mutate(r = N(sep2)) %>%
mutate(V.A = N(mag1)) %>%
mutate(V.B = N(mag2)) %>%
mutate(dV = V.B-V.A) %>%
select(designation=WDS, observed.at, PA, r, V.A, V.B, dV) %>%
filter(complete.cases(.))
nrow(WDS)
## [1] 141006
Print excerpt from WDS catalog
## WDS Date PA [°] r [arcsec] V A [mag] V B [mag] ΔV [mag]
## 00000+7530 1982 235 0.6 10.27 11.5 1.23
## 00000+4004 2014 253 4.2 12.10 13.1 1.00
## 00000+4004 2014 66 14.5 12.10 14.1 2.00
## 00000+3852 2015 108 6.6 6.62 11.4 4.78
## 00000+0044 2007 30 1.2 17.80 17.8 0.00
## ---
## 23599+1413 2007 173 1.6 16.30 16.6 0.30
## 23599+0048 2013 268 20.0 11.93 13.9 1.97
## 23599-0412 2010 151 77.0 16.50 17.1 0.60
## 23599-1818 1958 272 2.0 18.60 18.7 0.10
## 23599-3112 1991 248 0.7 11.54 11.9 0.36
Where WDS is the WDS designation, Date is the year of observation, PA is the position angle, r is the seperation, V A magnitude of component A, V B magnitude of component B.
Number of double stars visible
filter.ds <- function(r, V, dV) { sum(WDS$r > r & WDS$V.A < V & WDS$V.B < V & WDS$dV < dV)}
n.ds <- apply(T, 1, function(t) filter.ds(N(t['limiting.r']), N(t['limiting.V']), dV.ds))
Ranking of telescopes with respect to double star observing
ranking.ds <- data.frame(telescope = T$name,
n = n.ds,
npp = round(n.ds/T$price,2),
p = round(100*n.ds/nrow(WDS),2)) %>%
mutate(rank = min_rank(desc(npp))) %>%
arrange(rank) %>%
mutate(points = POINTS[1:nrow(T)]) %>%
select(rank, points, telescope, n, npp, p)
Plot of telescope ranking
p <- ggplot(data=ranking.ds) +
geom_bar(aes(x=telescope, y=npp), stat="identity") +
labs(title="Double Stars Criterion", x="Telescope", y=paste0("N/Price [1/",currency,"]")) +
theme_bw()
print(p)
Save plot to file
png(filename="output/telescopes+wds.png", width=1024, height=800)
print(p + theme(text = element_text(size=20)))
invisible(dev.off())
Print telescope ranking
DF <- ranking.ds
names(DF) <- c("Rank", "Points", "Telescope", "N", paste0("N/Price [1/",currency,"]"), "P [%]")
print(data.table(DF), row.names=FALSE)
## Rank Points Telescope N N/Price [1/€] P [%]
## 1 10 N 200 27813 13.91 19.72
## 2 7 R 130 21351 11.24 15.14
## 3 5 MN 152 22909 9.55 16.25
## 4 4 N 254 31675 8.45 22.46
## 5 3 R 140 21560 6.95 15.29
## 6 2 MC 809 27807 6.70 19.72
## 7 1 R 152 23081 6.59 16.37
write.csv(DF, "output/telescopes+wds.csv", row.names=FALSE)
Where N is the number of splittable double stars, N/Price is the number of splittable double stars per €, and P is the fraction of splittable double stars with respect to the WDS catalog.
Read raw GCVS catalog data from file and preprocess
raw <- read.csv("data/GCVS.csv", sep=";", comment="#", as.is=TRUE, strip.white=TRUE)
GCVS <- raw %>%
mutate(V.max = N(magMax)) %>%
mutate(V.min = N(Min1)) %>%
mutate(type = F(VarTypeII)) %>%
select(designation=GCVS, V.max, V.min, type) %>%
filter(complete.cases(.))
nrow(GCVS)
## [1] 51597
Print excerpt from GCVS variable star catalog
## GCVS V max [mag] V min [mag] Type
## R And 5.80 15.20 M
## S And 5.80 16.00 SN
## T And 7.70 14.50 M
## U And 9.00 15.00 M
## V And 9.00 15.20 M
## ---
## V0548 Vul 14.90 0.50
## V0549 Vul 15.98 0.22
## V0550 Vul 13.00 14.10
## V0551 Vul 11.95 12.40
## V0552 Vul 14.68 15.28
Where GCVS is the GCVS designation, V max is the star’s magnitude at maximum, and V min is the star’s magnitude at minimum.
Select Mira type variables only
GCVS.M <- GCVS %>% filter(type %in% "M")
nrow(GCVS.M)
## [1] 1644
Number of variable stars of Mira type visible during their complete light curve cycle allowing a magnitude safty of 0.5 mag for clean brightness estimates.
filter.vs <- function(V, dV) { sum(GCVS.M$V.min < V-dV) }
n.vs <- apply(T, 1, function(t) filter.vs(N(t['limiting.V']), dV.vs))
Telescope ranking with respect to variable star observing
ranking.vs <- data.frame(telescope = T$name,
n = n.vs,
npp = round(n.vs/T$price,2),
p = round(n.vs/nrow(GCVS.M)*100,2)) %>%
mutate(rank = min_rank(desc(npp))) %>%
arrange(rank) %>%
mutate(points = POINTS[1:nrow(T)]) %>%
select(rank, points, telescope, n, npp, p)
Plot of telescope ranking
p <- ggplot(data=ranking.vs) +
geom_bar(aes(x=telescope, y=npp), stat="identity") +
labs(title="Variable Stars Criterion", x="Telescope", y=paste0("N/Price [1/",currency,"]")) +
theme_bw()
print(p)
Save plot to file
png(filename="output/telescopes+gcvs.png", width=1024, height=800)
print(p + theme(text = element_text(size=20)))
invisible(dev.off())
Print of telescope ranking
DF <- ranking.vs
names(DF) <- c("Rank", "Points", "Telescope", "N", paste0("N/Price [1/",currency,"]"), "P [%]")
print(data.table(DF), row.names=FALSE)
## Rank Points Telescope N N/Price [1/€] P [%]
## 1 10 N 200 327 0.16 19.89
## 2 7 N 254 474 0.13 28.83
## 3 5 MC 809 327 0.08 19.89
## 4 4 MN 152 175 0.07 10.64
## 4 3 R 130 138 0.07 8.39
## 6 2 R 152 193 0.06 11.74
## 7 1 R 140 146 0.05 8.88
write.csv(DF, "output/telescopes+gcvs.csv", row.names=FALSE)
Where N is the number of Mira type variable stars observable during their complete light curve cycle (this is from minimum to maximum), N/Price is the number of observable variable stars per €, and P is the fraction of observable variable stars with respect to the number of Mira type stars found in GCVS catalog.
weight.ds <- 1
weight.vs <- 3
For total ranking double star observing to variable star observing is weighted 1: 3.
Calculate overall ranking
ranking.overall <- merge(ranking.ds, ranking.vs, by="telescope", suffixes=c(".ds", ".vs")) %>%
mutate(points = (points.ds * weight.ds + points.vs * weight.vs) / (weight.ds + weight.vs)) %>%
mutate(rank = row_number(desc(points))) %>%
arrange(rank) %>%
select(rank, points, telescope)
Print overall ranking
DF <- ranking.overall
names(DF) <- c("Rank", "Points", "Telescope")
print(data.table(DF), row.names=FALSE)
## Rank Points Telescope
## 1 10.00 N 200
## 2 6.25 N 254
## 3 4.25 MC 809
## 4 4.25 MN 152
## 5 4.00 R 130
## 6 1.75 R 152
## 7 1.50 R 140
Where Rank is total ranking, Points weighted average points from double star observing and variable star observing ranking, and Telescope is telescope’s name.
Simple helper functions to convert a value to a number and factor
N <- function(x) suppressWarnings(as.numeric(x))
F <- function(x) suppressWarnings(as.factor(x))