+ "List of *local* times of sunrise, sunset, and daylight on Gregorian DATE.
+
+Corresponding value is nil if there is no sunrise/sunset."
+ (let* (; first, get the exact moment of local noon.
+ (exact-local-noon (solar-exact-local-noon date))
+ ; get the time from the 2000 epoch.
+ (t0 (solar-julian-ut-centuries (car exact-local-noon)))
+ ; store the sidereal time at Greenwich at midnight of UT time.
+ ; find if summer or winter slightly above the equator
+ (equator-rise-set
+ (progn (setq solar-sidereal-time-greenwich-midnight
+ (solar-sidereal-time t0))
+ (solar-sunrise-and-sunset
+ (list t0 (car (cdr exact-local-noon)))
+ 1.0
+ (calendar-longitude) 0)))
+ ; store the spring/summer information,
+ ; compute sunrise and sunset (two first components of rise-set).
+ ; length of day is the third component (it is only the difference
+ ; between sunset and sunrise when there is a sunset and a sunrise)
+ (rise-set
+ (progn
+ (setq solar-northern-spring-or-summer-season
+ (if (> (car (cdr (cdr equator-rise-set))) 12) t nil))
+ (solar-sunrise-and-sunset
+ (list t0 (car (cdr exact-local-noon)))
+ (calendar-latitude)
+ (calendar-longitude) -0.61)))
+ (rise (car rise-set))
+ (adj-rise (if rise (dst-adjust-time date rise) nil))
+ (set (car (cdr rise-set)))
+ (adj-set (if set (dst-adjust-time date set) nil))
+ (length (car (cdr (cdr rise-set)))) )
+ (list
+ (and rise (calendar-date-equal date (car adj-rise)) (cdr adj-rise))
+ (and set (calendar-date-equal date (car adj-set)) (cdr adj-set))
+ (solar-daylight length))))
+
+(defun solar-sunrise-sunset-string (date)
+ "String of *local* times of sunrise, sunset, and daylight on Gregorian DATE."
+ (let ((l (solar-sunrise-sunset date)))
+ (format
+ "%s, %s at %s (%s hours daylight)"
+ (if (car l)
+ (concat "Sunrise " (apply 'solar-time-string (car l)))
+ "No sunrise")
+ (if (car (cdr l))
+ (concat "sunset " (apply 'solar-time-string (car (cdr l))))
+ "no sunset")
+ (eval calendar-location-name)
+ (car (cdr (cdr l))))))
+
+(defun solar-julian-ut-centuries (date)
+ "Number of Julian centuries elapsed since 1 Jan, 2000 at noon U.T. for Gregorian DATE."
+ (/ (- (calendar-absolute-from-gregorian date)
+ (calendar-absolute-from-gregorian '(1 1.5 2000)))
+ 36525.0))
+
+(defun solar-ephemeris-time(time)
+ "Ephemeris Time at moment TIME.
+
+TIME is a pair with the first component being the number of Julian centuries
+elapsed at 0 Universal Time, and the second component being the universal
+time. For instance, the pair corresponding to November 28, 1995 at 16 UT is
+\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
+Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT.
+
+Result is in julian centuries of ephemeris time."
+ (let* ((t0 (car time))
+ (ut (car (cdr time)))
+ (t1 (+ t0 (/ (/ ut 24.0) 36525)))
+ (y (+ 2000 (* 100 t1)))
+ (dt (* 86400 (solar-ephemeris-correction (floor y)))))
+ (+ t1 (/ (/ dt 86400) 36525))))
+
+(defun solar-date-next-longitude (d l)
+ "First moment on or after Julian day number D when sun's longitude is a
+multiple of L degrees at calendar-location-name with that location's
+local time (including any daylight savings rules).
+
+L must be an integer divisor of 360.
+
+Result is in local time expressed astronomical (Julian) day numbers.
+
+The values of calendar-daylight-savings-starts,
+calendar-daylight-savings-starts-time, calendar-daylight-savings-ends,
+calendar-daylight-savings-ends-time, calendar-daylight-time-offset, and
+calendar-time-zone are used to interpret local time."
+ (let* ((long)
+ (start d)
+ (start-long (solar-longitude d))
+ (next (mod (* l (1+ (floor (/ start-long l)))) 360))
+ (end (+ d (* (/ l 360.0) 400)))
+ (end-long (solar-longitude end)))
+ (while ;; bisection search for nearest minute
+ (< 0.00001 (- end start))
+ ;; start <= d < end
+ ;; start-long <= next < end-long when next != 0
+ ;; when next = 0, we look for the discontinuity (start-long is near 360
+ ;; and end-long is small (less than l).
+ (setq d (/ (+ start end) 2.0))
+ (setq long (solar-longitude d))
+ (if (or (and (/= next 0) (< long next))
+ (and (= next 0) (< l long)))
+ (progn
+ (setq start d)
+ (setq start-long long))
+ (setq end d)
+ (setq end-long long)))
+ (/ (+ start end) 2.0)))
+
+(defun solar-horizontal-coordinates
+ (time latitude longitude for-sunrise-sunset)
+ "Azimuth and height of the sun at TIME, LATITUDE, and LONGITUDE.
+
+TIME is a pair with the first component being the number of Julian centuries
+elapsed at 0 Universal Time, and the second component being the universal
+time. For instance, the pair corresponding to November 28, 1995 at 16 UT is
+\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
+Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT.
+
+The azimuth is given in degrees as well as the height (between -180 and 180)."
+ (let* ((ut (car (cdr time)))
+ (ec (solar-equatorial-coordinates time for-sunrise-sunset))
+ (st (+ solar-sidereal-time-greenwich-midnight
+ (* ut 1.00273790935)))
+ (ah (- (* st 15) (* 15 (car ec)) (* -1 (calendar-longitude))))
+ ; hour angle (in degrees)
+ (de (car (cdr ec)))
+ (azimuth (solar-atn2 (- (* (solar-cosine-degrees ah)
+ (solar-sin-degrees latitude))
+ (* (solar-tangent-degrees de)
+ (solar-cosine-degrees latitude)))
+ (solar-sin-degrees ah)))
+ (height (solar-arcsin
+ (+ (* (solar-sin-degrees latitude) (solar-sin-degrees de))
+ (* (solar-cosine-degrees latitude)
+ (solar-cosine-degrees de)
+ (solar-cosine-degrees ah))))))
+ (if (> height 180) (setq height (- height 360)))
+ (list azimuth height)))
+
+(defun solar-equatorial-coordinates (time for-sunrise-sunset)
+ "Right ascension (in hours) and declination (in degrees) of the sun at TIME.
+
+TIME is a pair with the first component being the number of Julian centuries
+elapsed at 0 Universal Time, and the second component being the universal
+time. For instance, the pair corresponding to November 28, 1995 at 16 UT is
+\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
+Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT."
+ (let* ((tm (solar-ephemeris-time time))
+ (ec (solar-ecliptic-coordinates tm for-sunrise-sunset)))
+ (list (solar-right-ascension (car ec) (car (cdr ec)))
+ (solar-declination (car ec) (car (cdr ec))))))
+
+(defun solar-ecliptic-coordinates (time for-sunrise-sunset)
+ "Apparent longitude of the sun, ecliptic inclination, (both in degrees)
+equation of time (in hours) and nutation in longitude (in seconds)
+at moment `time', expressed in julian centuries of Ephemeris Time
+since January 1st, 2000, at 12 ET."
+ (let* ((l (+ 280.46645
+ (* 36000.76983 time)
+ (* 0.0003032 time time))) ; sun mean longitude
+ (ml (+ 218.3165
+ (* 481267.8813 time))) ; moon mean longitude
+ (m (+ 357.52910
+ (* 35999.05030 time)
+ (* -0.0001559 time time)
+ (* -0.00000048 time time time))) ; sun mean anomaly
+ (i (+ 23.43929111 (* -0.013004167 time)
+ (* -0.00000016389 time time)
+ (* 0.0000005036 time time time))); mean inclination
+ (c (+ (* (+ 1.914600
+ (* -0.004817 time)
+ (* -0.000014 time time))
+ (solar-sin-degrees m))
+ (* (+ 0.019993 (* -0.000101 time))
+ (solar-sin-degrees (* 2 m)))
+ (* 0.000290
+ (solar-sin-degrees (* 3 m))))) ; center equation
+ (L (+ l c)) ; total longitude
+ (omega (+ 125.04
+ (* -1934.136 time))) ; longitude of moon's ascending node
+ ; on the ecliptic
+ (nut (if (not for-sunrise-sunset)
+ (+ (* -17.20 (solar-sin-degrees omega))
+ (* -1.32 (solar-sin-degrees (* 2 l)))
+ (* -0.23 (solar-sin-degrees (* 2 ml)))
+ (* 0.21 (solar-sin-degrees (* 2 omega))))
+ nil))
+ ; nut = nutation in longitude, measured in seconds of angle.
+ (ecc (if (not for-sunrise-sunset)
+ (+ 0.016708617
+ (* -0.000042037 time)
+ (* -0.0000001236 time time)) ; eccentricity of earth's orbit
+ nil))
+ (app (+ L
+ -0.00569
+ (* -0.00478
+ (solar-sin-degrees omega)))) ; apparent longitude of sun
+ (y (if (not for-sunrise-sunset)
+ (* (solar-tangent-degrees (/ i 2))
+ (solar-tangent-degrees (/ i 2)))
+ nil))
+ (time-eq (if (not for-sunrise-sunset)
+ (/ (* 12 (+ (* y (solar-sin-degrees (* 2 l)))
+ (* -2 ecc (solar-sin-degrees m))
+ (* 4 ecc y (solar-sin-degrees m)
+ (solar-cosine-degrees (* 2 l)))
+ (* -0.5 y y (solar-sin-degrees (* 4 l)))
+ (* -1.25 ecc ecc (solar-sin-degrees (* 2 m)))))
+ 3.1415926535)
+ nil)))
+ ; equation of time, in hours
+ (list app i time-eq nut)))
+
+(defconst solar-data-list
+ '((403406 4.721964 1.621043)
+ (195207 5.937458 62830.348067)
+ (119433 1.115589 62830.821524)
+ (112392 5.781616 62829.634302)
+ (3891 5.5474 125660.5691)
+ (2819 1.5120 125660.984)
+ (1721 4.1897 62832.4766)
+ (0 1.163 0.813)
+ (660 5.415 125659.31)
+ (350 4.315 57533.85)
+ (334 4.553 -33.931)
+ (314 5.198 777137.715)
+ (268 5.989 78604.191)
+ (242 2.911 5.412)
+ (234 1.423 39302.098)
+ (158 0.061 -34.861)
+ (132 2.317 115067.698)
+ (129 3.193 15774.337)
+ (114 2.828 5296.670)
+ (99 0.52 58849.27)
+ (93 4.65 5296.11)
+ (86 4.35 -3980.70)
+ (78 2.75 52237.69)
+ (72 4.50 55076.47)
+ (68 3.23 261.08)
+ (64 1.22 15773.85)
+ (46 0.14 188491.03)
+ (38 3.44 -7756.55)
+ (37 4.37 264.89)
+ (32 1.14 117906.27)
+ (29 2.84 55075.75)
+ (28 5.96 -7961.39)
+ (27 5.09 188489.81)
+ (27 1.72 2132.19)
+ (25 2.56 109771.03)
+ (24 1.92 54868.56)
+ (21 0.09 25443.93)
+ (21 5.98 -55731.43)
+ (20 4.03 60697.74)
+ (18 4.47 2132.79)
+ (17 0.79 109771.63)
+ (14 4.24 -7752.82)
+ (13 2.01 188491.91)
+ (13 2.65 207.81)
+ (13 4.98 29424.63)
+ (12 0.93 -7.99)
+ (10 2.21 46941.14)
+ (10 3.59 -68.29)
+ (10 1.50 21463.25)
+ (10 2.55 157208.40)))
+
+(defun solar-longitude (d)
+ "Longitude of sun on astronomical (Julian) day number D.
+Accurary is about 0.0006 degree (about 365.25*24*60*0.0006/360 = 1 minutes).
+
+The values of calendar-daylight-savings-starts,
+calendar-daylight-savings-starts-time, calendar-daylight-savings-ends,
+calendar-daylight-savings-ends-time, calendar-daylight-time-offset, and
+calendar-time-zone are used to interpret local time."
+ (let* ((a-d (calendar-absolute-from-astro d))
+ ;; get Universal Time
+ (date (calendar-astro-from-absolute
+ (- a-d
+ (if (dst-in-effect a-d)
+ (/ calendar-daylight-time-offset 24.0 60.0) 0)
+ (/ calendar-time-zone 60.0 24.0))))
+ ;; get Ephemeris Time
+ (date (+ date (solar-ephemeris-correction
+ (extract-calendar-year
+ (calendar-gregorian-from-absolute
+ (floor
+ (calendar-absolute-from-astro
+ date)))))))
+ (U (/ (- date 2451545) 3652500))
+ (longitude
+ (+ 4.9353929
+ (* 62833.1961680 U)
+ (* 0.0000001
+ (apply '+
+ (mapcar '(lambda (x)
+ (* (car x)
+ (sin (mod
+ (+ (car (cdr x))
+ (* (car (cdr (cdr x))) U))
+ (* 2 pi)))))
+ solar-data-list)))))
+ (aberration
+ (* 0.0000001 (- (* 17 (cos (+ 3.10 (* 62830.14 U)))) 973)))
+ (A1 (mod (+ 2.18 (* U (+ -3375.70 (* 0.36 U)))) (* 2 pi)))
+ (A2 (mod (+ 3.51 (* U (+ 125666.39 (* 0.10 U)))) (* 2 pi)))
+ (nutation (* -0.0000001 (+ (* 834 (sin A1)) (* 64 (sin A2))))))
+ (mod (radians-to-degrees (+ longitude aberration nutation)) 360.0)))