The Matterhorn at Night

While I’m still working on processing multiple deep-sky images as well as the new book, let’s take a brief detour to Switzerland.  My vacation to Zermatt was unfortunately timed with the full moon, so the skies were not as dark as I’d like, but the rising moon did nicely light up the east face.

Samyang 135mm on Canon 6D (2 seconds at f/2.0, ISO 3200)

I had hoped to capture the full moon with the mountain, but to get the right angle, I would have had to have hiked up in the mountains around 3 am. Needless to say, that did not happen.

 

What part of the sky spends the most time in darkness?

SUMMARY:  I set out to answer the question: considering the perspective of everyone on Earth and accounting for the seasonal changes in daylight, what spots in the night sky are least and most visible? The incremental steps to get there illustrate the effect of latitude and the seasons on the night sky.  The final map is shown above, but for those interested, a  full explanation with code snippets is shown.

When I moved to the London area from the US last year, I noticed that not only were the summer nights shorter than in New Jersey, but twilight seemed to go on forever.  London is at about 52° N, so for a few days around the summer solstice, the sky never reaches astronomical twilight (when the sun is 15° below the horizon, which is also when it is dark enough for imaging).

So I wondered how much imaging time I get each year for each area of the sky. (Ignoring weather of course.)  And beyond my view of the sky, what does the picture looks like when you consider the view from everywhere on the planet.  Obviously, the summer sky isn’t in darkness as long as the winter sky, but how big is the difference?  And to take the question to a larger level, what happens when you consider where people live on the earth?

I’m sure there is a way to solve this analytically for someone really talented in spherical trigonometry, but I set out to solve it numerically.  I had a lot of R code that I’d written to create the maps in The Astrophotography Sky Atlas, so how hard could it be?  It turns out that it was a good bit of effort for such an unremarkable question.  (At least 30 hours, including this write up.)

Let’s start with my assumptions.  A point in the sky is considered “image-able” when:

• It’s fully dark: the sun is at least 15° below the horizon
• The point is at least 20° above the horizon. Lower than that, and you’re looking through over three airmasses.

First I had to write and test the code in R to do all the positional astronomy calculations.  This turned out to non-trivial, but after a couple of evenings reading and testing, I had functions that returned comparable values to the online ephemerides.  The sunrise and sunset values I get differ by a few minutes from various online sources, so there are probably more precise equations available, but let’s call these good enough.  I’ve posted snippets of them here for critique and transparency. These are all of the helper functions I wrote.

SunCoords<-read.csv("SunRADec2017.csv") #get sun coordinates for any year, expressed in decimal degrees

CalcSunRiseSetTime <- function(lat, day, altitude, sunset=TRUE){ #calculates time sun is at altitude. Default is sunset=descending sun.
d<-day+6208.5 #6208.5 days from Jan 1 2000 noon to Jan 1 2017 0:00, which is the year I'm using here
ST0<-GMT2GMSThours(d)*15 #calculate sidereal time at long=0 at midnight, expressed in degrees

#use ifelse so function is vectorized...
#modulus handles discontinuity near equinox
EoT<-((SunCoords[day,]$RADegrees - ST0 - 0) / 15)%%24 #"Equation of time" at long=0, expressed in hours

#semi-diurnal arc
cosH<-(sin(altitude/57.2957795) - sin(lat/57.2957795)*sin(SunCoords[day,]$DecDegrees/57.2957795)) / (cos(lat/57.2957795) * cos(SunCoords[day,]$DecDegrees/57.2957795)) #handle cases near poles:&nbsp; if cosH >1, return ~11:59pm, as sun never rises.&nbsp; if <-1, return ~12:01am, as sun never sets
if(sunset){
r<-ifelse(cosH > 1, 0.001, ifelse(cosH < -1, 23.999, EoT + (acos(cosH)*57.2957795/15)))
}
else{
r<-ifelse(cosH > 1, 23.999, ifelse(cosH < -1, 0.001, EoT - (acos(cosH)*57.2957795/15)))
}
r #time in decimal hours
}

