Load external libraries

library(data.table)
library(dplyr)
library(ggplot2)
library(knitr)
library(png)
source("R/vizor.R")

Telescopes

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)

Telescope Properties

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)))

Double Star Observing

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.

Variable Star Observing

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.

Summary

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.

Appendix

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))