function solstice,yyyy,lat

;------------------------------------------------------------------------------------------------
; Routine computes the day of year for summer solstice in either hemisphere. 
;
; Accounts for leap year oscillation with a fit to data from:
; http://www.usno.navy.mil/USNO/astronomical-applications/data-services/earth-seasons.
;
; This routine does not capture the slow change over the years as seen on:
; http://scienceworld.wolfram.com/astronomy/SummerSolstice.html
; 
; Input:
;   yyyy........year (e.g., 2007), scalar
;   lat.........1 = northern hemisphere, -1 = southern hemisphere, scalar
;
; Output:
;   sols........day of year of summer solstice   
;
; Source: Mark Hervig
;------------------------------------------------------------------------------------------------

;- linear fit coeffs. to the Navy data, solstice (as day of year) vs. years since leap year, 
;  see "solstice_look.pro"

   rf_n = [171.995, 0.246490]  ; for north
   rf_s = [355.487, 0.249108]  ; for south
   
;- Find the solstice
   
   sols = -999.
     
   leapyr = yyyy - fix(yyyy/4) * 4  ; leapyr = 0 if it is a leap year
     
   if (leapyr eq 0) then leapyr = 4
     
   if (lat gt 0) then sols = poly(leapyr,rf_n)  ; Northen Hemisphere

   if (lat lt 0) then sols = poly(leapyr,rf_s)  ; Southern Hemisphere
          
;- done

   return, sols
   end