CalcObjectRiseSetTime <- function(lat, day, RA, Dec, altitude, set=TRUE){ #RA and Dec of object in degrees
cosLHA<-(sin(altitude/57.2957795) - sin(lat/57.2957795)*sin(Dec/57.2957795)) / (cos(lat/57.2957795) * cos(Dec/57.2957795)) if(set){ if(cosLHA > 1) return(0.001)
if(cosLHA < -1) return(23.999)
SiderealAnswer<-RA + (acos(cosLHA)*57.2957795) #in degrees } else{ if(cosLHA > 1) return(23.999)
if(cosLHA < -1) return(0.001)
SiderealAnswer<-RA - (acos(cosLHA)*57.2957795) #in degrees
}

#convert Greenwich sidereal time to GMT
GMST2GMT(day, SiderealAnswer/360*24)
}

GMST2GMT<- function(day, ST){#input day of 2017 and sidereal time in hours ***why is it up to 8 minutes off?
#return value in hours
((ST/24-0.277-day*0.00273042)*24)%%24 #2017 starts at 0.277 day sidereal time, declines by 3m56s each day
}

GMT2GMSThours<- function(GMTday){ #input days and fraction since J2000.0
(18.697374558 + 24.06570982441908 *GMTday)%%24 #return sidereal time at long=0, expressed in hours
}

TimeInNightSky <- function(lat, day, RA, Dec, MinObjectAltitude, SunAltitude){ #
ObjectRise<-CalcObjectRiseSetTime(lat, day, RA, Dec, MinObjectAltitude, FALSE)
ObjectSet<-CalcObjectRiseSetTime(lat, day, RA, Dec, MinObjectAltitude, TRUE)
SunRise<-CalcSunRiseSetTime(lat, day, SunAltitude, sunset=FALSE)
SunSet<-CalcSunRiseSetTime(lat, day, SunAltitude, sunset=TRUE) if(ObjectRise > ObjectSet){
time<-min(ObjectSet, SunRise) + 24 - max(ObjectRise, SunSet)
}
else{
time<-max(0, SunRise-ObjectRise) + max(0, ObjectSet-SunSet)
}
time
}

CalcEoT <- function(day){
d<-day+6208.5 #6208.5 days from Jan 1 2000 noon to Jan 1 2017 0:00, which is the year I'm using here
ST0<-GMT2GMSThours(d)*15 #calculate sidereal time at long=0 at midnight, expressed in degrees

#use ifelse so function is vectorized...
EoT<-((SunCoords[day,]$RADegrees - ST0 - 0) / 15)%%24 #"Equation of time" at long=0, expressed in hours
#EoT<-ifelse(abs(EoT) > 12, EoT+24, EoT) #handle discontinuity near vernal equinox
EoT
}

#calculate map data for one night at one location
lati<- 25 #latitude to use
minAlt<-20 #lowest imageable altitude
lowDec<-max(-89, lati-90+minAlt) #only compute values for Dec within range of imageability for latitude
highDec<-min(89, 90+lati-minAlt)
decSequence<-seq(lowDec,highDec,2)
RASequence<-seq(1,365,2)
timegrid=data.frame(matrix(0, nrow=180, ncol=365))

for (r in RASequence){
for(d in decSequence){
timegrid[d+90, r]<-timegrid[d+90, r] + TimeInNightSky(lati, 105, r, d, minAlt, -15) #-15 is position of sun at astronomical twilight
}
}

#now convert the dataframe above into a format readable by the plotting code
timemap2 <- data.frame(RA=numeric(), Dec=numeric(), time=numeric() )
for (x in RASequence){ #only look where values were loaded into the matrix .. old version was "1:ncol(timegrid)){"
for(y in decSequence+90){ #1:nrow(timegrid)){
temp<-data.frame(x, y-90, timegrid[y,x]) #the +/-90 is to handle the negative Dec values, which are invalid indexes in a matrix
cols <- c("RA","Dec","time")
colnames(temp) <- cols
timemap2<-rbind(timemap2, temp)
}
}
#flip the RA axis
timemap2$RA<-365-timemap2$RA

With these functions, R can easily plot a few interesting things to verify the code makes sense.  One is the Equation of Time.

Sunrise and sunsets (astronomical twilight times, actually) for a year look pretty good, too.  This is sunrise and sunset for every day of the year at 25° N, not accounting for Daylight Savings Time of course.

The functions appear to work well, so I iterated those functions to calculate a single night sky (April 15^th) from London, to see what parts get the most image-able hours.

