Meade 4000 Filter

Set this flag to TRUE when you want to write out plots.

write.plot.enabled = F

library(data.table)
library(dplyr)
library(ggplot2)
library(jsonlite)

Filter Transmission

Read Meade 4000 filter set data from files

filter <- data.frame(wavelength=c(), transmission=c(), name=c())
for (name in filter.names) {
  filename <- paste0("Meade 4000 Filter ",name,".csv")
  filter <- rbind(filter, preprocess.filter(read.filter(filename), seq(200,900), name))
}

Plot transmission of Meade 4000 filter set

p <- ggplot(data=filter) +
  geom_line(aes(wavelength,transmission,color=factor(name, levels=c(filter.names, "None")))) +
  scale_colour_manual(name="Filter", values=filter.colors, labels=filter.labels) +
  labs(title="Meade 4000 Filter", x="Wavelength [nm]", y="Transmission [%]") +
  theme_bw() +
  coord_cartesian(xlim=c(250,800))
p

if (write.plot.enabled) write.plot(p, "Transmission")

Scotopic/Photopic Vision

Read data for scotopic and photopic vision from file

vision <- rbind(read.vision("scvle.csv", "Scotopic"), read.vision("vljve.csv", "Photopic"))
vision$type = factor(vision$type, levels=c("Scotopic", "Photopic")) # Reorder vision type

Wavelength limits for plots

w.min <- max(min(vision$wavelength), 380)
w.max <- min(max(vision$wavelength), 720)

Plot relative luminous efficency of photopic and scotopic vision

p <- ggplot(data=vision) +
  geom_line(aes(wavelength,rle,color=type)) +
  scale_color_manual(name="Vision", values=c("black","grey70")) +
  labs(title="Vision Efficiency", x="Wavelength [nm]", y="Relative Luminous Efficiency") +
  coord_cartesian(xlim=c(w.min,w.max)) +
  theme_bw()
p

if (write.plot.enabled) write.plot(p)

Visual Observing of Sunlit Object

Read solar spectral distribution from file

sun <- read.csv("ASTMG173.csv", header=F, skip=2)[,c(1,3)]
colnames(sun) <- c("wavelength","irradiance")

Set the average solar irradiance between 630 and 690 nm to one

N <- mean(sun$irradiance[which(sun$wavelength == 630):which(sun$wavelength == 690)])
sun$irradiance <- sun$irradiance / N

Plot solar spectral distribution

p <- ggplot(data=sun) +
  geom_line(aes(wavelength,irradiance)) +
  labs(title="Solar Spectral Irradiance / Air Mass 1.5", x="Wavelength [nm]", y="Irradiance") +
  theme_bw() +
  coord_cartesian(xlim=c(w.min,w.max))
p

Calculate flux for visual observing of sunlit object under scotopic and photopic conditions

sunlit.data <- rbind(
  make.sunlit.data(filter, vision, "Scotopic", sun), 
  make.sunlit.data(filter, vision, "Photopic", sun)
)

sunlit.flux <- sunlit.data %>%

  # Reorder vision.
  mutate(vision = factor(vision, levels=c("Scotopic","Photopic"))) %>%
    
  # Group by vision and filter.
  group_by(vision, filter) %>%
  
  # Add vision specific constant.
  mutate(k = ifelse(vision == "Scotopic", 1699, 683)) %>%
  
  # Calculate flux.
  summarize(flux = sum(fle * irradiance * reflectance * k)) %>%
  
  # Sort by flux.
  arrange(-flux)

print(data.table(sunlit.flux))

##       vision filter        flux
##  1: Scotopic   None 180828.5546
##  2: Scotopic    82A 163391.6940
##  3: Scotopic      8 152820.0533
##  4: Scotopic    80A 104859.6375
##  5: Scotopic     11 104451.7505
##  6: Scotopic     56  62305.8152
##  7: Scotopic     58  45900.9903
##  8: Scotopic    38A  23159.3168
##  9: Scotopic     21   8272.4053
## 10: Scotopic     47   5343.3646
## 11: Scotopic    23A   2684.8845
## 12: Scotopic    25A    795.5685
## 13: Photopic   None  79609.2418
## 14: Photopic      8  73531.9428
## 15: Photopic    82A  68834.3033
## 16: Photopic     11  58546.4630
## 17: Photopic     21  31978.4646
## 18: Photopic    80A  30406.1984
## 19: Photopic     56  28906.3474
## 20: Photopic    23A  20450.6883
## 21: Photopic     58  17103.7880
## 22: Photopic    25A  12034.1782
## 23: Photopic    38A   1524.4048
## 24: Photopic     47    301.1459
##       vision filter        flux

