;; Copyright (C) 1992, 1993, 1995 Free Software Foundation, Inc.
;; Author: Edward M. Reingold <reingold@cs.uiuc.edu>
+;; Denis B. Roegel <Denis.Roegel@loria.fr>
;; Keywords: calendar
;; Human-Keywords: sunrise, sunset, equinox, solstice, calendar, diary,
;; holidays
;; Based on the ``Almanac for Computers 1984,'' prepared by the Nautical
;; Almanac Office, United States Naval Observatory, Washington, 1984, on
;; ``Astronomical Formulae for Calculators,'' 3rd ed., by Jean Meeus,
-;; Willmann-Bell, Inc., 1985, and on ``Astronomical Algorithms'' by Jean
-;; Meeus, Willmann-Bell, Inc., 1991.
+;; Willmann-Bell, Inc., 1985, on ``Astronomical Algorithms'' by Jean Meeus,
+;; Willmann-Bell, Inc., 1991, and on ``Planetary Programs and Tables from
+;; -4000 to +2800'' by Pierre Bretagnon and Jean-Louis Simon, Willmann-Bell,
+;; Inc., 1986.
;;
-;; WARNINGS:
-;; 1. SUNRISE/SUNSET calculations will be accurate only to +/- 2 minutes.
-;; Locations should be between +/- 65 degrees of latitude.
-;; Dates should be in the latter half of the 20th century.
+;; Accuracy:
+;; 1. Sunrise/sunset times will be accurate to the minute for years
+;; 1951--2050. For other years the times will be within +/- 2 minutes.
;;
-;; 2. Equinox/solstice times will be accurate only to +/- 15 minutes.
-
-;; The author would be delighted to have an astronomically more sophisticated
-;; person rewrite the code for the solar calculations in this file!
+;; 2. Equinox/solstice times will be accurate to the minute for years
+;; 1951--2050. For other years the times will be within +/- 1 minute.
;; Comments, corrections, and improvements should be sent to
;; Edward M. Reingold Department of Computer Science
(error "Solar/lunar calculations impossible since floating point is unavailable."))
(require 'cal-dst)
+(require 'cal-julian)
;;;###autoload
(defvar calendar-time-display-form
can be a vector [degrees minutes east/west] such as [73 55 west] for New
York City.
-This variable should be set in site-local.el.")
+This variable should be set in site-start.el.")
(defsubst calendar-latitude ()
"Convert calendar-latitude to a signed decimal fraction, if needed."
For example, \"New York City\". Default value is just the latitude, longitude
pair.
-This variable should be set in site-local.el.")
+This variable should be set in site-start.el.")
+
+(defvar solar-error 0.5
+"*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
+times will be correct to the minute.
+
+It is useless to set the value smaller than 4*delta, where delta is the
+accuracy in the longitude of the sun (given by the function
+`solar-ecliptic-coordinates') in degrees since (delta/360) x (86400/60) = 4 x
+delta. At present, delta = 0.01 degrees, so the value of the variable
+`solar-error' should be at least 0.04 minutes (about 2.5 seconds).")
(defvar solar-n-hemi-seasons
'("Vernal Equinox" "Summer Solstice" "Autumnal Equinox" "Winter Solstice")
'("Autumnal Equinox" "Winter Solstice" "Vernal Equinox" "Summer Solstice")
"List of season changes for the southern hemisphere.")
+(defvar solar-sidereal-time-greenwich-midnight
+ nil
+ "Sidereal time at Greenwich at midnight (universal time).")
+
+(defvar solar-spring-or-summer-season nil
+ "T if spring or summer and nil otherwise.
+Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.")
+
(defun solar-setup ()
"Prompt user for latitude, longitude, and time zone."
(beep)
(if (not (string-equal x ""))
(string-to-int x))))
-(defsubst solar-sin-degrees (x)
- (sin (degrees-to-radians (mod x 360.0))))
-
-(defsubst solar-cosine-degrees (x)
- (cos (degrees-to-radians (mod x 360.0))))
-
-(defsubst solar-tangent-degrees (x)
- (tan (degrees-to-radians (mod x 360.0))))
-
+;; 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)))
+(defun solar-cosine-degrees (x)
+ (condition-case nil
+ (cos (degrees-to-radians (mod x 360.0)))
+ (solar-cosine-degrees x)))
+(defun solar-tangent-degrees (x)
+ (condition-case nil
+ (tan (degrees-to-radians (mod x 360.0)))
+ (solar-tangent-degrees x)))
+
(defun solar-xy-to-quadrant (x y)
"Determines the quadrant of the point X, Y."
(if (> x 0)
((equal quad 4) (+ deg 360))
(t deg))))
+(defun solar-atn2 (x y)
+ "Arctan of point X, Y."
+ (if (= y 0)
+ (if (> x 0) 90 270)
+ (solar-arctan (/ x y) y)))
+
(defun solar-arccos (x)
+ "Arcos of X."
(let ((y (sqrt (- 1 (* x x)))))
(solar-arctan (/ y x) (solar-xy-to-quadrant x y))))
(defun solar-arcsin (y)
+ "Arcsin of Y."
(let ((x (sqrt (- 1 (* y y)))))
(solar-arctan (/ y x) (solar-xy-to-quadrant x y))))
-(defconst solar-earth-inclination 23.441884
- "Inclination of earth's equator to its solar orbit in degrees.")
-
-(defconst solar-cos-inclination (solar-cosine-degrees solar-earth-inclination)
- "Cosine of earth's inclination.")
-
-(defconst solar-sin-inclination (solar-sin-degrees solar-earth-inclination)
- "Sine of earth's inclination.")
-
-(defconst solar-earth-orbit-eccentricity 0.016718
- "Eccentricity of orbit of the earth around the sun.")
-
-(defsubst solar-degrees-to-hours (deg)
- (/ deg 15.0))
+(defsubst solar-degrees-to-hours (degrees)
+ "Convert DEGREES to hours."
+ (/ degrees 15.0))
(defsubst solar-hours-to-days (hour)
+ "Convert HOUR to decimal fraction of a day."
(/ hour 24.0))
-(defun solar-longitude-of-sun (day)
- "Longitude of the sun at DAY in the year."
- (let ((mean-anomaly (- (* 0.9856 day) 3.289)))
- (mod (+ mean-anomaly
- (* 1.916 (solar-sin-degrees mean-anomaly))
- (* 0.020 (solar-sin-degrees (* 2 mean-anomaly)))
- 282.634)
- 360)))
-
-(defun solar-right-ascension (longitude)
- "Right ascension of the sun, given its LONGITUDE."
+(defun solar-right-ascension (longitude obliquity)
+ "Right ascension of the sun, in hours, given LONGITUDE and OBLIQUITY.
+Both arguments are in degrees."
(solar-degrees-to-hours
(solar-arctan
- (* solar-cos-inclination (solar-tangent-degrees longitude))
+ (* (solar-cosine-degrees obliquity) (solar-tangent-degrees longitude))
(solar-degrees-to-quadrant longitude))))
-(defun solar-declination (longitude)
- "Declination of the sun, given its LONGITUDE."
+(defun solar-declination (longitude obliquity)
+ "Declination of the sun, in degrees, given LONGITUDE and OBLIQUITY.
+Both arguments are in degrees."
(solar-arcsin
- (* solar-sin-inclination
+ (* (solar-sin-degrees obliquity)
(solar-sin-degrees longitude))))
-(defun solar-sunrise (date)
- "Calculates the *standard* time of sunrise for Gregorian DATE.
-Calculation is for location given by `calendar-latitude' and
-`calendar-longitude'.
-
-Returns a decimal fraction of hours. Returns nil if the sun does not rise at
-that location on that day."
- (let* ((day-of-year (calendar-day-number date))
- (approx-sunrise
- (+ day-of-year
- (solar-hours-to-days
- (- 6 (solar-degrees-to-hours (calendar-longitude))))))
- (solar-longitude-of-sun-at-sunrise
- (solar-longitude-of-sun approx-sunrise))
- (solar-right-ascension-at-sunrise
- (solar-right-ascension solar-longitude-of-sun-at-sunrise))
- (solar-declination-at-sunrise
- (solar-declination solar-longitude-of-sun-at-sunrise))
- (cos-local-sunrise
- (/ (- (solar-cosine-degrees (+ 90 (/ 50.0 60.0)))
- (* (solar-sin-degrees solar-declination-at-sunrise)
- (solar-sin-degrees (calendar-latitude))))
- (* (solar-cosine-degrees solar-declination-at-sunrise)
- (solar-cosine-degrees (calendar-latitude))))))
- (if (<= (abs cos-local-sunrise) 1);; otherwise, no sunrise that day
- (let* ((local-sunrise (solar-degrees-to-hours
- (- 360 (solar-arccos cos-local-sunrise))))
- (local-mean-sunrise
- (mod (- (+ local-sunrise solar-right-ascension-at-sunrise)
- (+ (* 0.065710 approx-sunrise)
- 6.622))
- 24)))
- (+ (- local-mean-sunrise (solar-degrees-to-hours (calendar-longitude)))
- (/ calendar-time-zone 60.0))))))
-
-(defun solar-sunset (date)
- "Calculates the *standard* time of sunset for Gregorian DATE.
-Calculation is for location given by `calendar-latitude' and
-`calendar-longitude'.
-
-Returns a decimal fractions of hours. Returns nil if the sun does not set at
-that location on that day."
- (let* ((day-of-year (calendar-day-number date))
- (approx-sunset
- (+ day-of-year
- (solar-hours-to-days
- (- 18 (solar-degrees-to-hours (calendar-longitude))))))
- (solar-longitude-of-sun-at-sunset
- (solar-longitude-of-sun approx-sunset))
- (solar-right-ascension-at-sunset
- (solar-right-ascension solar-longitude-of-sun-at-sunset))
- (solar-declination-at-sunset
- (solar-declination solar-longitude-of-sun-at-sunset))
- (cos-local-sunset
- (/ (- (solar-cosine-degrees (+ 90 (/ 50.0 60.0)))
- (* (solar-sin-degrees solar-declination-at-sunset)
- (solar-sin-degrees (calendar-latitude))))
- (* (solar-cosine-degrees solar-declination-at-sunset)
- (solar-cosine-degrees (calendar-latitude))))))
- (if (<= (abs cos-local-sunset) 1);; otherwise, no sunset that day
- (let* ((local-sunset (solar-degrees-to-hours
- (solar-arccos cos-local-sunset)))
- (local-mean-sunset
- (mod (- (+ local-sunset solar-right-ascension-at-sunset)
- (+ (* 0.065710 approx-sunset) 6.622))
- 24)))
- (+ (- local-mean-sunset (solar-degrees-to-hours (calendar-longitude)))
- (/ calendar-time-zone 60.0))))))
+(defun solar-sunrise-and-sunset (time latitude longitude)
+ "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.
+
+Coordinates are included because this function is called with latitude=10
+degrees to find out if polar regions have 24 hours of sun or only night."
+ (let* ((rise-time (solar-moment -1 latitude longitude time))
+ (set-time (solar-moment 1 latitude longitude time))
+ (day-length))
+ (if (not (and rise-time set-time))
+ (if (or (and (> latitude 0) solar-spring-or-summer-season)
+ (and (< latitude 0) (not solar-spring-or-summer-season)))
+ (setq day-length 24)
+ (setq day-length 0))
+ (setq day-length (- set-time rise-time)))
+ (list (+ rise-time (/ calendar-time-zone 60.0))
+ (+ set-time (/ calendar-time-zone 60.0)) day-length)))
+
+(defun solar-moment (direction latitude longitude time)
+ "Sunrise/sunset at location.
+Sunrise if DIRECTION =-1 or sunset if =1 at LATITUDE, LONGITUDE, with midday
+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
+Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT.
+
+Uses binary search."
+ (let* ((ut (car (cdr time)))
+ (possible 1) ; 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
+ (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 -0.61)
+ (if (> hmax -0.61)
+ (while ;(< i 20) ; we perform a simple dichotomy
+ ; (> (abs (+ hut 0.61)) 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 -0.61) (setq utmin utmoment))
+ (if (> hut -0.61) (setq utmax utmoment))
+ )
+ (setq possible 0)) the sun never rises
+ (setq possible 0)) ; the sun never sets
+ (if (equal possible 0) nil utmoment)))
(defun solar-time-string (time time-zone)
"Printable form for decimal fraction TIME in TIME-ZONE.
(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-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)
+ (ut (- 12.0 (/ calendar-longitude 15)))
+ (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))))
+ (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
+ (list nd ut)))
+
(defun solar-sunrise-sunset (date)
- "String giving local times of sunrise and sunset on Gregorian DATE."
- (let* ((rise (solar-sunrise date))
- (adj-rise (if rise (dst-adjust-time date rise)))
- (set (solar-sunset date))
- (adj-set (if set (dst-adjust-time date set))))
- (format "%s, %s at %s"
- (if (and rise (calendar-date-equal date (car adj-rise)))
- (concat "Sunrise " (apply 'solar-time-string (cdr adj-rise)))
- "No sunrise")
- (if (and set (calendar-date-equal date (car adj-set)))
- (concat "sunset " (apply 'solar-time-string (cdr adj-set)))
- "no sunset")
- (eval calendar-location-name))))
+ "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 the time from the 2000 epoch.
+ (t0 (solar-julian-ut-centuries (car exact-local-noon)))
+ ; store the sidereal time at Greenwich at midnight of UT time.
+ ; find if summer or winter slightly above the equator
+ (equator-rise-set
+ (progn (setq solar-sidereal-time-greenwich-midnight
+ (solar-sidereal-time t0))
+ (solar-sunrise-and-sunset
+ (list t0 (car (cdr exact-local-noon)))
+ 10.0
+ calendar-longitude)))
+ ; 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-spring-or-summer-season
+ (if (> (car (cdr (cdr equator-rise-set))) 12) 1 0))
+ (solar-sunrise-and-sunset
+ (list t0 (car (cdr exact-local-noon)))
+ calendar-latitude
+ calendar-longitude)))
+ (rise (car rise-set))
+ (adj-rise (if rise (dst-adjust-time date rise) nil))
+ (set (car (cdr rise-set)))
+ (adj-set (if set (dst-adjust-time date set) nil))
+ (length (car (cdr (cdr rise-set)))) )
+ (list
+ (and rise (calendar-date-equal date (car adj-rise)) (cdr adj-rise))
+ (and set (calendar-date-equal date (car adj-set)) (cdr adj-set))
+ (solar-daylight length))))
+
+(defun solar-sunrise-sunset-string (date)
+ "String of *local* times of sunrise, sunset, and daylight on Gregorian DATE."
+ (let ((l (solar-sunrise-sunset date)))
+ (format
+ "%s, %s at %s (%s hours daylight)"
+ (if (car l)
+ (concat "Sunrise " (apply 'solar-time-string (car l)))
+ "No sunrise")
+ (if (car (cdr l))
+ (concat "sunset " (apply 'solar-time-string (car (cdr l))))
+ "no sunset")
+ (eval calendar-location-name)
+ (car (cdr (cdr l))))))
+
+(defun solar-julian-ut-centuries (date)
+ "Number of Julian centuries elapsed since 1 Jan, 2000 at noon U.T. for Gregorian DATE."
+ (/ (- (calendar-absolute-from-gregorian date)
+ (calendar-absolute-from-gregorian '(1 1.5 2000)))
+ 36525.0))
+
+(defun solar-ephemeris-time(time)
+ "Ephemeris Time at moment TIME.
+
+TIME is a pair with the first component being the number of Julian centuries
+elapsed at 0 Universal Time, and the second component being the universal
+time. For instance, the pair corresponding to November 28, 1995 at 16 UT is
+(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
+Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT.
+
+Result is in julian centuries of ephemeris time."
+ (let* ((t0 (car time))
+ (ut (car (cdr time)))
+ (t1 (+ t0 (/ (/ ut 24.0) 36525)))
+ (y (+ 2000 (* 100 t1)))
+ (dt (* 86400 (solar-ephemeris-correction (floor y)))))
+ (+ t1 (/ (/ dt 86400) 36525))))
(defun solar-date-next-longitude (d l)
"First moment on or after Julian day number D when sun's longitude is a
(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-sin-degrees ah)
+ (- (* (solar-cosine-degrees ah)
+ (solar-sin-degrees latitude))
+ (* (solar-tangent-degrees de)
+ (solar-cosine-degrees latitude)))))
+ (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)))
+
(defun solar-longitude (d)
"Longitude of sun on astronomical (Julian) day number D.
-Accuracy is about 0.01 degree (about 365.25*24*60*0.01/360 = 15 minutes).
+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))
- (date (calendar-gregorian-from-absolute (floor a-d)))
- (time (* 24 (- a-d (truncate a-d))))
- (rounded-abs-date (+ (calendar-absolute-from-gregorian date)
- (/ (round (* 60 time)) 60.0 24.0)))
- ;; get local standard time
- (a-d (- rounded-abs-date
- (if (dst-in-effect rounded-abs-date)
- (/ calendar-daylight-time-offset 24.0 60.0) 0)))
;; get Universal Time
- (a-d (- a-d (/ calendar-time-zone 60.0 24.0)))
- (date (calendar-astro-from-absolute a-d))
+ (date (calendar-astro-from-absolute
+ (- a-d
+ (if (dst-in-effect a-d)
+ (/ calendar-daylight-time-offset 24.0 60.0) 0)
+ (/ calendar-time-zone 60.0 24.0))))
;; get Ephemeris Time
(date (+ date (solar-ephemeris-correction
(extract-calendar-year
(floor
(calendar-absolute-from-astro
date)))))))
- (T (/ (- date 2451545.0) 36525.0))
- (Lo (mod (+ 280.46645 (* 36000.76983 T) (* 0.0003032 T T)) 360.0))
- (M (mod (+ 357.52910
- (* 35999.05030 T)
- (* -0.0001559 T T)
- (* -0.00000048 T T T))
- 360.0))
- (e (+ 0.016708617 (* -0.000042037 T) (* -0.0000001236 T T)))
- (C (+ (* (+ 1.914600 (* -0.004817 T) (* -0.000014 T T))
- (solar-sin-degrees M))
- (* (+ 0.019993 (* -0.000101 T)) (solar-sin-degrees (* 2 M)))
- (* 0.000290 (solar-sin-degrees (* 3 M)))))
- (true-longitude (+ Lo C))
- (omega (+ 125.04 (* -1934.136 T)))
- (apparent-longitude (mod
- (+ true-longitude
- -0.00569
- (* -0.00478 (solar-sin-degrees omega)))
- 360.0)))
- apparent-longitude))
+ (U (/ (- date 2451545) 3652500))
+ (longitude
+ (+ 4.9353929
+ (* 62833.1961680 U)
+ (* 0.0000001
+ (apply '+
+ (mapcar '(lambda (x)
+ (* (car x)
+ (sin (mod
+ (+ (car (cdr x))
+ (* (car (cdr (cdr x))) U))
+ (* 2 pi)))))
+ solar-data-list)))))
+ (aberration
+ (* 0.0000001 (- (* 17 (cos (+ 3.10 (* 62830.14 U)))) 973)))
+ (A1 (mod (+ 2.18 (* U (+ -3375.70 (* 0.36 U)))) (* 2 pi)))
+ (A2 (mod (+ 3.51 (* U (+ 125666.39 (* 0.10 U)))) (* 2 pi)))
+ (nutation (* -0.0000001 (+ (* 834 (sin A1)) (* 64 (sin A2))))))
+ (mod (radians-to-degrees (+ longitude aberration nutation)) 360.0)))
+
+(defconst solar-data-list
+ '((403406 4.721964 1.621043)
+ (195207 5.937458 62830.348067)
+ (119433 1.115589 62830.821524)
+ (112392 5.781616 62829.634302)
+ (3891 5.5474 125660.5691)
+ (2819 1.5120 125660.984)
+ (1721 4.1897 62832.4766)
+ (0 1.163 0.813)
+ (660 5.415 125659.31)
+ (350 4.315 57533.85)
+ (334 4.553 -33.931)
+ (314 5.198 777137.715)
+ (268 5.989 78604.191)
+ (242 2.911 5.412)
+ (234 1.423 39302.098)
+ (158 0.061 -34.861)
+ (132 2.317 115067.698)
+ (129 3.193 15774.337)
+ (114 2.828 5296.670)
+ (99 0.52 58849.27)
+ (93 4.65 5296.11)
+ (86 4.35 -3980.70)
+ (78 2.75 52237.69)
+ (72 4.50 55076.47)
+ (68 3.23 261.08)
+ (64 1.22 15773.85)
+ (46 0.14 188491.03)
+ (38 3.44 -7756.55)
+ (37 4.37 264.89)
+ (32 1.14 117906.27)
+ (29 2.84 55075.75)
+ (28 5.96 -7961.39)
+ (27 5.09 188489.81)
+ (27 1.72 2132.19)
+ (25 2.56 109771.03)
+ (24 1.92 54868.56)
+ (21 0.09 25443.93)
+ (21 5.98 -55731.43)
+ (20 4.03 60697.74)
+ (18 4.47 2132.79)
+ (17 0.79 109771.63)
+ (14 4.24 -7752.82)
+ (13 2.01 188491.91)
+ (13 2.65 207.81)
+ (13 4.98 29424.63)
+ (12 0.93 -7.99)
+ (10 2.21 46941.14)
+ (10 3.59 -68.29)
+ (10 1.50 21463.25)
+ (10 2.55 157208.40)))
(defun solar-ephemeris-correction (year)
- "Ephemeris time minus Universal Time at astronomical (Julian) day D.
-Result is in days For the years 1800-1987, the maximum error is 1.9 seconds.
+ "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))
(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))))
+
;;;###autoload
(defun sunrise-sunset (&optional arg)
- "Local time of sunrise and sunset for today. Accurate to +/- 2 minutes.
+ "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,
(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 date))
+ (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))
(defun calendar-sunrise-sunset ()
"Local time of sunrise and sunset for date under cursor.
-Accurate to +/- 2 minutes."
+Accurate to a few seconds."
(interactive)
(if (not (and calendar-latitude calendar-longitude calendar-time-zone))
(solar-setup))
(let ((date (calendar-cursor-to-date t)))
(message "%s: %s"
(calendar-date-string date t t)
- (solar-sunrise-sunset date))))
+ (solar-sunrise-sunset-string date))))
(defun diary-sunrise-sunset ()
"Local time of sunrise and sunset as a diary entry.
-Accurate to +/- 2 minutes."
+Accurate to a few seconds."
(if (not (and calendar-latitude calendar-longitude calendar-time-zone))
(solar-setup))
- (solar-sunrise-sunset date))
+ (solar-sunrise-sunset-string date))
(defun diary-sabbath-candles ()
"Local time of candle lighting diary entry--applies if date is a Friday.
(if (not (and calendar-latitude calendar-longitude calendar-time-zone))
(solar-setup))
(if (= (% (calendar-absolute-from-gregorian date) 7) 5);; Friday
- (let* ((sunset (solar-sunset date))
+ (let* ((sunset (car (cdr (solar-sunrise-sunset date))))
(light (if sunset
(dst-adjust-time
date
- (- sunset (/ 18.0 60.0))))))
+ (- (car sunset) (/ 18.0 60.0))))))
(if (and light (calendar-date-equal date (car light)))
(format "%s Sabbath candle lighting"
(apply 'solar-time-string (cdr light)))))))
+(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."
+ (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))))))
+ solar-seasons-data)))
+ (JDE (+ JDE0 (/ (* 0.00001 S) Delta-lambda)))
+ (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))))))
+
+; 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."
+ (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)))))))
+
+; from Meeus, 1991, page 167
+(defconst solar-seasons-data
+ '((485 324.96 1934.136)
+ (203 337.23 32964.467)
+ (199 342.08 20.186)
+ (182 27.85 445267.112)
+ (156 73.14 45036.886)
+ (136 171.52 22518.443)
+ (77 222.54 65928.934)
+ (74 296.72 3034.906)
+ (70 243.58 9037.513)
+ (58 119.81 33718.147)
+ (52 297.17 150.678)
+ (50 21.02 2281.226)
+ (45 247.54 29929.562)
+ (44 325.15 31555.956)
+ (29 60.93 4443.417)
+ (18 155.12 67555.328)
+ (17 288.79 4562.452)
+ (16 198.04 62894.029)
+ (14 199.76 31436.921)
+ (12 95.39 14577.848)
+ (12 287.11 31931.756)
+ (12 320.81 34777.259)
+ (9 227.73 1222.114)
+ (8 15.45 16859.074)))
+
;;;###autoload
(defun solar-equinoxes-solstices ()
- "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 window.
Requires floating point."
(let ((m displayed-month)
(y displayed-year))
(if calendar-time-zone calendar-daylight-savings-ends))
(calendar-time-zone (if calendar-time-zone calendar-time-zone 0))
(k (1- (/ m 3)))
- (d (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)))
+ (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 d1) (+ (car (cdr d1))
+ (/ (car (cdr adj)) 24.0))
+ (car (cdr (cdr d1)))))
+ ; 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"
HCoop / bpt / emacs.git / blobdiff