Solar eclipse
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Introduction
Today there was a solar eclipse that was not visible on my side of the
Atlantic, but was seen on the European side, either as a partial eclipse,
towards the south, or a total one, towards the north [1]. Eclipses being rare
and solar power being a new thing, this event caused unprecedented reduction of
solar power [2].
A good spot for viewing the total eclipse was Longyearbyen, on Svalbard, and
readers seeking good images might want to include that name in web searches. (I
could not find open-source images at the time of writing, but of course that
was not really my goal in this blog posting…)
I thought it might be interesting to see whether the sun and moon functions in
the oce package could reproduce the eclipse timing, so I constructed a function
to measure the mismatch between sun and moon position in the sky, and set up an
optimization problem to find the time of least mismatch.
The oce functions sunAngle()
and moonAngle()
are at the heart of the
work. Each returns a list that contains, among other things, altitude
and
azimuth
. I set up a mismatch function to calculate a combination of these,
with a scale factor (line 11 of the code) to account for a sort of
converging-meridian effect of azimuth. This function is unity at the horizon,
where a degree change in azimuth is the same angular distance as a degree
change of altitude, and it falls to 0 at the zenith.
Below is the code, and the graph it makes. The black line is the estimated time
when sun and moon centres are nearest each other in the sky, and the time is
written in black near the top of the graph. The red line is the estimate of
eclipse maximum from e.g. [3] The red line is the estimate of eclipse maximum
from [3].
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
library(oce) angle <- function(t, lon=15+40/60, lat=78+12/60) { ## Default location is Longyearbyen, Svalbard sa <- sunAngle(t, lon, lat) ma <- moonAngle(t, lon, lat) saz <- sa$azimuth sal <- sa$altitude maz <- ma$azimuth mal <- ma$altitude C <- cos(0.5*(ma$altitude+sa$altitude) * pi / 180) sqrt((C*(saz-maz))^2 + (sal-mal)^2) } # time is from reference 3 nasa <- as.POSIXct("2015-03-20 9:45:39", tz="UTC") times <- nasa + seq(-1800, 1800, 30) misfit <- NULL for (i in seq_along(times)) { misfit <- c(misfit, angle(times[i])) } oce.plot.ts(times, misfit, ylab="Angular misfit [deg]") o <- optimize(function(t) angle(nasa+t), lower=-3600, upper=3600) eclipse <- nasa + o$minimum abline(v=eclipse, col='black') abline(v=nasa, col='red') mtext(sprintf("Here: %s ", format(eclipse)), line=-1, adj=1, col="black") mtext(sprintf("NASA: %s ", format(nasa)), line=-2, adj=1, col="red") |
Results and discussion
The present method suggests the maximum eclipse occured 3 minutes after the NASA estimate (see [3]).
I’m not too sure why this is, so I’ve put some thoughts in the exercises.
Exercises
-
Verify or correct the “C” factor in line 11 of the code.
-
Determine why the times disagree by a few minutes.
-
Determine whether oce code is now out of synch with UTC. I notice on the
NASA site [3] they show TD (terrestrial dynamical time) being about a minute
after UTC [4]. I wonder whether TD may be a replacement for the time I used
in the oce code, which was based on algorithms from the 1970s?
References and resources
-
Overview of eclipse (wikipedia).
-
Effect of eclipse on power grids (reuters).
-
Image from
NASA
showing eclipse details, including timing. -
Time types at NASA
-
R source code used here: 2015-03-20-eclipse.R.
-
Jekyll source code for this blog entry: 2015-03-20-eclipse.Rmd
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.