Plot flux for visual observing of a sunlit object

p <- ggplot(sunlit.flux) +
  geom_bar(aes(x=filter,y=flux,fill=factor(filter, levels=c(filter.names, "None"))),stat="identity") +
  scale_fill_manual(name="Filter", values=filter.colors, labels=filter.labels) +
  facet_grid(. ~ vision) +
  labs(title="Sunlit Object", x="Wavelength [nm]", y="Flux [p.d.u.]") +
  theme_bw()
p

if (write.plot.enabled) write.plot(p)

Visual Observing of Jupiter

Read Jupiter NEB/NTrZ reflectance data from file

neb <- read.csv("NEB.csv")
ntrz <- read.csv("NTrZ.csv")

Plot of Jupiter NEB/NTrZ reflectance

p <- ggplot(data=rbind(neb,ntrz)) +
  geom_line(aes(wavelength, reflectance, color=feature)) +
  scale_color_discrete(name="Feature") +
  labs(title="Jupiter Reflectance", x="Wavelength [nm]", y="Reflectance [p.d.u.]") +
  theme_bw() +
  coord_cartesian(xlim=c(w.min,w.max))
p

if (write.plot.enabled) write.plot(p)

Calculate flux and contrast for visual observing Jupiter under scotopic and photopic conditions

jupiter.data <- rbind(
  make.jupiter.data(filter, vision, "Scotopic", neb, ntrz), 
  make.jupiter.data(filter, vision, "Photopic", neb, ntrz)
)

jupiter.flux <- jupiter.data %>%
  
  # Adjust vision order.
  mutate(vision=factor(vision, levels=c("Scotopic","Photopic"))) %>%
  
    # Add vision specific constant.
  mutate(k = ifelse(vision=="Scotopic",1699,683)) %>%

  # Group by vision and filter.
  group_by(vision, filter) %>%
  
  # Calculate flux.
  summarize(flux.neb=sum(k*fle*neb), flux.ntrz=sum(k*fle*ntrz)) %>%
  
  # Calculate Michelson contrast.
  mutate(michelson=(flux.ntrz-flux.neb)/(flux.ntrz+flux.neb)) %>%
  
  # Calculate Weber contrast.
  mutate(weber=(flux.ntrz-flux.neb)/(flux.ntrz)) %>%
  
  # Sort by contrast.
  arrange(-michelson)

print(data.table(jupiter.flux))

##       vision filter    flux.neb   flux.ntrz  michelson      weber
##  1: Scotopic     47   4814.8283   6332.0067 0.13610845 0.23960467
##  2: Scotopic    38A  17855.1740  23450.3667 0.13545865 0.23859724
##  3: Scotopic    80A  80226.5267 101919.6436 0.11909730 0.21284530
##  4: Scotopic    82A 125118.8880 156832.6806 0.11247957 0.20221418
##  5: Scotopic   None 138643.9680 173499.9945 0.11166651 0.20089929
##  6: Scotopic      8 115699.2222 143990.6788 0.10894323 0.19648117
##  7: Scotopic     58  34817.7088  42982.7677 0.10494870 0.18996122
##  8: Scotopic     11  79757.4754  97906.3743 0.10215302 0.18536994
##  9: Scotopic     56  47548.7074  58356.4484 0.10205113 0.18520217
## 10: Scotopic     21   7000.1660   7713.6329 0.04848964 0.09249427
## 11: Scotopic    23A   2356.7722   2534.4421 0.03632430 0.07010218
## 12: Scotopic    25A    722.4019    763.5347 0.02768135 0.05387146
## 13: Photopic     47    286.4265    371.8131 0.12971971 0.22964938
## 14: Photopic    38A   1175.9206   1523.1293 0.12864106 0.22795744
## 15: Photopic     58  13424.7239  15935.4852 0.08551579 0.15755789
## 16: Photopic     56  23106.5216  26939.2607 0.07658466 0.14227336
## 17: Photopic    80A  24616.6349  28679.8316 0.07623764 0.14167436
## 18: Photopic    82A  56606.6272  64794.3680 0.06744377 0.12636501
## 19: Photopic   None  65622.7091  74958.6850 0.06640976 0.12454829
## 20: Photopic     11  48036.3705  54831.0920 0.06605316 0.12392096
## 21: Photopic      8  60616.4020  69057.0836 0.06509181 0.12222760
## 22: Photopic     21  28280.5154  30351.8320 0.03532720 0.06824355
## 23: Photopic    23A  18583.9462  19636.1271 0.02752954 0.05358393
## 24: Photopic    25A  11211.1833  11699.2875 0.02130485 0.04172085
##       vision filter    flux.neb   flux.ntrz  michelson      weber

