X-Git-Url: https://code.delx.au/gnu-emacs/blobdiff_plain/fb97eeac20f11e5aa692851472c9f81ede40a71d..a01a7932080e8a6e7bc8472c58cefabcc2c37df3:/lisp/calendar/solar.el

diff --git a/lisp/calendar/solar.el b/lisp/calendar/solar.el
index d0b02b4b11..b7a728461f 100644
--- a/lisp/calendar/solar.el
+++ b/lisp/calendar/solar.el
@@ -1,20 +1,21 @@
 ;;; solar.el --- calendar functions for solar events
 
 ;; Copyright (C) 1992, 1993, 1995, 1997, 2001, 2002, 2003, 2004, 2005,
-;;   2006, 2007, 2008  Free Software Foundation, Inc.
+;;   2006, 2007, 2008, 2009, 2010  Free Software Foundation, Inc.
 
 ;; Author: Edward M. Reingold <reingold@cs.uiuc.edu>
 ;;         Denis B. Roegel <Denis.Roegel@loria.fr>
 ;; Maintainer: Glenn Morris <rgm@gnu.org>
 ;; Keywords: calendar
 ;; Human-Keywords: sunrise, sunset, equinox, solstice, calendar, diary, holidays
+;; Package: calendar
 
 ;; This file is part of GNU Emacs.
 
-;; GNU Emacs is free software; you can redistribute it and/or modify
+;; GNU Emacs is free software: you can redistribute it and/or modify
 ;; it under the terms of the GNU General Public License as published by
-;; the Free Software Foundation; either version 3, or (at your option)
-;; any later version.
+;; the Free Software Foundation, either version 3 of the License, or
+;; (at your option) any later version.
 
 ;; GNU Emacs is distributed in the hope that it will be useful,
 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -22,15 +23,12 @@
 ;; GNU General Public License for more details.
 
 ;; You should have received a copy of the GNU General Public License
-;; along with GNU Emacs; see the file COPYING.  If not, write to the
-;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
-;; Boston, MA 02110-1301, USA.
+;; along with GNU Emacs.  If not, see <http://www.gnu.org/licenses/>.
 
 ;;; Commentary:
 
-;; This collection of functions implements the features of calendar.el,
-;; diary.el, and holiday.el that deal with times of day, sunrise/sunset, and
-;; equinoxes/solstices.
+;; See calendar.el.  This file implements features that deal with
+;; times of day, sunrise/sunset, and equinoxes/solstices.
 
 ;; Based on the ``Almanac for Computers 1984,'' prepared by the Nautical
 ;; Almanac Office, United States Naval Observatory, Washington, 1984, on
@@ -48,21 +46,12 @@
 ;;    2. Equinox/solstice times will be accurate to the minute for years
 ;;       1951--2050.  For other years the times will be within +/- 1 minute.
 
-;; Technical details of all the calendrical calculations can be found in
-;; ``Calendrical Calculations: The Millennium Edition'' by Edward M. Reingold
-;; and Nachum Dershowitz, Cambridge University Press (2001).
-
 ;;; Code:
 