#calculate map data for one night at one location
lati<- 25 #latitude to use
minAlt<-20 #lowest imageable altitude
lowDec<-max(-89, lati-90+minAlt) #only compute values for Dec within range of imageability for latitude
highDec<-min(89, 90+lati-minAlt)
decSequence<-seq(lowDec,highDec,2)
RASequence<-seq(1,365,2)
timegrid=data.frame(matrix(0, nrow=180, ncol=365))

for (r in RASequence){
for(d in decSequence){
timegrid[d+90, r]<-timegrid[d+90, r] + TimeInNightSky(lati, 105, r, d, minAlt, -15) #-15 is position of sun at astronomical twilight
}
}

#now convert the dataframe above into a format readable by the plotting code
timemap2 <- data.frame(RA=numeric(), Dec=numeric(), time=numeric() )
for (x in RASequence){ #only look where values were loaded into the matrix .. old version was "1:ncol(timegrid)){"
for(y in decSequence+90){ #1:nrow(timegrid)){
temp<-data.frame(x, y-90, timegrid[y,x]) #the +/-90 is to handle the negative Dec values, which are invalid indexes in a matrix
cols <- c("RA","Dec","time")
colnames(temp) <- cols
timemap2<-rbind(timemap2, temp)
}
}
#flip the RA axis
timemap2$RA<-365-timemap2$RA

Then I plotted the resulting data (using code not shown):

For comparison, I also plotted London’s sky at the winter and summer solstices:

You can see that in the winter, London provides plenty of dark sky time (though it’s rarely clear), but near the summer solstice the sun never gets 15° below the horizon!

For comparison, let’s calculate the same chart for Miami (at 25° N) on April 15^th to see how things change.  (Because the Earth is rotating, longitude is irrelevant to our calculations.)

You can see how the longer spring nights in Miami, compared to London, give you more dark sky time this time of year.  The opposite is true for winter, where the nights are looooong in London.

Now that I had one night, it was time to sum all of the nights for a year.  Even on a fairly fast computer, this took some time, so I had to go back and see if I could make the code more efficient.  R is an interpreted language, but highly optimized for vectors.  Unfortunately, this code required several “for” loops, which R does not handle quickly.  There is probably an elegant and fast solution with dplyr or something, but I couldn’t make it work. Using a resolution of 2×2° squares, it still took 40-80 minutes to crunch the numbers for each map, depending on latitude.

Here is the map from 50° N.

How do the maps look near the equator?  And what about near the poles?  Here are maps at 70°N, 20°N, and the equator.

You can see that the seasonal changes are greater the further you move away from the equator.

Now that everything was working, I generated full-sky data for latitudes of 80° N down to 80° S in 2.5° increments.  This took about 60 hours of computer time (mostly while I was sleeping or at work).  I’m sure there is a more elegant and faster way to do it in R, but again, I used for loops. This gave me the data I needed to calculate the “average” view from everywhere on Earth.

###avoid rbind by using a matrix, then re-shape the matrix back into a dataframe the plotting code can read
for(lati in seq(-80, 80, 2.5))
{
minAlt<-20 #lowest imageable altitude
lowDec<-max(-89, lati-90+minAlt) #only compute values for Dec within range of imageability for latitude
highDec<-min(89, 90+lati-minAlt)
decSequence<-seq(lowDec,highDec,2)
RASequence<-seq(1,365,2)
timegrid=data.frame(matrix(0, nrow=180, ncol=365))

for (r in RASequence){
for(d in decSequence){
for(date in 1:365){
timegrid[d+90, r]<-timegrid[d+90, r] + TimeInNightSky(lati, date, r, d, minAlt, -15) #-15 is position of sun at astronomical twilight
}
}
}
#get the average
timegrid<-timegrid/365

#covert the grid into a dataframe of the right format to be read by mapping code below
timemap2 <- data.frame(RA=numeric(), Dec=numeric(), time=numeric() )
for (x in RASequence){ #only look where values were loaded into the matrix .. old version was "1:ncol(timegrid)){"
for(y in decSequence+90){ #1:nrow(timegrid)){
temp<-data.frame(x, y-90, timegrid[y,x]) #the +/-90 is to handle the negative Dec values, which are invalid indexes in a matrix
cols <- c("RA","Dec","time")
colnames(temp) <- cols
timemap2<-rbind(timemap2, temp)
}
}
#flip the RA axis, since map is printed from 24 to 0 hrs
timemap2$RA<-365-timemap2$RA

#save as a separate file
fname=paste("NightMap Lat=", lati, ".Rda", sep="")
save(timemap2, file=fname)
}

