Set Up

Load R libraries

library(astroxo)
library(FITSio)
library(ggplot2)

Create data directories

analytical.data.dir <- "Analytical Data/"
output.data.dir <- "Output Data/"
raw.data.dir <- "Raw Data/"

for (d in c(analytical.data.dir, output.data.dir)) {
  if (!dir.exists(d)) dir.create(d)
}

rm(d)

Read Moon edge detected image and Earth’s shadow edge detected image from file

m.1923 <- read.fits(paste0(raw.data.dir, "1923 Moon Edges.fit"))
s.1923 <- read.fits(paste0(raw.data.dir, "1923 Shadow Edges.fit"))

m.2072 <- read.fits(paste0(raw.data.dir, "2072 Moon Edges.fit"))
s.2072 <- read.fits(paste0(raw.data.dir, "2072 Shadow Edges.fit"))

Extract Moon and Earth’s shadow edge data from images

m.1923$edges <- which(m.1923$image>0, arr.ind=TRUE)
s.1923$edges <- which(s.1923$image>0, arr.ind=TRUE)

m.2072$edges <- which(m.2072$image>0, arr.ind=TRUE)
s.2072$edges <- which(s.2072$image>0, arr.ind=TRUE)

Fit circle to Moon and Earth’s shadow edge data

m.1923$circle <- fit.circle(m.1923$edges)
print(m.1923$circle$params)
## $x
## [1] 474.3077
## 
## $y
## [1] 490.9698
## 
## $r
## [1] 398.3104
s.1923$circle <- fit.circle(s.1923$edges)
print(s.1923$circle$params)
## $x
## [1] -250.7275
## 
## $y
## [1] 21.19255
## 
## $r
## [1] 897.515
m.2072$circle <- fit.circle(m.2072$edges)
print(m.2072$circle$params)
## $x
## [1] 513.3296
## 
## $y
## [1] 497.2333
## 
## $r
## [1] 395.1695
s.2072$circle <- fit.circle(s.2072$edges)
print(s.2072$circle$params)
## $x
## [1] 1444.992
## 
## $y
## [1] 216.4748
## 
## $r
## [1] 870.4163

Normalize circle points and edge coordinates such that shadow cent is at 0,0

m.1923.px <- m.1923$circle$points[,1] - s.1923$circle$params$x
m.1923.py <- m.1923$circle$points[,2] - s.1923$circle$params$y

m.1923.ex <- m.1923$edges[,1] - s.1923$circle$params$x
m.1923.ey <- m.1923$edges[,2] - s.1923$circle$params$y

s.1923.px <- s.1923$circle$points[,1] - s.1923$circle$params$x
s.1923.py <- s.1923$circle$points[,2] - s.1923$circle$params$y

s.1923.ex <- s.1923$edges[,1] - s.1923$circle$params$x
s.1923.ey <- s.1923$edges[,2] - s.1923$circle$params$y

m.2072.px <- m.2072$circle$points[,1] - s.2072$circle$params$x
m.2072.py <- m.2072$circle$points[,2] - s.2072$circle$params$y

m.2072.ex <- m.2072$edges[,1] - s.2072$circle$params$x
m.2072.ey <- m.2072$edges[,2] - s.2072$circle$params$y

s.2072.px <- s.2072$circle$points[,1] - s.2072$circle$params$x
s.2072.py <- s.2072$circle$points[,2] - s.2072$circle$params$y

s.2072.ex <- s.2072$edges[,1] - s.2072$circle$params$x
s.2072.ey <- s.2072$edges[,2] - s.2072$circle$params$y

Estimate error for Earth’s shadow radius

m.1923$circle$params$r.err <- sd(with (m.1923, {
  circle$params$r - sqrt((circle$params$x - edges[,1])^2 + (circle$params$y - edges[,2])^2)
}))

s.1923$circle$params$r.err <- sd(with (s.1923, {
  circle$params$r - sqrt((circle$params$x - edges[,1])^2 + (circle$params$y - edges[,2])^2)
}))
s.1923$circle$params$r.err
## [1] 7.048667
m.2072$circle$params$r.err <- sd(with (m.2072, {
  circle$params$r - sqrt((circle$params$x - edges[,1])^2 + (circle$params$y - edges[,2])^2)
}))