-(defvar displayed-month)
-(defvar displayed-year)
-
-(if (fboundp 'atan)
-    (require 'lisp-float-type)
-  (error "Solar calculations impossible since floating point is unavailable"))
-
+(require 'calendar)
 (require 'cal-dst)
-(require 'cal-julian)
+;; calendar-astro-to-absolute and v versa are cal-autoloads.
+;;;(require 'cal-julian)
 
 
 (defcustom calendar-time-display-form
@@ -128,14 +117,14 @@ This variable should be set in `site-start'.el."
                   (/ (aref calendar-latitude 1) 60.0)))
              (if (numberp calendar-latitude)
                  (if (> calendar-latitude 0) "N" "S")
-               (if (equal (aref calendar-latitude 2) 'north) "N" "S"))
+               (if (eq (aref calendar-latitude 2) 'north) "N" "S"))
              (if (numberp calendar-longitude)
                  (abs calendar-longitude)
                (+ (aref calendar-longitude 0)
                   (/ (aref calendar-longitude 1) 60.0)))
              (if (numberp calendar-longitude)
                  (if (> calendar-longitude 0) "E" "W")
-               (if (equal (aref calendar-longitude 2) 'east) "E" "W"))))
+               (if (eq (aref calendar-longitude 2) 'east) "E" "W"))))
   "Expression evaluating to the name of the calendar location.
 For example, \"New York City\".  The default value is just the
 variable `calendar-latitude' paired with the variable `calendar-longitude'.
@@ -159,26 +148,29 @@ delta.  At present, delta = 0.01 degrees, so the value of the variable
   :type 'number
   :group 'calendar)
 
-(defcustom diary-sabbath-candles-minutes 18
-  "Number of minutes before sunset for sabbath candle lighting."
-  :group 'diary
-  :type 'integer
-  :version "21.1")
-
-
-;;; End of user options.
-
-
-(defvar solar-n-hemi-seasons
+(defcustom solar-n-hemi-seasons
   '("Vernal Equinox" "Summer Solstice" "Autumnal Equinox" "Winter Solstice")
-  "List of season changes for the northern hemisphere.")
+  "List of season changes for the northern hemisphere."
+  :type '(list
+          (string :tag "Vernal Equinox")
+          (string :tag "Summer Solstice")
+          (string :tag "Autumnal Equinox")
+          (string :tag "Winter Solstice"))
+  :group 'calendar)
 
-(defvar solar-s-hemi-seasons
+(defcustom solar-s-hemi-seasons
   '("Autumnal Equinox" "Winter Solstice" "Vernal Equinox" "Summer Solstice")
-  "List of season changes for the southern hemisphere.")
+  "List of season changes for the southern hemisphere."
+  :type '(list
+          (string :tag "Autumnal Equinox")
+          (string :tag "Winter Solstice")
+          (string :tag "Vernal Equinox")
+          (string :tag "Summer Solstice"))
+  :group 'calendar)
+
+;;; End of user options.
 
-(defvar solar-sidereal-time-greenwich-midnight
-  nil
+(defvar solar-sidereal-time-greenwich-midnight nil
   "Sidereal time at Greenwich at midnight (universal time).")
 
 (defvar solar-northern-spring-or-summer-season nil
@@ -192,7 +184,7 @@ Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.")
       calendar-latitude
     (let ((lat (+ (aref calendar-latitude 0)
                   (/ (aref calendar-latitude 1) 60.0))))
-      (if (equal (aref calendar-latitude 2) 'north)
+      (if (eq (aref calendar-latitude 2) 'north)
           lat
         (- lat)))))
 
@@ -202,10 +194,17 @@ Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.")
       calendar-longitude
     (let ((long (+ (aref calendar-longitude 0)
                    (/ (aref calendar-longitude 1) 60.0))))
-      (if (equal (aref calendar-longitude 2) 'east)
+      (if (eq (aref calendar-longitude 2) 'east)
           long
         (- long)))))
 
+(defun solar-get-number (prompt)
+  "Return a number from the minibuffer, prompting with PROMPT.
+Returns nil if nothing was entered."
+  (let ((x (read-string prompt "")))
+    (unless (string-equal x "")
+      (string-to-number x))))
+
 (defun solar-setup ()
   "Prompt for `calendar-longitude', `calendar-latitude', `calendar-time-zone'."
   (beep)
@@ -223,13 +222,6 @@ Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.")
              "Enter difference from Coordinated Universal Time (in minutes): ")
             )))
 
-(defun solar-get-number (prompt)
-  "Return a number from the minibuffer, prompting with PROMPT.
-Returns nil if nothing was entered."
-  (let ((x (read-string prompt "")))
-    (if (not (string-equal x ""))
-        (string-to-number x))))
-
 (defun solar-sin-degrees (x)
   "Return sin of X degrees."
   (sin (degrees-to-radians (mod x 360.0))))
@@ -255,10 +247,10 @@ Returns nil if nothing was entered."
 (defun solar-arctan (x quad)
   "Arctangent of X in quadrant QUAD."
   (let ((deg (radians-to-degrees (atan x))))
-    (cond ((equal quad 2) (+ deg 180))
-          ((equal quad 3) (+ deg 180))
-          ((equal quad 4) (+ deg 360))
-          (t              deg))))
+    (cond ((= quad 2) (+ deg 180))
+          ((= quad 3) (+ deg 180))
+          ((= quad 4) (+ deg 360))
+          (t          deg))))
 
 (defun solar-atn2 (x y)
   "Arctangent of point X, Y."
@@ -299,36 +291,182 @@ Both arguments are in degrees."
    (* (solar-sin-degrees obliquity)
       (solar-sin-degrees longitude))))
 
-(defun solar-sunrise-and-sunset (time latitude longitude height)
-  "Sunrise, sunset and length of day.
-Parameters are the midday TIME and the LATITUDE, LONGITUDE of the location.
+(defun solar-ecliptic-coordinates (time sunrise-flag)
+  "Return solar longitude, ecliptic inclination, equation of time, nutation.
+Values are for TIME in Julian centuries of Ephemeris Time since
+January 1st, 2000, at 12 ET.  Longitude and inclination are in
+degrees, equation of time in hours, and nutation in seconds of longitude.
+If SUNRISE-FLAG is non-nil, only calculate longitude and inclination."
+  (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
+         ;; Longitude of moon's ascending node on the ecliptic.
+         (omega (+ 125.04
+                   (* -1934.136 time)))
+         ;; nut = nutation in longitude, measured in seconds of angle.
+         (nut (unless sunrise-flag
+                (+ (* -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))))))
+         (ecc (unless sunrise-flag     ; eccentricity of earth's orbit
+                (+ 0.016708617
+                   (* -0.000042037 time)
+                   (* -0.0000001236 time time))))
+         (app (+ L                      ; apparent longitude of sun
+                 -0.00569
+                 (* -0.00478
+                    (solar-sin-degrees omega))))
+         (y (unless sunrise-flag
+              (* (solar-tangent-degrees (/ i 2))
+                 (solar-tangent-degrees (/ i 2)))))
+         ;; Equation of time, in hours.
+         (time-eq (unless sunrise-flag
+                    (/ (* 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))))
+    (list app i time-eq nut)))
+
+(defun solar-ephemeris-correction (year)
+  "Ephemeris time minus Universal Time during Gregorian YEAR.
+Result is in days.  For the years 1800-1987, the maximum error is
+1.9 seconds.  For the other years, the maximum error is about 30 seconds."
+  (cond ((and (<= 1988 year) (< year 2020))
+         (/ (+ year -2000 67.0) 60.0 60.0 24.0))
+        ((and (<= 1900 year) (< year 1988))
+         (let* ((theta (/ (- (calendar-astro-from-absolute
+                              (calendar-absolute-from-gregorian
+                               (list 7 1 year)))
+                             (calendar-astro-from-absolute
+                              (calendar-absolute-from-gregorian
+                               '(1 1 1900))))
+                          36525.0))
+                (theta2 (* theta theta))
+                (theta3 (* theta2 theta))
+                (theta4 (* theta2 theta2))
+                (theta5 (* theta3 theta2)))
+           (+ -0.00002
+              (* 0.000297 theta)
+              (* 0.025184 theta2)
+              (* -0.181133 theta3)
+              (* 0.553040 theta4)
+              (* -0.861938 theta5)
+              (* 0.677066 theta3 theta3)
+              (* -0.212591 theta4 theta3))))
+        ((and (<= 1800 year) (< year 1900))
+         (let* ((theta (/ (- (calendar-astro-from-absolute
+                              (calendar-absolute-from-gregorian
+                               (list 7 1 year)))
+                             (calendar-astro-from-absolute
+                              (calendar-absolute-from-gregorian
+                               '(1 1 1900))))
+                          36525.0))
+                (theta2 (* theta theta))
+                (theta3 (* theta2 theta))
+                (theta4 (* theta2 theta2))
+                (theta5 (* theta3 theta2)))
+           (+ -0.000009
+              (* 0.003844 theta)
+              (* 0.083563 theta2)
+              (* 0.865736 theta3)
+              (* 4.867575 theta4)
+              (* 15.845535 theta5)
+              (* 31.332267 theta3 theta3)
+              (* 38.291999 theta4 theta3)
+              (* 28.316289 theta4 theta4)
+              (* 11.636204 theta4 theta5)
+              (* 2.043794 theta5 theta5))))
+        ((and (<= 1620 year) (< year 1800))
+         (let ((x (/ (- year 1600) 10.0)))
+           (/ (+ (* 2.19167 x x) (* -40.675 x) 196.58333) 60.0 60.0 24.0)))
+        (t (let* ((tmp (- (calendar-astro-from-absolute
+                           (calendar-absolute-from-gregorian
+                            (list 1 1 year)))
+                          2382148))
+                  (second (- (/ (* tmp tmp) 41048480.0) 15)))
+             (/ second 60.0 60.0 24.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.
 
-HEIGHT is the angle the center of the sun has over the horizon for the contact
-we are trying to find.  For sunrise and sunset, it is usually -0.61 degrees,
-accounting for the edge of the sun being on the horizon.
+Result is in Julian centuries of ephemeris time."
+  (let* ((t0 (car time))
+         (ut (cadr time))
+         (t1 (+ t0 (/ (/ ut 24.0) 36525)))
+         (y (+ 2000 (* 100 t1)))
+         (dt (* 86400 (solar-ephemeris-correction (floor y)))))
+    (+ t1 (/ (/ dt 86400) 36525))))
 
-Coordinates are included because this function is called with latitude=1
-degrees to find out if polar regions have 24 hours of sun or only night."
-  (let* ((rise-time (solar-moment -1 latitude longitude time height))
-         (set-time (solar-moment 1 latitude longitude time height))
-         (day-length))
-    (if (not (and rise-time set-time))
-        (if (or (and (> latitude 0)
-                     solar-northern-spring-or-summer-season)
-                (and (< latitude 0)
-                     (not solar-northern-spring-or-summer-season)))
-            (setq day-length 24)
-          (setq day-length 0))
-      (setq day-length (- set-time rise-time)))
-    (list (if rise-time (+ rise-time (/ calendar-time-zone 60.0)) nil)
-          (if set-time (+ set-time (/ calendar-time-zone 60.0)) nil)
-          day-length)))
+(defun solar-equatorial-coordinates (time sunrise-flag)
+  "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.  SUNRISE-FLAG is passed
+to `solar-ecliptic-coordinates'."
+  (let ((ec (solar-ecliptic-coordinates (solar-ephemeris-time time)
+                                        sunrise-flag)))
+    (list (solar-right-ascension (car ec) (cadr ec))
+          (solar-declination (car ec) (cadr ec)))))
+
+(defun solar-horizontal-coordinates (time latitude longitude sunrise-flag)
+  "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.  SUNRISE-FLAG
+is passed to `solar-ecliptic-coordinates'.  Azimuth and
+height (between -180 and 180) are both in degrees."
+  (let* ((ut (cadr time))
+         (ec (solar-equatorial-coordinates time sunrise-flag))
+         (st (+ solar-sidereal-time-greenwich-midnight
+                (* ut 1.00273790935)))
+         ;; Hour angle (in degrees).
+         (ah (- (* st 15) (* 15 (car ec)) (* -1 (calendar-longitude))))
+         (de (cadr 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-moment (direction latitude longitude time height)
   "Sunrise/sunset at location.
@@ -377,6 +515,37 @@ Uses binary search."
       (setq possible nil))              ; the sun never sets
     (if possible utmoment)))
 
+(defun solar-sunrise-and-sunset (time latitude longitude height)
+  "Sunrise, sunset and length of day.
+Parameters are the midday TIME and the LATITUDE, LONGITUDE of the location.
+
+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.
+
+HEIGHT is the angle the center of the sun has over the horizon for the contact
+we are trying to find.  For sunrise and sunset, it is usually -0.61 degrees,
+accounting for the edge of the sun being on the horizon.
+
+Coordinates are included because this function is called with latitude=1
+degrees to find out if polar regions have 24 hours of sun or only night."
+  (let ((rise-time (solar-moment -1 latitude longitude time height))
+        (set-time (solar-moment 1 latitude longitude time height))
+        day-length)
+    (if (not (and rise-time set-time))
+        (if (or (and (> latitude 0)
+                     solar-northern-spring-or-summer-season)
+                (and (< latitude 0)
+                     (not solar-northern-spring-or-summer-season)))
+            (setq day-length 24)
+          (setq day-length 0))
+      (setq day-length (- set-time rise-time)))
+    (list (if rise-time (+ rise-time (/ calendar-time-zone 60.0)) nil)
+          (if set-time (+ set-time (/ calendar-time-zone 60.0)) nil)
+          day-length)))
+
 (defun solar-time-string (time time-zone)
   "Printable form for decimal fraction TIME in TIME-ZONE.
 Format used is given by `calendar-time-display-form'."
@@ -388,13 +557,27 @@ Format used is given by `calendar-time-display-form'."
          (24-hours (format "%02d" 24-hours)))
     (mapconcat 'eval calendar-time-display-form "")))
 
-
 (defun solar-daylight (time)
   "Printable form for TIME expressed in hours."
   (format "%d:%02d"
           (floor time)
           (floor (* 60 (- time (floor time))))))
 
+(defun solar-julian-ut-centuries (date)
+  "Number of Julian centuries since 1 Jan, 2000 at noon UT for Gregorian DATE."
+  (/ (- (calendar-absolute-from-gregorian date)
+        (calendar-absolute-from-gregorian '(1 1.5 2000)))
+     36525.0))
+
+(defun solar-date-to-et (date ut)
+  "Ephemeris Time at Gregorian DATE at Universal Time UT (in hours).
+Expressed in Julian centuries of Ephemeris Time."
+  (solar-ephemeris-time (list (solar-julian-ut-centuries date) ut)))
+
+(defun solar-time-equation (date ut)
+  "Equation of time expressed in hours at Gregorian DATE at Universal time UT."
+  (nth 2 (solar-ecliptic-coordinates (solar-date-to-et date ut) nil)))
+
 (defun solar-exact-local-noon (date)
   "Date and Universal Time of local noon at *local date* DATE.
 The date may be different from the one asked for, but it will be the right
@@ -415,12 +598,28 @@ local date.  The second component of date should be an integer."
               (calendar-absolute-from-gregorian nd)))
     (list nd ut)))
 
-(defun solar-sunrise-sunset (date)
-  "List of *local* times of sunrise, sunset, and daylight on Gregorian DATE.
-Corresponding value is nil if there is no sunrise/sunset."
-  ;; First, get the exact moment of local noon.
-  (let* ((exact-local-noon (solar-exact-local-noon date))
-         ;; Get the time from the 2000 epoch.
+(defun solar-sidereal-time (t0)
+  "Sidereal time (in hours) in Greenwich at T0 Julian centuries.
+T0 must correspond to 0 hours UT."
+  (let* ((mean-sid-time (+ 6.6973746
+                           (* 2400.051337 t0)
+                           (* 0.0000258622 t0 t0)
+                           (* -0.0000000017222 t0 t0 t0)))
+         (et (solar-ephemeris-time (list t0 0.0)))
+         (nut-i (solar-ecliptic-coordinates et nil))
+         (nut (nth 3 nut-i))            ; nutation
+         (i (cadr nut-i)))              ; inclination
+    (mod (+ (mod (+ mean-sid-time
+                    (/ (/ (* nut (solar-cosine-degrees i)) 15) 3600)) 24.0)
+            24.0)
+         24.0)))
+
+(defun solar-sunrise-sunset (date)
+  "List of *local* times of sunrise, sunset, and daylight on Gregorian DATE.
+Corresponding value is nil if there is no sunrise/sunset."
+  ;; First, get the exact moment of local noon.
+  (let* ((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.
@@ -443,185 +642,32 @@ Corresponding value is nil if there is no sunrise/sunset."
              (list t0 (cadr exact-local-noon))
              (calendar-latitude)
              (calendar-longitude) -0.61)))
-         (rise (car rise-set))
-         (adj-rise (if rise (dst-adjust-time date rise)))
-         (set (cadr rise-set))          ; FIXME ?
-         (adj-set (if set (dst-adjust-time date set)))
+         (rise-time (car rise-set))
+         (adj-rise (if rise-time (dst-adjust-time date rise-time)))
+         (set-time (cadr rise-set))
+         (adj-set (if set-time (dst-adjust-time date set-time)))
          (length (nth 2 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))
+     (and rise-time (calendar-date-equal date (car adj-rise)) (cdr adj-rise))
+     (and set-time (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."
+(defun solar-sunrise-sunset-string (date &optional nolocation)
+  "String of *local* times of sunrise, sunset, and daylight on Gregorian DATE.
+Optional NOLOCATION non-nil means do not print the location."
   (let ((l (solar-sunrise-sunset date)))
     (format
-     "%s, %s at %s (%s hours daylight)"
+     "%s, %s%s (%s hours daylight)"
      (if (car l)
          (concat "Sunrise " (apply 'solar-time-string (car l)))
        "No sunrise")
      (if (cadr l)
          (concat "sunset " (apply 'solar-time-string (cadr l)))
        "no sunset")
-     (eval calendar-location-name)
+     (if nolocation ""
+       (format " at %s" (eval calendar-location-name)))
      (nth 2 l))))
 
-(defun solar-julian-ut-centuries (date)
-  "Number of Julian centuries since 1 Jan, 2000 at noon UT 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 (cadr 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 time after day D when solar longitude is a multiple of L degrees.
-D is a Julian day number.  L must be an integer divisor of 360.
-The result is for `calendar-location-name', and is in local time
-\(including any daylight saving rules) expressed in 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)
-            long (solar-longitude d))
-      (if (or (and (not (zerop next)) (< long next))
-              (and (zerop next) (< l long)))
-          (setq start d
-                start-long long)
-        (setq end d
-              end-long long)))
-    (/ (+ start end) 2.0)))
-
-(defun solar-horizontal-coordinates (time latitude longitude sunrise-flag)
-  "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.  SUNRISE-FLAG
-is passed to `solar-ecliptic-coordinates'.  Azimuth and
-height (between -180 and 180) are both in degrees."
-  (let* ((ut (cadr time))
-         (ec (solar-equatorial-coordinates time sunrise-flag))
-         (st (+ solar-sidereal-time-greenwich-midnight
-                (* ut 1.00273790935)))
-         ;; Hour angle (in degrees).
-         (ah (- (* st 15) (* 15 (car ec)) (* -1 (calendar-longitude))))
-         (de (cadr 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 sunrise-flag)
-  "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.  SUNRISE-FLAG is passed
-to `solar-ecliptic-coordinates'."
-  (let* ((tm (solar-ephemeris-time time))
-         (ec (solar-ecliptic-coordinates tm sunrise-flag)))
-    (list (solar-right-ascension (car ec) (car (cdr ec)))
-          (solar-declination (car ec) (car (cdr ec))))))
-
-(defun solar-ecliptic-coordinates (time sunrise-flag)
-  "Return solar longitude, ecliptic inclination, equation of time, nutation.
-Values are for TIME in Julian centuries of Ephemeris Time since
-January 1st, 2000, at 12 ET.  Longitude and inclination are in
-degrees, equation of time in hours, and nutation in seconds of longitude.
-If SUNRISE-FLAG is non-nil, only calculate longitude and inclination."
-  (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
-         ;; Longitude of moon's ascending node on the ecliptic.
-         (omega (+ 125.04
-                   (* -1934.136 time)))
-         ;; nut = nutation in longitude, measured in seconds of angle.
-         (nut (unless sunrise-flag
-                (+ (* -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))))))
-         (ecc (unless sunrise-flag     ; eccentricity of earth's orbit
-                (+ 0.016708617
-                   (* -0.000042037 time)
-                   (* -0.0000001236 time time))))
-         (app (+ L                      ; apparent longitude of sun
-                 -0.00569
-                 (* -0.00478
-                    (solar-sin-degrees omega))))
-         (y (unless sunrise-flag
-              (* (solar-tangent-degrees (/ i 2))
-                 (solar-tangent-degrees (/ i 2)))))
-         ;; Equation of time, in hours.
-         (time-eq (unless sunrise-flag
-                    (/ (* 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))))
-    (list app i time-eq nut)))
-
 (defconst solar-data-list
   '((403406 4.721964 1.621043)
     (195207 5.937458 62830.348067)
@@ -682,7 +728,7 @@ 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))
+  (let* ((a-d (calendar-astro-to-absolute d))
          ;; Get Universal Time.
          (date (calendar-astro-from-absolute
                 (- a-d
@@ -691,10 +737,10 @@ The values of `calendar-daylight-savings-starts',
                    (/ calendar-time-zone 60.0 24.0))))
          ;; Get Ephemeris Time.
          (date (+ date (solar-ephemeris-correction
-                        (extract-calendar-year
+                        (calendar-extract-year
                          (calendar-gregorian-from-absolute
                           (floor
-                           (calendar-absolute-from-astro
+                           (calendar-astro-to-absolute
                             date)))))))
          (U (/ (- date 2451545) 3652500))
          (longitude
@@ -705,101 +751,49 @@ The values of `calendar-daylight-savings-starts',
                        (mapcar (lambda (x)
                                  (* (car x)
                                     (sin (mod
-                                          (+ (car (cdr x))
-                                             (* (car (cdr (cdr x))) U))
-                                          (* 2 pi)))))
+                                          (+ (cadr x)
+                                             (* (nth 2 x) U))
+                                          (* 2 float-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)))
+         (A1 (mod (+ 2.18 (* U (+ -3375.70 (* 0.36 U)))) (* 2 float-pi)))
+         (A2 (mod (+ 3.51 (* U (+ 125666.39 (* 0.10 U)))) (* 2 float-pi)))
          (nutation (* -0.0000001 (+ (* 834 (sin A1)) (* 64 (sin A2))))))
     (mod (radians-to-degrees (+ longitude aberration nutation)) 360.0)))
 
-(defun solar-ephemeris-correction (year)
-  "Ephemeris time minus Universal Time during Gregorian YEAR.
-Result is in days.  For the years 1800-1987, the maximum error is
-1.9 seconds.  For the other years, the maximum error is about 30 seconds."
-  (cond ((and (<= 1988 year) (< year 2020))
-         (/ (+ year -2000 67.0) 60.0 60.0 24.0))
-        ((and (<= 1900 year) (< year 1988))
-         (let* ((theta (/ (- (calendar-astro-from-absolute
-                              (calendar-absolute-from-gregorian
-                               (list 7 1 year)))
-                             (calendar-astro-from-absolute
-                              (calendar-absolute-from-gregorian
-                               '(1 1 1900))))
-                          36525.0))
-                (theta2 (* theta theta))
-                (theta3 (* theta2 theta))
-                (theta4 (* theta2 theta2))
-                (theta5 (* theta3 theta2)))
-           (+ -0.00002
-              (* 0.000297 theta)
-              (* 0.025184 theta2)
-              (* -0.181133 theta3)
-              (* 0.553040 theta4)
-              (* -0.861938 theta5)
-              (* 0.677066 theta3 theta3)
-              (* -0.212591 theta4 theta3))))
-        ((and (<= 1800 year) (< year 1900))
-         (let* ((theta (/ (- (calendar-astro-from-absolute
-                              (calendar-absolute-from-gregorian
-                               (list 7 1 year)))
-                             (calendar-astro-from-absolute
-                              (calendar-absolute-from-gregorian
-                               '(1 1 1900))))
-                          36525.0))
-                (theta2 (* theta theta))
-                (theta3 (* theta2 theta))
-                (theta4 (* theta2 theta2))
-                (theta5 (* theta3 theta2)))
-           (+ -0.000009
-              (* 0.003844 theta)
-              (* 0.083563 theta2)
-              (* 0.865736 theta3)
-              (* 4.867575 theta4)
-              (* 15.845535 theta5)
-              (* 31.332267 theta3 theta3)
-              (* 38.291999 theta4 theta3)
-              (* 28.316289 theta4 theta4)
-              (* 11.636204 theta4 theta5)
-              (* 2.043794 theta5 theta5))))
-        ((and (<= 1620 year) (< year 1800))
-         (let ((x (/ (- year 1600) 10.0)))
-           (/ (+ (* 2.19167 x x) (* -40.675 x) 196.58333) 60.0 60.0 24.0)))
-        (t (let* ((tmp (- (calendar-astro-from-absolute
-                           (calendar-absolute-from-gregorian
-                            (list 1 1 year)))
-                          2382148))
-                  (second (- (/ (* tmp tmp) 41048480.0) 15)))
-             (/ second 60.0 60.0 24.0)))))
-
-(defun solar-sidereal-time (t0)
-  "Sidereal time (in hours) in Greenwich at T0 Julian centuries.
-T0 must correspond to 0 hours UT."
-  (let* ((mean-sid-time (+ 6.6973746
-                           (* 2400.051337 t0)
-                           (* 0.0000258622 t0 t0)
-                           (* -0.0000000017222 t0 t0 t0)))
-         (et (solar-ephemeris-time (list t0 0.0)))
-         (nut-i (solar-ecliptic-coordinates et nil))
-         (nut (nth 3 nut-i))            ; nutation
-         (i (cadr nut-i)))              ; inclination
-    (mod (+ (mod (+ mean-sid-time
-                    (/ (/ (* nut (solar-cosine-degrees i)) 15) 3600)) 24.0)
-            24.0)
-         24.0)))
-
-(defun solar-time-equation (date ut)
-  "Equation of time expressed in hours at Gregorian DATE at Universal time UT."
-  (nth 2 (solar-ecliptic-coordinates (solar-date-to-et date ut) nil)))
-
-(defun solar-date-to-et (date ut)
-  "Ephemeris Time at Gregorian DATE at Universal Time UT (in hours).
-Expressed in Julian centuries of Ephemeris Time."
-  (solar-ephemeris-time (list (solar-julian-ut-centuries date) ut)))
+(defun solar-date-next-longitude (d l)
+  "First time after day D when solar longitude is a multiple of L degrees.
+D is a Julian day number.  L must be an integer divisor of 360.
+The result is for `calendar-location-name', and is in local time
+\(including any daylight saving rules) expressed in 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)
+            long (solar-longitude d))
+      (if (or (and (not (zerop next)) (< long next))
+              (and (zerop next) (< l long)))
+          (setq start d
+                start-long long)
+        (setq end d
+              end-long long)))
+    (/ (+ start end) 2.0)))
 
+;; FIXME but there already is solar-sunrise-sunset.
 ;;;###autoload
 (defun sunrise-sunset (&optional arg)
   "Local time of sunrise and sunset for today.  Accurate to a few seconds.
@@ -835,18 +829,18 @@ This function is suitable for execution in a .emacs file."
                            (/ (aref calendar-latitude 1) 60.0)))
                       (if (numberp calendar-latitude)
                           (if (> calendar-latitude 0) "N" "S")
-                        (if (equal (aref calendar-latitude 2) 'north) "N" "S"))
+                        (if (eq (aref calendar-latitude 2) 'north) "N" "S"))
                       (if (numberp calendar-longitude)
                           (abs calendar-longitude)
                         (+ (aref calendar-longitude 0)
                            (/ (aref calendar-longitude 1) 60.0)))
                       (if (numberp calendar-longitude)
                           (if (> calendar-longitude 0) "E" "W")
-                        (if (equal (aref calendar-longitude 2) 'east)
+                        (if (eq (aref calendar-longitude 2) 'east)
                             "E" "W"))))))
          (calendar-standard-time-zone-name
           (if (< arg 16) calendar-standard-time-zone-name
-            (cond ((= calendar-time-zone 0) "UTC")
+            (cond ((zerop calendar-time-zone) "UTC")
                   ((< calendar-time-zone 0)
                    (format "UTC%dmin" calendar-time-zone))
                   (t  (format "UTC+%dmin" calendar-time-zone)))))
@@ -873,20 +867,41 @@ This function is suitable for execution in a .emacs file."
 contents of temp window."))))))
 
 ;;;###cal-autoload
-(defun calendar-sunrise-sunset ()
+(defun calendar-sunrise-sunset (&optional event)
   "Local time of sunrise and sunset for date under cursor.
 Accurate to a few seconds."
-  (interactive)
+  (interactive (list last-nonmenu-event))
   (or (and calendar-latitude calendar-longitude calendar-time-zone)
       (solar-setup))
-  (let ((date (calendar-cursor-to-date t)))
+  (let ((date (calendar-cursor-to-date t event)))
     (message "%s: %s"
              (calendar-date-string date t t)
              (solar-sunrise-sunset-string date))))
 
+;;;###cal-autoload
+(defun calendar-sunrise-sunset-month (&optional event)
+  "Local time of sunrise and sunset for month under cursor or at EVENT."
+  (interactive (list last-nonmenu-event))
+  (or (and calendar-latitude calendar-longitude calendar-time-zone)
+      (solar-setup))
+  (let* ((date (calendar-cursor-to-date t event))
+         (month (car date))
+         (year (nth 2 date))
+         (last (calendar-last-day-of-month month year))
+         (title (format "Sunrise/sunset times for %s %d at %s"
+                        (calendar-month-name month) year
+                        (eval calendar-location-name))))
+    (calendar-in-read-only-buffer solar-sunrises-buffer
+      (calendar-set-mode-line title)
+      (insert title ":\n\n")
+      (dotimes (i last)
+        (setq date (list month (1+ i) year))
+        (insert (format "%s %2d: " (calendar-month-name month t) (1+ i))
+                (solar-sunrise-sunset-string date t) "\n")))))
+
 (defvar date)
 
-;; To be called from list-sexp-diary-entries, where DATE is bound.
+;; To be called from diary-list-sexp-entries, where DATE is bound.
 ;;;###diary-autoload
 (defun diary-sunrise-sunset ()
   "Local time of sunrise and sunset as a diary entry.
@@ -895,27 +910,6 @@ Accurate to a few seconds."
       (solar-setup))
   (solar-sunrise-sunset-string date))
 
-;; To be called from list-sexp-diary-entries, where DATE is bound.
-;;;###diary-autoload
-(defun diary-sabbath-candles (&optional mark)
-  "Local time of candle lighting diary entry--applies if date is a Friday.
-No diary entry if there is no sunset on that date.
-
-An optional parameter MARK specifies a face or single-character string to
-use when highlighting the day in the calendar."
-  (or (and calendar-latitude calendar-longitude calendar-time-zone)
-      (solar-setup))
-  (if (= (% (calendar-absolute-from-gregorian date) 7) 5) ; Friday
-      (let* ((sunset (cadr (solar-sunrise-sunset date)))
-             (light (if sunset
-                        (cons (- (car sunset)
-                                 (/ diary-sabbath-candles-minutes 60.0))
-                              (cdr sunset)))))
-        (if sunset
-            (cons mark
-                  (format "%s Sabbath candle lighting"
-                          (apply 'solar-time-string light)))))))
-
 ;; From Meeus, 1991, page 167.
 (defconst solar-seasons-data
   '((485 324.96 1934.136)
@@ -972,97 +966,101 @@ Accurate to within a minute between 1951 and 2050."
 (defun solar-mean-equinoxes/solstices (k year)
   "Julian day of mean equinox/solstice K for YEAR.
 K=0, spring equinox; K=1, summer solstice; K=2, fall equinox; K=3, winter
-solstice.  These formulas are only to be used between 1000 BC and 3000 AD."
+solstice.  These formulae are only to be used between 1000 BC and 3000 AD."
   (let ((y (/ year 1000.0))
         (z (/ (- year 2000) 1000.0)))
     (if (< year 1000)                ; actually between -1000 and 1000
-        (cond ((equal k 0) (+ 1721139.29189
-                              (*  365242.13740 y)
-                              (* 0.06134 y y)
-                              (* 0.00111 y y y)
-                              (* -0.00071 y y y y)))
-              ((equal k 1) (+ 1721233.25401
-                              (* 365241.72562 y)
-                              (* -0.05323 y y)
-                              (* 0.00907 y y y)
-                              (* 0.00025 y y y y)))
-              ((equal k 2) (+ 1721325.70455
-                              (* 365242.49558 y)
-                              (* -0.11677 y y)
-                              (* -0.00297 y y y)
-                              (* 0.00074 y y y y)))
-              ((equal k 3) (+ 1721414.39987
-                              (* 365242.88257 y)
-                              (* -0.00769 y y)
-                              (* -0.00933 y y y)
-                              (* -0.00006 y y y y))))
+        (cond ((= k 0) (+ 1721139.29189
+                          (*  365242.13740 y)
+                          (* 0.06134 y y)
+                          (* 0.00111 y y y)
+                          (* -0.00071 y y y y)))
+              ((= k 1) (+ 1721233.25401
+                          (* 365241.72562 y)
+                          (* -0.05323 y y)
+                          (* 0.00907 y y y)
+                          (* 0.00025 y y y y)))
+              ((= k 2) (+ 1721325.70455
+                          (* 365242.49558 y)
+                          (* -0.11677 y y)
+                          (* -0.00297 y y y)
+                          (* 0.00074 y y y y)))
+              ((= k 3) (+ 1721414.39987
+                          (* 365242.88257 y)
+                          (* -0.00769 y y)
+                          (* -0.00933 y y y)
+                          (* -0.00006 y y y y))))
                                         ; actually between 1000 and 3000
-      (cond ((equal k 0) (+ 2451623.80984
-                            (* 365242.37404  z)
-                            (* 0.05169 z z)
-                            (* -0.00411 z z z)
-                            (* -0.00057 z z z z)))
-            ((equal k 1) (+ 2451716.56767
-                            (* 365241.62603 z)
-                            (* 0.00325 z z)
-                            (* 0.00888 z z z)
-                            (* -0.00030 z z z z)))
-            ((equal k 2) (+ 2451810.21715
-                            (* 365242.01767 z)
-                            (* -0.11575 z z)
-                            (* 0.00337 z z z)
-                            (* 0.00078 z z z z)))
-            ((equal k 3) (+ 2451900.05952
-                            (* 365242.74049 z)
-                            (* -0.06223 z z)
-                            (* -0.00823 z z z)
-                            (* 0.00032 z z z z)))))))
+      (cond ((= k 0) (+ 2451623.80984
+                        (* 365242.37404  z)
+                        (* 0.05169 z z)
+                        (* -0.00411 z z z)
+                        (* -0.00057 z z z z)))
+            ((= k 1) (+ 2451716.56767
+                        (* 365241.62603 z)
+                        (* 0.00325 z z)
+                        (* 0.00888 z z z)
+                        (* -0.00030 z z z z)))
+            ((= k 2) (+ 2451810.21715
+                        (* 365242.01767 z)
+                        (* -0.11575 z z)
+                        (* 0.00337 z z z)
+                        (* 0.00078 z z z z)))
+            ((= k 3) (+ 2451900.05952
+                        (* 365242.74049 z)
+                        (* -0.06223 z z)
+                        (* -0.00823 z z z)
+                        (* 0.00032 z z z z)))))))
+
+(defvar displayed-month)                ; from calendar-generate
+(defvar displayed-year)
 
 ;;;###holiday-autoload
 (defun solar-equinoxes-solstices ()
   "Local date and time of equinoxes and solstices, if visible in the calendar.
 Requires floating point."
-  (let ((m displayed-month)
-        (y displayed-year))
-    (increment-calendar-month m y (cond ((= 1 (% m 3)) -1)
-                                        ((= 2 (% m 3))  1)
-                                        (t              0)))
-    (let* ((calendar-standard-time-zone-name
-            (if calendar-time-zone calendar-standard-time-zone-name "UTC"))
-           (calendar-daylight-savings-starts
-            (if calendar-time-zone calendar-daylight-savings-starts))
-           (calendar-daylight-savings-ends
-            (if calendar-time-zone calendar-daylight-savings-ends))
-           (calendar-time-zone (if calendar-time-zone calendar-time-zone 0))
-           (k (1- (/ m 3)))
-           (d0 (solar-equinoxes/solstices k y))
-           (d1 (list (car d0) (floor (car (cdr d0))) (car (cdr (cdr d0)))))
-           (h0 (* 24 (- (car (cdr d0)) (floor (car (cdr d0))))))
-           (adj (dst-adjust-time d1 h0))
-           (d (list (caar adj)
-                    (+ (car (cdar adj))
-                       (/ (cadr adj) 24.0))
-                    (cadr (cdar adj))))
-           ;; The following is nearly as accurate, but not quite:
-           ;; (d0 (solar-date-next-longitude
-           ;;     (calendar-astro-from-absolute
-           ;;      (calendar-absolute-from-gregorian
-           ;;       (list (+ 3 (* k 3)) 15 y)))
-           ;;     90))
-           ;; (abs-day (calendar-absolute-from-astro d)))
-           (abs-day (calendar-absolute-from-gregorian d)))
-      (list
-       (list (calendar-gregorian-from-absolute (floor abs-day))
-             (format "%s %s"
-                     (nth k (if (and calendar-latitude
-                                     (< (calendar-latitude) 0))
-                                solar-s-hemi-seasons
-                              solar-n-hemi-seasons))
-                     (solar-time-string
-                      (* 24 (- abs-day (floor abs-day)))
-                      (if (dst-in-effect abs-day)
-                          calendar-daylight-time-zone-name
-                        calendar-standard-time-zone-name))))))))
+  (let* ((m displayed-month)
+         (y displayed-year)
+         (calendar-standard-time-zone-name
+          (if calendar-time-zone calendar-standard-time-zone-name "UTC"))
+         (calendar-daylight-savings-starts
+          (if calendar-time-zone calendar-daylight-savings-starts))
+         (calendar-daylight-savings-ends
+          (if calendar-time-zone calendar-daylight-savings-ends))
+         (calendar-time-zone (if calendar-time-zone calendar-time-zone 0))
+         (k (progn
+              (calendar-increment-month m y (cond ((= 1 (% m 3)) -1)
+                                                  ((= 2 (% m 3))  1)
+                                                  (t              0)))
+              (1- (/ m 3))))
+         (d0 (solar-equinoxes/solstices k y))
+         (d1 (list (car d0) (floor (cadr d0)) (nth 2 d0)))
+         (h0 (* 24 (- (cadr d0) (floor (cadr d0)))))
+         (adj (dst-adjust-time d1 h0))
+         (d (list (caar adj)
+                  (+ (car (cdar adj))
+                     (/ (cadr adj) 24.0))
+                  (cadr (cdar adj))))
+         ;; The following is nearly as accurate, but not quite:
+         ;; (d0 (solar-date-next-longitude
+         ;;     (calendar-astro-from-absolute
+         ;;      (calendar-absolute-from-gregorian
+         ;;       (list (+ 3 (* k 3)) 15 y)))
+         ;;     90))
+         ;; (abs-day (calendar-astro-to-absolute d)))
+         (abs-day (calendar-absolute-from-gregorian d)))
+    (list
+     (list (calendar-gregorian-from-absolute (floor abs-day))
+           (format "%s %s"
+                   (nth k (if (and calendar-latitude
+                                   (< (calendar-latitude) 0))
+                              solar-s-hemi-seasons
+                            solar-n-hemi-seasons))
+                   (solar-time-string
+                    (* 24 (- abs-day (floor abs-day)))
+                    (if (dst-in-effect abs-day)
+                        calendar-daylight-time-zone-name
+                      calendar-standard-time-zone-name)))))))
 
 
 (provide 'solar)