Position of the sun given time of day, and lat/long

Question

Anyone want to share some nice code that gives the position of the sun (elevation and azimuth) given latitude and longitude and time of day ?

And date, of course.

Answer 1

You'll not only need the time of day, but also the date. The position of the sun depends on the time of year.

Answer 2

How precise an answer do you need?

(For example, are your latitude and longitude expressed in WGS84 or some other coordinate system?)

P. J. Naughter's code probably has everything you need and then some...

http://www.naughter.com/aa.html

Answer 3

"Practical Astronomy with your Calculator" by Peter Duffett-Smith provides well-illustrated explanations and equations for this sort of thing. From these, its easy to code solutions in whatever language you choose.

Answer 4

PHP actually has some built-in functions that deal with similar tasks: http://au2.php.net/manual/en/function.date-sun-info.php

Unfortunately, this only tells you the time of sun rise, sun set, twilight, and some other things: not the location of the sun at a particular time.

I did find some formulae for this sort of thing here:

• http://www.qarlos.free.fr/navegacion/Sextante/formulas.htm (FYI: Firefox reports this as a dangerous site)
• http://www.kevinboone.com/suntimes.html
• http://williams.best.vwh.net/sunrise_sunset_algorithm.htm

Perhaps you could adapt them to your purposes..?

Answer 5

Doing this correctly isn't easy. It depends on the orbit of the Earth, which is probably beyond what you were looking to do. The easiest way to get a completely correct answer is to pull it from JPL Horizons. Just get the az/el readings once an hour, and linear interpolation. That'll look pretty good visually.

If you need to figure it out for a lot of positions on the Earth, get Horizons to give you ECI coordinates (which you can again linearly interpolate between), and then convert the ECI coordinates to az/el. See Celestrak for information on that conversion. You want the "Orbital Coordinate Systems" articles. Astronomical Algorithms by Meeus is also useful.

Answer 6

The National Renewable Energy Lab has a Solar Position and Intensity Calculator called Solpos. I used it years ago to track sun position to avoid "sunstrikes" against our optical communication systems.

See the C source code at: http://rredc.nrel.gov/solar/codesandalgorithms/solpos/

You can test drive the algorithm at: http://www.nrel.gov/midc/solpos/solpos.html

If you need real precise data, fed from NASA-JPL, check out JPL's HORIZONS data service which has various gateways to the data used by NASA to track celestial bodies. They have "431015 asteroids, 2966 comets, 168 planetary satellites, 8 planets, the Sun, L1, L2, select spacecraft, and system barycenters". HORIZONS is at http://ssd.jpl.nasa.gov/?horizons

Answer 7

Here is my implementation in R of the Astronomer's Almanac algorithm (Joseph J. Michalsky. The astronomical almanac’s algorithm for approximate solar position (1950–2050). Solar Energy, 40(3):227–235, 1988.) It's straightforward to translate it in any other language:

sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
                        lat=46.5, long=6.5) {
  twopi <- 2 * pi
  deg2rad <- pi / 180

  # Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
  month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
  day <- day + cumsum(month.days)[month]
  leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
  day[leapdays] <- day[leapdays] + 1

  # Get Julian date - 2400000
  hour <- hour + min / 60 + sec / 3600 # hour plus fraction
  delta <- year - 1949
  leap <- trunc(delta / 4) # former leapyears
  jd <- 32916.5 + delta * 365 + leap + day + hour / 24

  # The input to the Atronomer's almanach is the difference between
  # the Julian date and JD 2451545.0 (noon, 1 January 2000)
  time <- jd - 51545.

  # Ecliptic coordinates

  # Mean longitude
  mnlong <- 280.460 + .9856474 * time
  mnlong <- mnlong %% 360
  mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360

  # Mean anomaly
  mnanom <- 357.528 + .9856003 * time
  mnanom <- mnanom %% 360
  mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
  mnanom <- mnanom * deg2rad

  # Ecliptic longitude and obliquity of ecliptic
  eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
  eclong <- eclong %% 360
  eclong[eclong < 0] <- eclong[eclong < 0] + 360
  oblqec <- 23.429 - 0.0000004 * time
  eclong <- eclong * deg2rad
  oblqec <- oblqec * deg2rad

  # Celestial coordinates
  # Right ascension and declination
  num <- cos(oblqec) * sin(eclong)
  den <- cos(eclong)
  ra <- atan(num / den)
  ra[den < 0] <- ra[den < 0] + pi
  ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
  dec <- asin(sin(oblqec) * sin(eclong))

  # Local coordinates
  # Greenwich mean sidereal time
  gmst <- 6.697375 + .0657098242 * time + hour
  gmst <- gmst %% 24
  gmst[gmst < 0] <- gmst[gmst < 0] + 24.

  # Local mean sidereal time
  lmst <- gmst + long / 15.
  lmst <- lmst %% 24.
  lmst[lmst < 0] <- lmst[lmst < 0] + 24.
  lmst <- lmst * 15. * deg2rad

  # Hour angle
  ha <- lmst - ra
  ha[ha < -pi] <- ha[ha < -pi] + twopi
  ha[ha > pi] <- ha[ha > pi] - twopi

  # Latitude to radians
  lat <- lat * deg2rad

  # Azimuth and elevation
  el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
  az <- asin(-cos(dec) * sin(ha) / cos(el))
  elc <- asin(sin(dec) / sin(lat))
  az[el >= elc] <- pi - az[el >= elc]
  az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi

  el <- el / deg2rad
  az <- az / deg2rad
  lat <- lat / deg2rad

  return(list(elevation=el, azimuth=az))
}

Answer 8

If you can get access to IDL (I think trial versions are free) and the astronomy user's library (which is free) questions like this are very easy. Take a look at sunpos.

Answer 9

Here's code that uses ephem a python package for performing high-precision astronomy computations:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import datetime
import ephem # to install, run `pip install pyephem`

o = ephem.Observer()
# Los Angeles, Calif. 34°3'N, 118°15'W
o.lat, o.long, o.date = '34:3', '-118:15', datetime.datetime.utcnow()
sun = ephem.Sun(o)

print o.date
print sun.az, sun.alt    

Output

2010/3/29 20:52:25 UTC
205:40:44.9 56:59:27.7

, where:

• az — azimuth east of north
• alt — altitude above horizon