s.2072$circle$params$r.err <- sd(with (s.2072, {
  circle$params$r - sqrt((circle$params$x - edges[,1])^2 + (circle$params$y - edges[,2])^2)
}))
s.2072$circle$params$r.err
## [1] 8.814526

Calculate Earth’s shadow path boundaries

index <- round(314/4)
x1 <- c(m.2072.px[index], m.1923.px[index])
y1 <- c(m.2072.py[index], m.1923.py[index])
lm(y1 ~ x1)
## 
## Call:
## lm(formula = y1 ~ x1)
## 
## Coefficients:
## (Intercept)           x1  
##     782.711        0.116
index <- round(3*314/4)
x2 <- c(m.2072.px[index], m.1923.px[index])
y2 <- c(m.2072.py[index], m.1923.py[index])
lm(y2 ~ x2)
## 
## Call:
## lm(formula = y2 ~ x2)
## 
## Coefficients:
## (Intercept)           x2  
##    -10.0990       0.1122

Plot of Moon and Earth’s shadow

p <- ggplot() +
  geom_smooth(aes(x=x1, y=y1), method="lm", formula=y~x, se=FALSE) +
  geom_smooth(aes(x=x2, y=y2), method="lm", formula=y~x, se=FALSE) +
  geom_path(aes(x=m.1923.px, y=m.1923.py), color="red") +
  geom_point(aes(x=m.1923.ex, y=m.1923.ey), color="red") +
  geom_path(aes(x=s.1923.px, y=s.1923.py)) +
  geom_point(aes(x=s.1923.ex, y=s.1923.ey)) +
  geom_path(aes(x=m.2072.px, y=m.2072.py), color="blue") +
  geom_point(aes(x=m.2072.ex, y=m.2072.ey), color="blue") +
  geom_path(aes(x=s.2072.px, y=s.2072.py)) +
  geom_point(aes(x=s.2072.ex, y=s.2072.ey)) +
  coord_cartesian(xlim=c(-1500,1500), ylim=c(-1500,1500)) +
  xlab("X") +
  ylab("Y") +
  theme_bw()
print(p)

png(paste0(output.data.dir,"Lunar Eclipse 2015-09-28 Fit.png"))
print(p)
dev.off()
## quartz_off_screen 
##                 2
rm(p)

Size of Earth’s Moon Shadow

As Moon’s diameter is

moon.diameter <- 2 * 1737 # km
moon.diameter
## [1] 3474

my imaging system (TMB 80/600 + Canon 10D) has a resolution of

resolution <- moon.diameter / (2 * m.1923$circle$params$r) # km/pixel
signif(resolution,4)
## [1] 4.361

Earth’s shadow diameter

s.1923$diameter = 2 * s.1923$circle$params$r * resolution # km
signif(s.1923$diameter,4)
## [1] 7828
s.1923$diameter.err <- 2 * s.1923$circle$params$r.err * resolution # km
round(s.1923$diameter.err)
## [1] 61
s.2072$diameter = 2 * s.2072$circle$params$r * resolution # km
signif(s.2072$diameter,4)
## [1] 7592
s.2072$diameter.err <- 2 * s.2072$circle$params$r.err * resolution # km
round(s.2072$diameter.err)
## [1] 77
delta <- abs(s.1923$diameter - s.2072$diameter) # km
round(delta)
## [1] 236

Appendix

fit.circle <- function(xy, a0=mean(xy[,1]), b0=mean(xy[,2]), r0=mean(sqrt((xy[,1]-a0)^2+(xy[,2]-b0)^2)), ...) {
  sos <- function(abr){
    sum((abr[3] - sqrt((xy[,1]-abr[1])^2 + (xy[,2]-abr[2])^2))^2)
  }
  fit <- optim(c(a0,b0,r0), sos, ...)
  
  result <- list()
  result$params <- list(x=fit$par[1],y=fit$par[2],r=fit$par[3])
  result$points <- circle(c(fit$par[1], fit$par[2]), 2*fit$par[3])
  result
}
read.fits <- function(filename) {
  data <- list()
  fits <- readFITS(filename)
  data$header <- fits$header 
  data$image <- fits$imDat
  data
}