Calculate sunrise on Dec 12, 2014, at Halifax.

Abstract. Demo that Oce gives “correct” sunrise times, if the useRefraction argument is used.

Goal. Check to see if we can recover the stated sunrise at 7:34 local time.

Method. Use uniroot. There are other ways, if we want to search through a longer time interval, but we know the longitude so we know roughly what time to examine. THe problem is of course that uniroot will pick a sunset as happily as it picks a sunrise. A simple fix is to add d(altitude)/dt to the claculation.

1. Set some things up.

library(oce)

## Loading required package: mapproj
## Loading required package: maps

t1 <- as.numeric(ISOdatetime(2014, 12, 2,  8, 0, 0, tz="UTC"))
t2 <- as.numeric(ISOdatetime(2014, 12, 2, 12, 0, 0, tz="UTC"))
longitude <- -(63 + 43.2/60)
latitude <-44 + 55.2/60

2. Calculate without accounting for refraction. This is not the right way to do it.

root <- uniroot(function(t) {
                sa <- sunAngle(t, longitude, latitude)
                sa$altitude
                },
                lower=t1, upper=t2)$root
root

## [1] 1417520351

numberAsPOSIXct(root)

## [1] "2014-12-02 11:39:10 UTC"

3. Calculate with refraction of the sun’s rays in the lower sky. This is still not quite the way to do it.

root <- uniroot(function(t) {
                sa <- sunAngle(t, longitude, latitude, useRefraction=TRUE)
                sa$altitude
                },
                tol=1e-9,
                lower=t1, upper=t2)$root
sprintf("%.2f", root)

## [1] "1417520046.15"

numberAsPOSIXct(root)

## [1] "2014-12-02 11:34:06 UTC"

4. Calculate with refraction, and also add half the sun diameter to the altitude, so that we find when the sun’s edge just breaks the horizon.

root <- uniroot(function(t) {
                sa <- sunAngle(t, longitude, latitude, useRefraction=TRUE)
                sa$altitude+sa$diameter/2
                },
                tol=1e-9,
                lower=t1, upper=t2)$root
sprintf("%.2f", root)

## [1] "1417520046.15"

numberAsPOSIXct(root)

## [1] "2014-12-02 11:34:06 UTC"

NOTES

I had expected calc 4 should be a few seconds later than calc 3, because it takes the sun a few seconds to rise over the half-disk height, but I suppose what is happening is that the useRefraction calculation may be handling that. If you put 3*sa$diameter/2 you’ll see a time shift, though. What happens with the refraction calculation is that the sun jumps position – hard to explain without graphing it, but in any case the formulae are from standard sources, listed in the Oce docs.