Bringing all of my maps together, this is the view from everywhere on earth (technically, everywhere between latitudes 80° N and 80° S):

Two things stand out:

1. The contours are a little jagged. This is a result of running the calculations for latitudes at 2.5° increments.  Increasing the “latitude resolution” would further smooth them.
2. The hemispheres aren’t completely symmetrical. Summer in the southern hemisphere seems to have a slightly greater impact on the southern night sky than the northern hemisphere’s summer does on the northern sky.  At first I thought it could be an error in my equations (and it could indeed be).  But I think it’s the effect of Earth’s elliptical orbit.  The changes in sunlight on the days near the solstices are not symmetrical (see the shape of the analemma for proof).  Someone more skilled in positional astronomy could surely offer a clearer explanation or show me I’m wrong.

All in all, I think it’s a pretty neat map.  The winter sky clearly gets more time in darkness, and the southern hemisphere has more extreme differences between the summer and winter. Surprisingly, the area of sky with the least good imaging time on average is just south of the most-imaged object in the night sky:  M42 in Orion.  And the area with most imaging time available on average is the stinger of Scorpius, where there are a slew of beautiful objects.

But we can take this even further.  You may object that people are not evenly distributed by latitude.  And indeed, that is why you see so many images of M42 compared to the Cat’s Paw Nebula, which gets more hours in darkness.  A proper map of “image-able” time for the night sky would incorporate where people actually live by latitude.  So let’s account for where people live.  First, I had to get the data, which had to be distilled from a much larger data set.

The data on population from SEDAC* come in a raster format, so it took a little more code to get a simple percentage of population at each latitude.  (Note that in the code here, I cut the data slightly differently, so the bins would be centered on 5° increments. )

library(rgdal)
library(sp)
library(raster)
GPW<-raster("gpw-v4-population-count-adjusted-to-2015-unwpp-country-totals_2015.tif")
popPerLatDegree<-aggregate(GPW, c(43200,1200), fun=sum) #aggregate data across longitudes and into 2.5 degree groups by latitude
popPerLatDegree<-as.data.frame(popPerLatDegree, xy=TRUE)
names(popPerLatDegree)[3] <- "pop" #shorten column name
popPerLatDegree$x<-NULL #we don't need the longitude, since there's only one value now
popPerLatDegree[is.na(popPerLatDegree$pop),]$pop<-0 #replace the NA's with zeros
popPerLatDegreeTrimmed<-popPerLatDegree[3:142,] #trim top two degrees off so center of each group of five degrees (approximately) aligns with 0,5,10,etc.

#sum every five degrees into one value
popPerFiveDegree <- data.frame(centralLat=numeric(), pop=numeric() )
for(n in seq(1, nrow(popPerLatDegreeTrimmed), 5))
{
total<-0
for(y in 0:4)
{
total<-total + popPerLatDegreeTrimmed[n+y,]$pop
}
tmp<-data.frame(popPerLatDegreeTrimmed[n+2,]$y, total)
cols <- c("centralLat","pop")
colnames(tmp) <- cols
popPerFiveDegree<-rbind(popPerFiveDegree, tmp)
}

#calculate proportion in each band
worldPop<-sum(popPerFiveDegree$pop)
popPerFiveDegree$pop<-popPerFiveDegree$pop/worldPop

The following table shows the percent of the world’s population by latitude:

 

A couple of things are immediately apparent:

• Only 13% of the world’s population lives in the southern hemisphere.
• 49% of the world’s people live in the range from 20 to 40° N.

Okay, now that we know where people live, we can re-create the sky map weighted by where people live.

Sadly, it turned out to be pretty uninteresting. The results are simpler than I first thought because the distribution of human population is so skewed that its effect swamps any seasonality.  The final plot mostly consists of data from the mid-northern latitudes, with only a tiny contribution from other latitudes.   The further south you go in declination, the fewer people there are to see it.  The short answer to my original question is that the most and least image-able parts of the sky is basically, “near the north and south celestial poles, respectively.”

