X-Git-Url: https://code.delx.au/gnu-emacs/blobdiff_plain/4586442a5abbb9ccd6e7c4de0730763c0170cc12..a01a7932080e8a6e7bc8472c58cefabcc2c37df3:/lisp/calendar/solar.el diff --git a/lisp/calendar/solar.el b/lisp/calendar/solar.el index 4281c0e130..b7a728461f 100644 --- a/lisp/calendar/solar.el +++ b/lisp/calendar/solar.el @@ -1,21 +1,21 @@ ;;; solar.el --- calendar functions for solar events ;; Copyright (C) 1992, 1993, 1995, 1997, 2001, 2002, 2003, 2004, 2005, -;; 2006, 2007 Free Software Foundation, Inc. +;; 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc. ;; Author: Edward M. Reingold -;; Denis B. Roegel +;; Denis B. Roegel ;; Maintainer: Glenn Morris ;; Keywords: calendar -;; Human-Keywords: sunrise, sunset, equinox, solstice, calendar, diary, -;; holidays +;; 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 2, 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 @@ -23,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 . ;;; 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 @@ -49,28 +46,18 @@ ;; 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 date) -(defvar displayed-month) -(defvar displayed-year) - -(if (fboundp 'atan) - (require 'lisp-float-type) - (error "Solar/lunar 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) + -;;;###autoload (defcustom calendar-time-display-form '(12-hours ":" minutes am-pm - (if time-zone " (") time-zone (if time-zone ")")) - "*The pseudo-pattern that governs the way a time of day is formatted. + (if time-zone " (") time-zone (if time-zone ")")) + "The pseudo-pattern that governs the way a time of day is formatted. A pseudo-pattern is a list of expressions that can involve the keywords `12-hours', `24-hours', and `minutes', all numbers in string form, @@ -85,10 +72,8 @@ would give military-style times like `21:07 (UTC)'." :type 'sexp :group 'calendar) -;;;###autoload (defcustom calendar-latitude nil - "*Latitude of `calendar-location-name' in degrees. - + "Latitude of `calendar-location-name' in degrees. The value can be either a decimal fraction (one place of accuracy is sufficient), + north, - south, such as 40.7 for New York City, or the value can be a vector [degrees minutes north/south] such as [40 50 north] for New @@ -96,19 +81,17 @@ York City. This variable should be set in `site-start'.el." :type '(choice (const nil) - (number :tag "Exact") - (vector :value [0 0 north] - (integer :tag "Degrees") - (integer :tag "Minutes") - (choice :tag "Position" - (const north) - (const south)))) + (number :tag "Exact") + (vector :value [0 0 north] + (integer :tag "Degrees") + (integer :tag "Minutes") + (choice :tag "Position" + (const north) + (const south)))) :group 'calendar) -;;;###autoload (defcustom calendar-longitude nil - "*Longitude of `calendar-location-name' in degrees. - + "Longitude of `calendar-location-name' in degrees. The value can be either a decimal fraction (one place of accuracy is sufficient), + east, - west, such as -73.9 for New York City, or the value can be a vector [degrees minutes east/west] such as [73 55 west] for New @@ -116,36 +99,15 @@ York City. This variable should be set in `site-start'.el." :type '(choice (const nil) - (number :tag "Exact") - (vector :value [0 0 west] - (integer :tag "Degrees") - (integer :tag "Minutes") - (choice :tag "Position" - (const east) - (const west)))) + (number :tag "Exact") + (vector :value [0 0 west] + (integer :tag "Degrees") + (integer :tag "Minutes") + (choice :tag "Position" + (const east) + (const west)))) :group 'calendar) -(defsubst calendar-latitude () - "Convert calendar-latitude to a signed decimal fraction, if needed." - (if (numberp calendar-latitude) - calendar-latitude - (let ((lat (+ (aref calendar-latitude 0) - (/ (aref calendar-latitude 1) 60.0)))) - (if (equal (aref calendar-latitude 2) 'north) - lat - (- lat))))) - -(defsubst calendar-longitude () - "Convert calendar-longitude to a signed decimal fraction, if needed." - (if (numberp calendar-longitude) - calendar-longitude - (let ((long (+ (aref calendar-longitude 0) - (/ (aref calendar-longitude 1) 60.0)))) - (if (equal (aref calendar-longitude 2) 'east) - long - (- long))))) - -;;;###autoload (defcustom calendar-location-name '(let ((float-output-format "%.1f")) (format "%s%s, %s%s" @@ -155,24 +117,24 @@ 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")))) - "*Expression evaluating to name of `calendar-longitude', `calendar-latitude'. -For example, \"New York City\". Default value is just the latitude, longitude -pair. + (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'. This variable should be set in `site-start'.el." :type 'sexp :group 'calendar) (defcustom solar-error 0.5 -"*Tolerance (in minutes) for sunrise/sunset calculations. + "Tolerance (in minutes) for sunrise/sunset calculations. A larger value makes the calculations for sunrise/sunset faster, but less accurate. The default is half a minute (30 seconds), so that sunrise/sunset @@ -186,94 +148,125 @@ delta. At present, delta = 0.01 degrees, so the value of the variable :type 'number :group 'calendar) -(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 - "Sidereal time at Greenwich at midnight (universal time).") +(defvar solar-sidereal-time-greenwich-midnight nil + "Sidereal time at Greenwich at midnight (universal time).") (defvar solar-northern-spring-or-summer-season nil "Non-nil if northern spring or summer and nil otherwise. Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.") + +(defsubst calendar-latitude () + "Ensure the variable `calendar-latitude' is a signed decimal fraction." + (if (numberp calendar-latitude) + calendar-latitude + (let ((lat (+ (aref calendar-latitude 0) + (/ (aref calendar-latitude 1) 60.0)))) + (if (eq (aref calendar-latitude 2) 'north) + lat + (- lat))))) + +(defsubst calendar-longitude () + "Ensure the variable `calendar-longitude' is a signed decimal fraction." + (if (numberp calendar-longitude) + calendar-longitude + (let ((long (+ (aref calendar-longitude 0) + (/ (aref calendar-longitude 1) 60.0)))) + (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 user for latitude, longitude, and time zone." + "Prompt for `calendar-longitude', `calendar-latitude', `calendar-time-zone'." (beep) - (if (not calendar-longitude) + (or calendar-longitude (setq calendar-longitude (solar-get-number "Enter longitude (decimal fraction; + east, - west): "))) - (if (not calendar-latitude) + (or calendar-latitude (setq calendar-latitude (solar-get-number "Enter latitude (decimal fraction; + north, - south): "))) - (if (not calendar-time-zone) + (or calendar-time-zone (setq calendar-time-zone (solar-get-number - "Enter difference from Coordinated Universal Time (in minutes): ")))) + "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)))) - -;; The condition-case stuff is needed to catch bogus arithmetic -;; exceptions that occur on some machines (like Sparcs) (defun solar-sin-degrees (x) - (condition-case nil - (sin (degrees-to-radians (mod x 360.0))) - (solar-sin-degrees x))) + "Return sin of X degrees." + (sin (degrees-to-radians (mod x 360.0)))) + (defun solar-cosine-degrees (x) - (condition-case nil - (cos (degrees-to-radians (mod x 360.0))) - (solar-cosine-degrees x))) + "Return cosine of X degrees." + (cos (degrees-to-radians (mod x 360.0)))) + (defun solar-tangent-degrees (x) - (condition-case nil - (tan (degrees-to-radians (mod x 360.0))) - (solar-tangent-degrees x))) + "Return tangent of X degrees." + (tan (degrees-to-radians (mod x 360.0)))) (defun solar-xy-to-quadrant (x y) - "Determines the quadrant of the point X, Y." + "Determine the quadrant of the point X, Y." (if (> x 0) (if (> y 0) 1 4) - (if (> y 0) 2 3))) + (if (> y 0) 2 3))) (defun solar-degrees-to-quadrant (angle) - "Determines the quadrant of ANGLE." + "Determine the quadrant of ANGLE degrees." (1+ (floor (mod angle 360) 90))) (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) - "Arctan of point X, Y." - (if (= x 0) - (if (> y 0) 90 270) - (solar-arctan (/ y x) (solar-xy-to-quadrant x y)))) + "Arctangent of point X, Y." + (if (zerop x) + (if (> y 0) 90 270) + (solar-arctan (/ y x) (solar-xy-to-quadrant x y)))) (defun solar-arccos (x) - "Arcos of X." - (let ((y (sqrt (- 1 (* x x))))) - (solar-atn2 x y))) + "Arccosine of X." + (let ((y (sqrt (- 1 (* x x))))) + (solar-atn2 x y))) (defun solar-arcsin (y) - "Arcsin of Y." - (let ((x (sqrt (- 1 (* y y))))) - (solar-atn2 x y) - )) + "Arcsin of Y." + (let ((x (sqrt (- 1 (* y y))))) + (solar-atn2 x y))) (defsubst solar-degrees-to-hours (degrees) "Convert DEGREES to hours." @@ -298,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 +\(-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. @@ -337,70 +476,110 @@ being 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 +\(-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, +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. Uses binary search." - (let* ((ut (car (cdr time))) - (possible t) ; we assume that rise or set are possible + (let* ((ut (cadr time)) + (possible t) ; we assume that rise or set are possible (utmin (+ ut (* direction 12.0))) - (utmax ut) ; the time searched is between utmin and utmax - ; utmin and utmax are in hours - (utmoment-old 0.0) ; rise or set approximation - (utmoment 1.0) ; rise or set approximation - (hut 0) ; sun height at utmoment + (utmax ut) ; the time searched is between utmin and utmax + ;; utmin and utmax are in hours. + (utmoment-old 0.0) ; rise or set approximation + (utmoment 1.0) ; rise or set approximation + (hut 0) ; sun height at utmoment (t0 (car time)) - (hmin (car (cdr - (solar-horizontal-coordinates (list t0 utmin) - latitude longitude t)))) - (hmax (car (cdr - (solar-horizontal-coordinates (list t0 utmax) - latitude longitude t))))) - ; -0.61 degrees is the height of the middle of the sun, when it rises - ; or sets. - (if (< hmin height) - (if (> hmax height) - (while ;(< i 20) ; we perform a simple dichotomy - ; (> (abs (- hut height)) epsilon) - (>= (abs (- utmoment utmoment-old)) - (/ solar-error 60)) - (setq utmoment-old utmoment) - (setq utmoment (/ (+ utmin utmax) 2)) - (setq hut (car (cdr - (solar-horizontal-coordinates - (list t0 utmoment) latitude longitude t)))) - (if (< hut height) (setq utmin utmoment)) - (if (> hut height) (setq utmax utmoment)) - ) - (setq possible nil)) ; the sun never rises - (setq possible nil)) ; the sun never sets - (if (not possible) nil utmoment))) + (hmin (cadr (solar-horizontal-coordinates (list t0 utmin) + latitude longitude t))) + (hmax (cadr (solar-horizontal-coordinates (list t0 utmax) + latitude longitude t)))) + ;; -0.61 degrees is the height of the middle of the sun, when it + ;; rises or sets. + (if (< hmin height) + (if (> hmax height) + (while ;;; (< i 20) ; we perform a simple dichotomy +;;; (> (abs (- hut height)) epsilon) + (>= (abs (- utmoment utmoment-old)) + (/ solar-error 60)) + (setq utmoment-old utmoment + utmoment (/ (+ utmin utmax) 2) + hut (cadr (solar-horizontal-coordinates + (list t0 utmoment) latitude longitude t))) + (if (< hut height) (setq utmin utmoment)) + (if (> hut height) (setq utmax utmoment))) + (setq possible nil)) ; the sun never rises + (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'." (let* ((time (round (* 60 time))) - (24-hours (/ time 60)) - (minutes (format "%02d" (% time 60))) - (12-hours (format "%d" (1+ (% (+ 24-hours 11) 12)))) - (am-pm (if (>= 24-hours 12) "pm" "am")) - (24-hours (format "%02d" 24-hours))) + (24-hours (/ time 60)) + (minutes (format "%02d" (% time 60))) + (12-hours (format "%d" (1+ (% (+ 24-hours 11) 12)))) + (am-pm (if (>= 24-hours 12) "pm" "am")) + (24-hours (format "%02d" 24-hours))) (mapconcat 'eval calendar-time-display-form ""))) - (defun solar-daylight (time) - "Printable form for time expressed in hours." + "Printable form for TIME expressed in hours." (format "%d:%02d" (floor time) (floor (* 60 (- time (floor time)))))) -(defun solar-exact-local-noon (date) - "Date and Universal Time of local noon at *local date* date. +(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 local date. The second component of date should be an integer." (let* ((nd date) @@ -408,236 +587,86 @@ local date. The second component of date should be an integer." (te (solar-time-equation date ut))) (setq ut (- ut te)) (if (>= ut 24) - (progn - (setq nd (list (car date) (+ 1 (car (cdr date))) - (car (cdr (cdr date))))) - (setq ut (- ut 24)))) + (setq nd (list (car date) (1+ (cadr date)) + (nth 2 date)) + ut (- ut 24))) (if (< ut 0) - (progn - (setq nd (list (car date) (- (car (cdr date)) 1) - (car (cdr (cdr date))))) - (setq ut (+ ut 24)))) - (setq nd (calendar-gregorian-from-absolute - (calendar-absolute-from-gregorian nd))) - ; date standardization + (setq nd (list (car date) (1- (cadr date)) + (nth 2 date)) + ut (+ ut 24))) + (setq nd (calendar-gregorian-from-absolute ; date standardization + (calendar-absolute-from-gregorian nd))) (list nd ut))) +(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." - (let* (; first, get the exact moment of local noon. - (exact-local-noon (solar-exact-local-noon date)) - ; get the time from the 2000 epoch. + ;; 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 + ;; 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))) + (list t0 (cadr 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) + ;; 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)) + (> (nth 2 equator-rise-set) 12)) (solar-sunrise-and-sunset - (list t0 (car (cdr exact-local-noon))) + (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) nil)) - (set (car (cdr rise-set))) - (adj-set (if set (dst-adjust-time date set) nil)) - (length (car (cdr (cdr rise-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 (car (cdr l)) - (concat "sunset " (apply 'solar-time-string (car (cdr l)))) + (if (cadr l) + (concat "sunset " (apply 'solar-time-string (cadr 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 saving 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))) + (if nolocation "" + (format " at %s" (eval calendar-location-name))) + (nth 2 l)))) (defconst solar-data-list '((403406 4.721964 1.621043) @@ -689,29 +718,29 @@ since January 1st, 2000, at 12 ET." (10 2.21 46941.14) (10 3.59 -68.29) (10 1.50 21463.25) - (10 2.55 157208.40))) + (10 2.55 157208.40)) + "Data used for calculation of solar longitude.") (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 +Accuracy 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-astro-to-absolute 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 + ;; 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 @@ -719,227 +748,169 @@ calendar-time-zone are used to interpret local time." (* 62833.1961680 U) (* 0.0000001 (apply '+ - (mapcar '(lambda (x) - (* (car x) - (sin (mod - (+ (car (cdr x)) - (* (car (cdr (cdr x))) U)) - (* 2 pi))))) + (mapcar (lambda (x) + (* (car x) + (sin (mod + (+ (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 of universal time. -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 (car (cdr (cdr (cdr nut-i))))) ; nutation - (i (car (cdr 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." - (let* ((et (solar-date-to-et date ut)) - (ec (solar-ecliptic-coordinates et nil))) - (car (cdr (cdr ec))))) - -(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." - (let ((t0 (solar-julian-ut-centuries date))) - (solar-ephemeris-time (list t0 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. -If called with an optional prefix argument, prompt for date. - -If called with an optional double prefix argument, prompt for longitude, -latitude, time zone, and date, and always use standard time. +If called with an optional prefix argument ARG, prompt for date. +If called with an optional double prefix argument, prompt for +longitude, latitude, time zone, and date, and always use standard time. This function is suitable for execution in a .emacs file." - (interactive "p") - (or arg (setq arg 1)) - (if (and (< arg 16) - (not (and calendar-latitude calendar-longitude calendar-time-zone))) - (solar-setup)) - (let* ((calendar-longitude - (if (< arg 16) calendar-longitude - (solar-get-number - "Enter longitude (decimal fraction; + east, - west): "))) - (calendar-latitude - (if (< arg 16) calendar-latitude - (solar-get-number - "Enter latitude (decimal fraction; + north, - south): "))) - (calendar-time-zone - (if (< arg 16) calendar-time-zone - (solar-get-number - "Enter difference from Coordinated Universal Time (in minutes): "))) - (calendar-location-name - (if (< arg 16) calendar-location-name - (let ((float-output-format "%.1f")) - (format "%s%s, %s%s" - (if (numberp calendar-latitude) - (abs calendar-latitude) - (+ (aref calendar-latitude 0) - (/ (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 (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")))))) - (calendar-standard-time-zone-name - (if (< arg 16) calendar-standard-time-zone-name - (cond ((= calendar-time-zone 0) "UTC") - ((< calendar-time-zone 0) - (format "UTC%dmin" calendar-time-zone)) - (t (format "UTC+%dmin" calendar-time-zone))))) - (calendar-daylight-savings-starts - (if (< arg 16) calendar-daylight-savings-starts)) - (calendar-daylight-savings-ends - (if (< arg 16) calendar-daylight-savings-ends)) - (date (if (< arg 4) (calendar-current-date) (calendar-read-date))) - (date-string (calendar-date-string date t)) - (time-string (solar-sunrise-sunset-string date)) - (msg (format "%s: %s" date-string time-string)) - (one-window (one-window-p t))) - (if (<= (length msg) (frame-width)) - (message "%s" msg) - (with-output-to-temp-buffer "*temp*" - (princ (concat date-string "\n" time-string))) - (message "%s" - (substitute-command-keys - (if one-window - (if pop-up-windows - "Type \\[delete-other-windows] to remove temp window." - "Type \\[switch-to-buffer] RET to remove temp window.") - "Type \\[switch-to-buffer-other-window] RET to restore old contents of temp window.")))))) - -(defun calendar-sunrise-sunset () + (interactive "p") + (or arg (setq arg 1)) + (if (and (< arg 16) + (not (and calendar-latitude calendar-longitude calendar-time-zone))) + (solar-setup)) + (let* ((calendar-longitude + (if (< arg 16) calendar-longitude + (solar-get-number + "Enter longitude (decimal fraction; + east, - west): "))) + (calendar-latitude + (if (< arg 16) calendar-latitude + (solar-get-number + "Enter latitude (decimal fraction; + north, - south): "))) + (calendar-time-zone + (if (< arg 16) calendar-time-zone + (solar-get-number + "Enter difference from Coordinated Universal Time (in minutes): "))) + (calendar-location-name + (if (< arg 16) calendar-location-name + (let ((float-output-format "%.1f")) + (format "%s%s, %s%s" + (if (numberp calendar-latitude) + (abs calendar-latitude) + (+ (aref calendar-latitude 0) + (/ (aref calendar-latitude 1) 60.0))) + (if (numberp calendar-latitude) + (if (> calendar-latitude 0) "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 (eq (aref calendar-longitude 2) 'east) + "E" "W")))))) + (calendar-standard-time-zone-name + (if (< arg 16) calendar-standard-time-zone-name + (cond ((zerop calendar-time-zone) "UTC") + ((< calendar-time-zone 0) + (format "UTC%dmin" calendar-time-zone)) + (t (format "UTC+%dmin" calendar-time-zone))))) + (calendar-daylight-savings-starts + (if (< arg 16) calendar-daylight-savings-starts)) + (calendar-daylight-savings-ends + (if (< arg 16) calendar-daylight-savings-ends)) + (date (if (< arg 4) (calendar-current-date) (calendar-read-date))) + (date-string (calendar-date-string date t)) + (time-string (solar-sunrise-sunset-string date)) + (msg (format "%s: %s" date-string time-string)) + (one-window (one-window-p t))) + (if (<= (length msg) (frame-width)) + (message "%s" msg) + (with-output-to-temp-buffer "*temp*" + (princ (concat date-string "\n" time-string))) + (message "%s" + (substitute-command-keys + (if one-window + (if pop-up-windows + "Type \\[delete-other-windows] to remove temp window." + "Type \\[switch-to-buffer] RET to remove temp window.") + "Type \\[switch-to-buffer-other-window] RET to restore old \ +contents of temp window.")))))) + +;;;###cal-autoload +(defun calendar-sunrise-sunset (&optional event) "Local time of sunrise and sunset for date under cursor. Accurate to a few seconds." - (interactive) - (if (not (and calendar-latitude calendar-longitude calendar-time-zone)) + (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 diary-list-sexp-entries, where DATE is bound. +;;;###diary-autoload (defun diary-sunrise-sunset () "Local time of sunrise and sunset as a diary entry. Accurate to a few seconds." - (if (not (and calendar-latitude calendar-longitude calendar-time-zone)) + (or (and calendar-latitude calendar-longitude calendar-time-zone) (solar-setup)) (solar-sunrise-sunset-string date)) -(defcustom diary-sabbath-candles-minutes 18 - "*Number of minutes before sunset for sabbath candle lighting." - :group 'diary - :type 'integer - :version "21.1") - -(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." - (if (not (and calendar-latitude calendar-longitude calendar-time-zone)) - (solar-setup)) - (if (= (% (calendar-absolute-from-gregorian date) 7) 5);; Friday - (let* ((sunset (car (cdr (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 +;; From Meeus, 1991, page 167. (defconst solar-seasons-data '((485 324.96 1934.136) (203 337.23 32964.467) @@ -964,134 +935,135 @@ use when highlighting the day in the calendar." (12 287.11 31931.756) (12 320.81 34777.259) (9 227.73 1222.114) - (8 15.45 16859.074))) + (8 15.45 16859.074)) + "Data for solar equinox/solstice calculations.") (defun solar-equinoxes/solstices (k year) "Date of equinox/solstice K for YEAR. K=0, spring equinox; K=1, summer solstice; K=2, fall equinox; -K=3, winter solstice. -RESULT is a gregorian local date. - -Accurate to less than a minute between 1951 and 2050." +K=3, winter solstice. RESULT is a Gregorian local date. +Accurate to within a minute between 1951 and 2050." (let* ((JDE0 (solar-mean-equinoxes/solstices k year)) (T (/ (- JDE0 2451545.0) 36525)) (W (- (* 35999.373 T) 2.47)) (Delta-lambda (+ 1 (* 0.0334 (solar-cosine-degrees W)) - (* 0.0007 (solar-cosine-degrees (* 2 W))))) - (S (apply '+ (mapcar '(lambda(x) - (* (car x) (solar-cosine-degrees - (+ (* (car (cdr (cdr x))) T) - (car (cdr x)))))) + (* 0.0007 (solar-cosine-degrees (* 2 W))))) + (S (apply '+ (mapcar (lambda(x) + (* (car x) (solar-cosine-degrees + (+ (* (nth 2 x) T) (cadr x))))) solar-seasons-data))) (JDE (+ JDE0 (/ (* 0.00001 S) Delta-lambda))) + ;; Ephemeris time correction. (correction (+ 102.3 (* 123.5 T) (* 32.5 T T))) - ; ephemeris time correction (JD (- JDE (/ correction 86400))) (date (calendar-gregorian-from-absolute (floor (- JD 1721424.5)))) - (time (- (- JD 0.5) (floor (- JD 0.5)))) - ) - (list (car date) (+ (car (cdr date)) time - (/ (/ calendar-time-zone 60.0) 24.0)) - (car (cdr (cdr date)))))) + (time (- (- JD 0.5) (floor (- JD 0.5))))) + (list (car date) (+ (cadr date) time + (/ (/ calendar-time-zone 60.0) 24.0)) + (nth 2 date)))) -; from Meeus, 1991, page 166 +;; From Meeus, 1991, page 166. (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)))) - ; 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))))))) + (if (< year 1000) ; actually between -1000 and 1000 + (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 ((= 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) -;;;###autoload +;;;###holiday-autoload (defun solar-equinoxes-solstices () - "*local* date and time of equinoxes and solstices, if visible in the calendar window. + "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 (car (car adj)) - (+ (car (cdr (car adj)) ) - (/ (car (cdr adj)) 24.0)) - (car (cdr (cdr (car 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) -;;; arch-tag: bc0ff693-df58-4666-bde4-2a7837ccb8fe +;; arch-tag: bc0ff693-df58-4666-bde4-2a7837ccb8fe ;;; solar.el ends here