title: Position of the sun given time of day, latitude and longitude
body: <p>This question has been asked <a href="https://stackoverflow.com/questions/257717/position-of-the-sun-given-time-of-day-and-lat-long/258106#comment10668948_258106">before</a> a little over three years ago. There was an answer given, however I've found a glitch in the solution.</p>
  
  <p>Code below is in R. I've ported it to another language, however have tested the original code directly in R to ensure the issue wasn't with my porting.</p>
  
  <pre><code>sunPosition &lt;- function(year, month, day, hour=12, min=0, sec=0,
                      lat=46.5, long=6.5) {
  
  
    twopi &lt;- 2 * pi
    deg2rad &lt;- pi / 180
  
    # Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
    month.days &lt;- c(0,31,28,31,30,31,30,31,31,30,31,30)
    day &lt;- day + cumsum(month.days)[month]
    leapdays &lt;- year %% 4 == 0 &amp; (year %% 400 == 0 | year %% 100 != 0) &amp; day &gt;= 60
    day[leapdays] &lt;- day[leapdays] + 1
  
    # Get Julian date - 2400000
    hour &lt;- hour + min / 60 + sec / 3600 # hour plus fraction
    delta &lt;- year - 1949
    leap &lt;- trunc(delta / 4) # former leapyears
    jd &lt;- 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 &lt;- jd - 51545.
  
    # Ecliptic coordinates
  
    # Mean longitude
    mnlong &lt;- 280.460 + .9856474 * time
    mnlong &lt;- mnlong %% 360
    mnlong[mnlong &lt; 0] &lt;- mnlong[mnlong &lt; 0] + 360
  
    # Mean anomaly
    mnanom &lt;- 357.528 + .9856003 * time
    mnanom &lt;- mnanom %% 360
    mnanom[mnanom &lt; 0] &lt;- mnanom[mnanom &lt; 0] + 360
    mnanom &lt;- mnanom * deg2rad
  
    # Ecliptic longitude and obliquity of ecliptic
    eclong &lt;- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
    eclong &lt;- eclong %% 360
    eclong[eclong &lt; 0] &lt;- eclong[eclong &lt; 0] + 360
    oblqec &lt;- 23.429 - 0.0000004 * time
    eclong &lt;- eclong * deg2rad
    oblqec &lt;- oblqec * deg2rad
  
    # Celestial coordinates
    # Right ascension and declination
    num &lt;- cos(oblqec) * sin(eclong)
    den &lt;- cos(eclong)
    ra &lt;- atan(num / den)
    ra[den &lt; 0] &lt;- ra[den &lt; 0] + pi
    ra[den &gt;= 0 &amp; num &lt; 0] &lt;- ra[den &gt;= 0 &amp; num &lt; 0] + twopi
    dec &lt;- asin(sin(oblqec) * sin(eclong))
  
    # Local coordinates
    # Greenwich mean sidereal time
    gmst &lt;- 6.697375 + .0657098242 * time + hour
    gmst &lt;- gmst %% 24
    gmst[gmst &lt; 0] &lt;- gmst[gmst &lt; 0] + 24.
  
    # Local mean sidereal time
    lmst &lt;- gmst + long / 15.
    lmst &lt;- lmst %% 24.
    lmst[lmst &lt; 0] &lt;- lmst[lmst &lt; 0] + 24.
    lmst &lt;- lmst * 15. * deg2rad
  
    # Hour angle
    ha &lt;- lmst - ra
    ha[ha &lt; -pi] &lt;- ha[ha &lt; -pi] + twopi
    ha[ha &gt; pi] &lt;- ha[ha &gt; pi] - twopi
  
    # Latitude to radians
    lat &lt;- lat * deg2rad
  
    # Azimuth and elevation
    el &lt;- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
    az &lt;- asin(-cos(dec) * sin(ha) / cos(el))
    elc &lt;- asin(sin(dec) / sin(lat))
    az[el &gt;= elc] &lt;- pi - az[el &gt;= elc]
    az[el &lt;= elc &amp; ha &gt; 0] &lt;- az[el &lt;= elc &amp; ha &gt; 0] + twopi
  
    el &lt;- el / deg2rad
    az &lt;- az / deg2rad
    lat &lt;- lat / deg2rad
  
    return(list(elevation=el, azimuth=az))
  }
  </code></pre>
  
  <p>The problem I'm hitting is that the azimuth it returns seems wrong. For example, if I run the function on the (southern) summer solstice at 12:00 for locations 0ºE and 41ºS, 3ºS, 3ºN and 41ºN:</p>
  
  <pre><code>&gt; sunPosition(2012,12,22,12,0,0,-41,0)
  $elevation
  [1] 72.42113
  
  $azimuth
  [1] 180.9211
  
  &gt; sunPosition(2012,12,22,12,0,0,-3,0)
  $elevation
  [1] 69.57493
  
  $azimuth
  [1] -0.79713
  
  Warning message:
  In asin(sin(dec)/sin(lat)) : NaNs produced
  &gt; sunPosition(2012,12,22,12,0,0,3,0)
  $elevation
  [1] 63.57538
  
  $azimuth
  [1] -0.6250971
  
  Warning message:
  In asin(sin(dec)/sin(lat)) : NaNs produced
  &gt; sunPosition(2012,12,22,12,0,0,41,0)
  $elevation
  [1] 25.57642
  
  $azimuth
  [1] 180.3084
  </code></pre>
  
  <p>These numbers just don't seem right. The elevation I'm happy with - the first two should be roughly the same, the third a touch lower, and the fourth much lower. However the first azimuth should be roughly due North, whereas the number it gives is the complete opposite. The remaining three should point roughly due South, however only the last one does. The two in the middle point just off North, again 180º out.</p>
  
  <p>As you can see there are also a couple of errors triggered with the low latitudes (close the equator)</p>
  
  <p>I believe the fault is in this section, with the error being triggered at the third line (starting with <code>elc</code>).</p>
  
  <pre><code>  # Azimuth and elevation
    el &lt;- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
    az &lt;- asin(-cos(dec) * sin(ha) / cos(el))
    elc &lt;- asin(sin(dec) / sin(lat))
    az[el &gt;= elc] &lt;- pi - az[el &gt;= elc]
    az[el &lt;= elc &amp; ha &gt; 0] &lt;- az[el &lt;= elc &amp; ha &gt; 0] + twopi
  </code></pre>
  
  <p>I googled around and found a similar chunk of code in C, converted to R the line it uses to calculate the azimuth would be something like</p>
  
  <pre><code>az &lt;- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))
  </code></pre>
  
  <p>The output here seems to be heading in the right direction, but I just can't get it to give me the right answer all the time when it's converted back to degrees.</p>
  
  <p>A correction of the code (suspect it's just the few lines above) to make it calculate the correct azimuth would be fantastic.</p>