* Population data from:  Center for International Earth Science Information Network – CIESIN – Columbia University. 2016. Gridded Population of the World, Version 4 (GPWv4): Population Count Adjusted to Match 2015 Revision of UN WPP Country Totals. Palisades, NY: NASA Socioeconomic Data and Applications Center (SEDAC). http://dx.doi.org/10.7927/H4SF2T42. Accessed 20 April 2017.

The Insight Astronomy Photographer of the Year Exhibition

It’s rare that you see much about astrophotography in the popular press or in museums, but this weekend I had the chance to see the exhibit for The Insight Astronomy Photographer of the Year contest. This contest is in its fifth year, drawing entries from around the world. The winners are not only published in a beautiful yearbook, but they are also exhibited at the Royal Observatory in Greenwich, UK. The exhibit is part of the free area, though if you are into astronomy, you may want to pay the fee to see the nearby observatory museum, where luminaries like Airy and Flamsteed worked. (Of course, you also get to stand on the prime meridian line, which is pretty cool too.)

As you can see the exhibit was popular with the public, which is heartening to see.

The winning images were printed and displayed on lightboxes, which really made them look fantastic. They also had a few video clips with interviews of the photographers playing. A few close-ups are below.

Overall, there seemed to be a bias toward nightscape images over deep-sky images, which makes sense, as it gives people a sense of their place in the universe, but I would personally like to have seen more. Two deep-sky images (not shown above), by Pavel Pech and Rolf Wahl Olsen,  were particularly impressive.

Each image listed the technical details and explained what you were looking at. Like the book, the exhibit quality was very high, and I know I came away inspired to get out with my camera more.

Lecture: Understanding Signal, Noise, and Resolution

Last night, I was honored to present a lecture to the New Jersey Astronomical Association’s imaging meeting. It’s a great group of people, and I’m sorry that I’m only connecting with them now.  Special thanks to Mike Franzyshen and Jim Roselli for the invitation.

In about an hour, I cover broad range of related topics:

• The statistics of shot noise, including the connections between Simeon-Denis Poisson and the particle theory of light
• Signal-to-noise ratio, with examples
• The effects of skyglow, with examples
• Resolution and sampling
• Aperture, focal length, and focal ratio

It was a lot of fun to put this lecture together, and during the research for it, I uncovered how Poisson’s life led to multiple insights about photons, none of which he was able to appreciate during his lifetime.

NOTE: At 40:30, I second-guessed myself, but the math on the slide is correct as written. The additional read noise is indeed 32. Since read noise is 10 per exposure, per the previous slide the additional noise is: SQRT(10^2 + 10^2 + 10^2 + 10^2 + 10^2 + 10^2 + 10^2 + 10^2 + 10^2 + 10^2) = 32.

Barnard’s Loop

When I was a kid, Barnard’s Loop was something that I saw on star charts, but it seemed so hopelessly dim, I never expected to actually see it.  And even when I started CCD imaging, it was still a somewhat elusive object: too large to capture unless you used a wide-angle lens, and even then you wouldn’t get decent resolution.  But the combination of a full-frame sensor and a very fast telephoto lens turns out to frame it nicely.

This image obviously has more in it than Barnard’s Loop.  M42/43, the Flame Nebula, the Horsehead Nebula, and M78 all sit nestled within the Loop.  But more interestingly for me, you can start to see the overall Orion Molecular Cloud complex in there: all the dim tendrils that connect each of these objects, some glowing, some blocking the view of the glow.  I regret stopping the lens down to f/2.8 now, as perhaps I would have captured more of the overall cloud that way.  I’d go back and retake the shot if I weren’t having so much fun with this new lens on other targets (and if I hadn’t spent five hours processing this one).  But I’ll consider this a success, as it’s another childhood dream accomplished.

(This image is reduced to 25% of full size, as the 6D’s output is over 20 megapixels.)

Image data:

• Exposures: 81×2 min at ISO800 – total exposure time:  2h 42m
• Telescope: Samyang 135 mm f/2 lens at f/2.8 (reviewed here)
• Camera: Canon 6D (modified) with Astronomik CLS clip-in filter
• Mount: Takahashi EM200
• Guiding: Orion Starshoot, guided using PHD2
• Conditions:  fair transparency, calm winds
• Processing: DeepSkyStacker -> PixInsight -> Photoshop
• Date: Feb 28, 2016

 