Plot contrast for visual observing Jupiter

p <- ggplot(jupiter.flux) +
  geom_bar(aes(x=filter,y=michelson,fill=factor(filter, levels=c(filter.names, "None"))),stat="identity") +
  scale_fill_manual(name="Filter", values=filter.colors, labels=filter.labels) +
  facet_grid(. ~ vision) +
  labs(title="Jupiter", x="Wavelength [nm]", y="Contrast") +
  theme_bw()
p

if (write.plot.enabled) write.plot(p)

Appendix

Smart approx function

approx.y <- function(x, y = NULL, xout, rule=2, ...) {
  approx(x,y=y,xout=xout,rule=rule,...)$y
}

Filter names

filter.names <- c(
  "8"  , "11" , "21", "23A", "25A", "38A", "47" , "56" , "58", "80A", "82A"
)

Filter colors

filter.colors <- c(
  "yellow", "yellowgreen", "orange", "red", "red3", "blue", "purple", "green",  "green3", "deepskyblue3",
  "skyblue",  "black"
)

Filter labels

filter.labels <- c(
  "#8", "#11", "#21", "#23A", "#25A", "#38A", "#47", "#56", "#58", "#80A", "#82A", "None"
)

Function to set up data frame for visual observing Jupiter

make.jupiter.data <- function(filter, vision, vision.type, neb, ntrz) {
  
  subvision <- subset(vision, type==vision.type)
  wavelength <- seq(min(filter$wavelength), max(filter$wavelength))
  
  # Create data frame with wavelength, transmission, filter.
  df <- data.frame(wavelength = filter$wavelength, transmission = filter$transmission, filter = filter$name)
  
  # Add plain vision data.
  df <- rbind(df, data.frame(wavelength = wavelength, transmission = 100, filter = "None"))
  
  # Add relative luminous efficiency.
  df$rle <- approx.y(subvision$wavelength, subvision$rle, df$wavelength)

  # Add filtered luminous efficiency.
  df$fle <- df$rle * df$transmission / 100

  # Add Jupiter NEB reflectance.
  df$neb <- approx.y(neb$wavelength, neb$reflectance, df$wavelength)

  # Add Jupiter NTrZ reflectance.
  df$ntrz <- approx.y(ntrz$wavelength, ntrz$reflectance, df$wavelength)

  # Add vision.
  df$vision <- vision.type

  df
}

Function to set up data frame for visual observing of a sunlight object

make.sunlit.data <- function(filter, vision, vision.type, sun, reflectance=1) {
  
  subvision <- subset(vision, type==vision.type)
  wavelength <- seq(min(filter$wavelength), max(filter$wavelength))
  
  # Create data frame with wavelength, transmission, filter.
  df <- data.frame(wavelength = filter$wavelength, transmission = filter$transmission, filter = filter$name)

  # Add plain vision data.
  df <- rbind(df, data.frame(wavelength = wavelength, transmission = 100, filter = "None"))

  # Add relative luminous efficiency.
  df$rle <- approx.y(subvision$wavelength, subvision$rle, df$wavelength)

  # Add filtered luminous efficiency
  df$fle <- df$rle * df$transmission/100 

  # Add solar spectral irrandiance
  df$irradiance <- approx.y(sun$wavelength, sun$irradiance, df$wavelength)

  # Add reflectance for a spectral uniform reflector.
  df$reflectance <- reflectance

  # Add vision.
  df$vision <- vision.type

  # Return.
  df
}

Function to preprocess filter data

preprocess.filter <- function(df, wavelength, name) {
  data.frame(
    wavelength = wavelength, 
    transmission = approx.y(df$wavelength, df$transmission, wavelength), 
    name=name
  )
}

Function to read in filter data

read.filter <- function(name) {
  read.csv(name, col.names=c("wavelength","transmission","X"))[,1:2]
}

Function to read vision data from file

read.vision <- function(filename, type) {
  df <- read.csv(filename)
  df$type <- type
  df
}

Function used to write plot to PNG file

write.plot <- function(p, name=NULL, width=1024, height=576, font.size=16) {
  title <- ifelse(!is.null(name), name, p$labels$title)
  if (is.null(title)) stop("Neither name nor title given")
  filename <- paste0("Meade 4000 - ", title, ".png")
  png(filename,width, height)
  print(p + theme_bw(base_size=font.size))
  dev.off()
}