Review: Samyang 135 mm f/2 Astrograph

By calling the Samyang 135 mm f/2 lens an astrograph in the title, I’m giving away a bit of the conclusion, so let me just state the conclusions up front: I found this lens to be sharp from corner to corner with a reasonably flat field across a full-frame sensor. This is better performance than all but the best prime telephoto lenses, and also better than many telescopes I’ve owned claiming to be astrographs. Even more impressive is the fact that it accomplishes this at a focal ratio four to eight times faster than “fast” refractors. Finally, it’s hard to beat the price: currently just over $500.

While the performance was adequate at f/2, I found that stopping down the lens one full stop to f/2.8 improved sharpness. This is true of any lens, and even stopped to f/2.8, that’s still four times the light gathered per sensor area than an f/5.6 telescope. It’s easy to forget how fast this is, but my first night using the lens reminded me. My usual telescopes are two f/5 William Optics Star 71s and a Takahashi FSQ-106ED, also f/5. I usually shoot narrowband exposures of 20 minutes with these. So the combination of broadband and f/2.8 put me in the realm of 30 to 120 second exposures—anything longer at ISO 800 overexposed the stars.

The infinity focus point is about 2 mm left of the mark on the lens barrel, so you’ll have to carefully dial in focus. There isn’t much tolerance for error at such fast focal ratios, as the zone of focus is very narrow.

Vignetting is substantial in the corners, but it is more reasonable if you move slightly inward.  For very fast optics, this is typical, though it does lead to lower SNR toward the corners. I was able to keep the full frame images without cropping by using good flat frames, but this is essential.

Master Flat created from 33 pictures (Average)

The quick 99×1 minute image below of the Rosette Nebula area gives you a sense of how wide the view is with a Canon 6D. In the center is the Christmas Tree/Foxfur/Cone Nebula area, with huge dark nebula Barnard 37 prominent. This was taken without a CLS filter, so light pollution prevented me from adequately revealing the Foxfur nebula well. The Rosette Nebula shines brightly to the left, though. (Note that this image is reduced to 25% of the actual image resolution.)

As you can see from the full resolution close-ups below, the lens is impressively sharp across the Canon 6D’s entire field of view, with very minimal distortion even in the extreme corners.

The 9-blade diaphragm of the Samyang results in a pleasing radiant around bright stars, but the lens does exhibit some internal reflections.

I look forward to using this lens as my (very) widefield astrograph. Depending on your sensor size, the ideal targets for this lens will vary, but I’m looking forward to shooting:

• The Orion Molecular Cloud Complex
• IC2177 and Thor’s Helmet area
• The California Nebula to Pleiades area
• Orion’s Head/ Meissa Nebula
• Heart and Soul Nebulae area
• The Rho Ophiuchi area
• Cygnus
• The IC405/IC410 area
• Taurus Molecular Cloud
• The Milky Way’s Pipe nebula region
• Sagittarius

After complex mosaics and multi-night narrowband CCD projects, it’s a joy to throw a simple setup like this onto the mount to grab bright widefield images in a few hours.

Pros:

• Fast focal ratio
• Sharpness
• Flatness of field
• Price (currently ~$529 USD)

Cons:

• Some internal reflections
• Limited targets available for this focal length
• Will require adapter to fit CCD cameras

Orion Rising Over Courthouse Butte

Sedona, Arizona has the clearest, darkest skies I may ever see (anywhere near civilization at least).  Who knew there was a winter Milky Way visible too?!  Not this suburbanite.

I got one clear night to test out both my new Samyang 14mm lens and the iOptron SkyTracker.  While this image is not the best example of the SkyTracker’s abilities, since I misaligned it, the foreground framing was better than the other shots I took.  So this is the one I chose to process first.

Image data:

• Exposures:  sky:  20×30 seconds, foreground: 1×30 seconds
• Telescope: Samyang 14mm f/2.8 lens at f/4
• Cameras: Canon 6D (modified)
• Mount: iOptron SkyTracker
• Guiding: none
• Conditions:  excellent transparency, passing clouds
• Processing: DeepSkyStacker -> PixInsight -> Photoshop
• Date: Jan 29, 2016