1 ;;; solar.el --- calendar functions for solar events.
3 ;; Copyright (C) 1992, 1993, 1995 Free Software Foundation, Inc.
5 ;; Author: Edward M. Reingold <reingold@cs.uiuc.edu>
7 ;; Human-Keywords: sunrise, sunset, equinox, solstice, calendar, diary,
10 ;; This file is part of GNU Emacs.
12 ;; GNU Emacs is free software; you can redistribute it and/or modify
13 ;; it under the terms of the GNU General Public License as published by
14 ;; the Free Software Foundation; either version 2, or (at your option)
17 ;; GNU Emacs is distributed in the hope that it will be useful,
18 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
19 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 ;; GNU General Public License for more details.
22 ;; You should have received a copy of the GNU General Public License
23 ;; along with GNU Emacs; see the file COPYING. If not, write to
24 ;; the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
28 ;; This collection of functions implements the features of calendar.el,
29 ;; diary.el, and holiday.el that deal with times of day, sunrise/sunset, and
30 ;; eqinoxes/solstices.
32 ;; Based on the ``Almanac for Computers 1984,'' prepared by the Nautical
33 ;; Almanac Office, United States Naval Observatory, Washington, 1984, on
34 ;; ``Astronomical Formulae for Calculators,'' 3rd ed., by Jean Meeus,
35 ;; Willmann-Bell, Inc., 1985, and on ``Astronomical Algorithms'' by Jean
36 ;; Meeus, Willmann-Bell, Inc., 1991.
40 ;; 1. SUNRISE/SUNSET calculations will be accurate only to +/- 2 minutes.
41 ;; Locations should be between +/- 65 degrees of latitude.
42 ;; Dates should be in the latter half of the 20th century.
44 ;; 2. Equinox/solstice times will be accurate only to +/- 15 minutes.
46 ;; The author would be delighted to have an astronomically more sophisticated
47 ;; person rewrite the code for the solar calculations in this file!
49 ;; Comments, corrections, and improvements should be sent to
50 ;; Edward M. Reingold Department of Computer Science
51 ;; (217) 333-6733 University of Illinois at Urbana-Champaign
52 ;; reingold@cs.uiuc.edu 1304 West Springfield Avenue
53 ;; Urbana, Illinois 61801
58 (require 'lisp-float-type)
59 (error "Solar/lunar calculations impossible since floating point is unavailable."))
64 (defvar calendar-time-display-form
65 '(12-hours ":" minutes am-pm
66 (if time-zone " (") time-zone (if time-zone ")"))
67 "*The pseudo-pattern that governs the way a time of day is formatted.
69 A pseudo-pattern is a list of expressions that can involve the keywords
70 `12-hours', `24-hours', and `minutes', all numbers in string form,
71 and `am-pm' and `time-zone', both alphabetic strings.
75 '(24-hours \":\" minutes
76 (if time-zone \" (\") time-zone (if time-zone \")\"))
78 would give military-style times like `21:07 (UTC)'.")
81 (defvar calendar-latitude nil
82 "*Latitude of `calendar-location-name' in degrees.
84 The value can be either a decimal fraction (one place of accuracy is
85 sufficient), + north, - south, such as 40.7 for New York City, or the value
86 can be a vector [degrees minutes north/south] such as [40 50 north] for New
89 This variable should be set in site-local.el.")
92 (defvar calendar-longitude nil
93 "*Longitude of `calendar-location-name' in degrees.
95 The value can be either a decimal fraction (one place of accuracy is
96 sufficient), + east, - west, such as -73.9 for New York City, or the value
97 can be a vector [degrees minutes east/west] such as [73 55 west] for New
100 This variable should be set in site-local.el.")
102 (defsubst calendar-latitude ()
103 "Convert calendar-latitude to a signed decimal fraction, if needed."
104 (if (numberp calendar-latitude)
106 (let ((lat (+ (aref calendar-latitude 0)
107 (/ (aref calendar-latitude 1) 60.0))))
108 (if (equal (aref calendar-latitude 2) 'north)
112 (defsubst calendar-longitude ()
113 "Convert calendar-longitude to a signed decimal fraction, if needed."
114 (if (numberp calendar-longitude)
116 (let ((long (+ (aref calendar-longitude 0)
117 (/ (aref calendar-longitude 1) 60.0))))
118 (if (equal (aref calendar-longitude 2) 'east)
123 (defvar calendar-location-name
124 '(let ((float-output-format "%.1f"))
126 (if (numberp calendar-latitude)
127 (abs calendar-latitude)
128 (+ (aref calendar-latitude 0)
129 (/ (aref calendar-latitude 1) 60.0)))
130 (if (numberp calendar-latitude)
131 (if (> calendar-latitude 0) "N" "S")
132 (if (equal (aref calendar-latitude 2) 'north) "N" "S"))
133 (if (numberp calendar-longitude)
134 (abs calendar-longitude)
135 (+ (aref calendar-longitude 0)
136 (/ (aref calendar-longitude 1) 60.0)))
137 (if (numberp calendar-longitude)
138 (if (> calendar-longitude 0) "E" "W")
139 (if (equal (aref calendar-longitude 2) 'east) "E" "W"))))
140 "*Expression evaluating to name of `calendar-longitude', calendar-latitude'.
141 For example, \"New York City\". Default value is just the latitude, longitude
144 This variable should be set in site-local.el.")
146 (defvar solar-n-hemi-seasons
147 '("Vernal Equinox" "Summer Solstice" "Autumnal Equinox" "Winter Solstice")
148 "List of season changes for the northern hemisphere.")
150 (defvar solar-s-hemi-seasons
151 '("Autumnal Equinox" "Winter Solstice" "Vernal Equinox" "Summer Solstice")
152 "List of season changes for the southern hemisphere.")
154 (defun solar-setup ()
155 "Prompt user for latitude, longitude, and time zone."
157 (if (not calendar-longitude)
158 (setq calendar-longitude
160 "Enter longitude (decimal fraction; + east, - west): ")))
161 (if (not calendar-latitude)
162 (setq calendar-latitude
164 "Enter latitude (decimal fraction; + north, - south): ")))
165 (if (not calendar-time-zone)
166 (setq calendar-time-zone
168 "Enter difference from Coordinated Universal Time (in minutes): "))))
170 (defun solar-get-number (prompt)
171 "Return a number from the minibuffer, prompting with PROMPT.
172 Returns nil if nothing was entered."
173 (let ((x (read-string prompt "")))
174 (if (not (string-equal x ""))
177 (defsubst solar-sin-degrees (x)
178 (sin (degrees-to-radians (mod x 360.0))))
180 (defsubst solar-cosine-degrees (x)
181 (cos (degrees-to-radians (mod x 360.0))))
183 (defsubst solar-tangent-degrees (x)
184 (tan (degrees-to-radians (mod x 360.0))))
186 (defun solar-xy-to-quadrant (x y)
187 "Determines the quadrant of the point X, Y."
192 (defun solar-degrees-to-quadrant (angle)
193 "Determines the quadrant of ANGLE."
194 (1+ (floor (mod angle 360) 90)))
196 (defun solar-arctan (x quad)
197 "Arctangent of X in quadrant QUAD."
198 (let ((deg (radians-to-degrees (atan x))))
199 (cond ((equal quad 2) (+ deg 180))
200 ((equal quad 3) (+ deg 180))
201 ((equal quad 4) (+ deg 360))
204 (defun solar-arccos (x)
205 (let ((y (sqrt (- 1 (* x x)))))
206 (solar-arctan (/ y x) (solar-xy-to-quadrant x y))))
208 (defun solar-arcsin (y)
209 (let ((x (sqrt (- 1 (* y y)))))
210 (solar-arctan (/ y x) (solar-xy-to-quadrant x y))))
212 (defconst solar-earth-inclination 23.441884
213 "Inclination of earth's equator to its solar orbit in degrees.")
215 (defconst solar-cos-inclination (solar-cosine-degrees solar-earth-inclination)
216 "Cosine of earth's inclination.")
218 (defconst solar-sin-inclination (solar-sin-degrees solar-earth-inclination)
219 "Sine of earth's inclination.")
221 (defconst solar-earth-orbit-eccentricity 0.016718
222 "Eccentricity of orbit of the earth around the sun.")
224 (defsubst solar-degrees-to-hours (deg)
227 (defsubst solar-hours-to-days (hour)
230 (defun solar-longitude-of-sun (day)
231 "Longitude of the sun at DAY in the year."
232 (let ((mean-anomaly (- (* 0.9856 day) 3.289)))
234 (* 1.916 (solar-sin-degrees mean-anomaly))
235 (* 0.020 (solar-sin-degrees (* 2 mean-anomaly)))
239 (defun solar-right-ascension (longitude)
240 "Right ascension of the sun, given its LONGITUDE."
241 (solar-degrees-to-hours
243 (* solar-cos-inclination (solar-tangent-degrees longitude))
244 (solar-degrees-to-quadrant longitude))))
246 (defun solar-declination (longitude)
247 "Declination of the sun, given its LONGITUDE."
249 (* solar-sin-inclination
250 (solar-sin-degrees longitude))))
252 (defun solar-sunrise (date)
253 "Calculates the *standard* time of sunrise for Gregorian DATE.
254 Calculation is for location given by `calendar-latitude' and
255 `calendar-longitude'.
257 Returns a decimal fraction of hours. Returns nil if the sun does not rise at
258 that location on that day."
259 (let* ((day-of-year (calendar-day-number date))
263 (- 6 (solar-degrees-to-hours (calendar-longitude))))))
264 (solar-longitude-of-sun-at-sunrise
265 (solar-longitude-of-sun approx-sunrise))
266 (solar-right-ascension-at-sunrise
267 (solar-right-ascension solar-longitude-of-sun-at-sunrise))
268 (solar-declination-at-sunrise
269 (solar-declination solar-longitude-of-sun-at-sunrise))
271 (/ (- (solar-cosine-degrees (+ 90 (/ 50.0 60.0)))
272 (* (solar-sin-degrees solar-declination-at-sunrise)
273 (solar-sin-degrees (calendar-latitude))))
274 (* (solar-cosine-degrees solar-declination-at-sunrise)
275 (solar-cosine-degrees (calendar-latitude))))))
276 (if (<= (abs cos-local-sunrise) 1);; otherwise, no sunrise that day
277 (let* ((local-sunrise (solar-degrees-to-hours
278 (- 360 (solar-arccos cos-local-sunrise))))
280 (mod (- (+ local-sunrise solar-right-ascension-at-sunrise)
281 (+ (* 0.065710 approx-sunrise)
284 (+ (- local-mean-sunrise (solar-degrees-to-hours (calendar-longitude)))
285 (/ calendar-time-zone 60.0))))))
287 (defun solar-sunset (date)
288 "Calculates the *standard* time of sunset for Gregorian DATE.
289 Calculation is for location given by `calendar-latitude' and
290 `calendar-longitude'.
292 Returns a decimal fractions of hours. Returns nil if the sun does not set at
293 that location on that day."
294 (let* ((day-of-year (calendar-day-number date))
298 (- 18 (solar-degrees-to-hours (calendar-longitude))))))
299 (solar-longitude-of-sun-at-sunset
300 (solar-longitude-of-sun approx-sunset))
301 (solar-right-ascension-at-sunset
302 (solar-right-ascension solar-longitude-of-sun-at-sunset))
303 (solar-declination-at-sunset
304 (solar-declination solar-longitude-of-sun-at-sunset))
306 (/ (- (solar-cosine-degrees (+ 90 (/ 50.0 60.0)))
307 (* (solar-sin-degrees solar-declination-at-sunset)
308 (solar-sin-degrees (calendar-latitude))))
309 (* (solar-cosine-degrees solar-declination-at-sunset)
310 (solar-cosine-degrees (calendar-latitude))))))
311 (if (<= (abs cos-local-sunset) 1);; otherwise, no sunset that day
312 (let* ((local-sunset (solar-degrees-to-hours
313 (solar-arccos cos-local-sunset)))
315 (mod (- (+ local-sunset solar-right-ascension-at-sunset)
316 (+ (* 0.065710 approx-sunset) 6.622))
318 (+ (- local-mean-sunset (solar-degrees-to-hours (calendar-longitude)))
319 (/ calendar-time-zone 60.0))))))
321 (defun solar-time-string (time time-zone)
322 "Printable form for decimal fraction TIME in TIME-ZONE.
323 Format used is given by `calendar-time-display-form'."
324 (let* ((time (round (* 60 time)))
325 (24-hours (/ time 60))
326 (minutes (format "%02d" (% time 60)))
327 (12-hours (format "%d" (1+ (% (+ 24-hours 11) 12))))
328 (am-pm (if (>= 24-hours 12) "pm" "am"))
329 (24-hours (format "%02d" 24-hours)))
330 (mapconcat 'eval calendar-time-display-form "")))
332 (defun solar-sunrise-sunset (date)
333 "String giving local times of sunrise and sunset on Gregorian DATE."
334 (let* ((rise (solar-sunrise date))
335 (adj-rise (if rise (dst-adjust-time date rise)))
336 (set (solar-sunset date))
337 (adj-set (if set (dst-adjust-time date set))))
338 (format "%s, %s at %s"
339 (if (and rise (calendar-date-equal date (car adj-rise)))
340 (concat "Sunrise " (apply 'solar-time-string (cdr adj-rise)))
342 (if (and set (calendar-date-equal date (car adj-set)))
343 (concat "sunset " (apply 'solar-time-string (cdr adj-set)))
345 (eval calendar-location-name))))
347 (defun solar-date-next-longitude (d l)
348 "First moment on or after Julian day number D when sun's longitude is a
349 multiple of L degrees at calendar-location-name with that location's
350 local time (including any daylight savings rules).
352 L must be an integer divisor of 360.
354 Result is in local time expressed astronomical (Julian) day numbers.
356 The values of calendar-daylight-savings-starts,
357 calendar-daylight-savings-starts-time, calendar-daylight-savings-ends,
358 calendar-daylight-savings-ends-time, calendar-daylight-time-offset, and
359 calendar-time-zone are used to interpret local time."
362 (start-long (solar-longitude d))
363 (next (mod (* l (1+ (floor (/ start-long l)))) 360))
364 (end (+ d (* (/ l 360.0) 400)))
365 (end-long (solar-longitude end)))
366 (while ;; bisection search for nearest minute
367 (< 0.00001 (- end start))
369 ;; start-long <= next < end-long when next != 0
370 ;; when next = 0, we look for the discontinuity (start-long is near 360
371 ;; and end-long is small (less than l).
372 (setq d (/ (+ start end) 2.0))
373 (setq long (solar-longitude d))
374 (if (or (and (/= next 0) (< long next))
375 (and (= next 0) (< l long)))
378 (setq start-long long))
380 (setq end-long long)))
381 (/ (+ start end) 2.0)))
383 (defun solar-longitude (d)
384 "Longitude of sun on astronomical (Julian) day number D.
385 Accurary is about 0.01 degree (about 365.25*24*60*0.01/360 = 15 minutes).
387 The values of calendar-daylight-savings-starts,
388 calendar-daylight-savings-starts-time, calendar-daylight-savings-ends,
389 calendar-daylight-savings-ends-time, calendar-daylight-time-offset, and
390 calendar-time-zone are used to interpret local time."
391 (let* ((a-d (calendar-absolute-from-astro d))
392 (date (calendar-gregorian-from-absolute (floor a-d)))
393 (time (* 24 (- a-d (truncate a-d))))
394 (rounded-abs-date (+ (calendar-absolute-from-gregorian date)
395 (/ (round (* 60 time)) 60.0 24.0)))
396 ;; get local standard time
397 (a-d (- rounded-abs-date
398 (if (dst-in-effect rounded-abs-date)
399 (/ calendar-daylight-time-offset 24.0 60.0) 0)))
400 ;; get Universal Time
401 (a-d (- a-d (/ calendar-time-zone 60.0 24.0)))
402 (date (calendar-astro-from-absolute a-d))
403 ;; get Ephemeris Time
404 (date (+ date (solar-ephemeris-correction
405 (extract-calendar-year
406 (calendar-gregorian-from-absolute
408 (calendar-absolute-from-astro
410 (T (/ (- date 2451545.0) 36525.0))
411 (Lo (mod (+ 280.46645 (* 36000.76983 T) (* 0.0003032 T T)) 360.0))
415 (* -0.00000048 T T T))
417 (e (+ 0.016708617 (* -0.000042037 T) (* -0.0000001236 T T)))
418 (C (+ (* (+ 1.914600 (* -0.004817 T) (* -0.000014 T T))
419 (solar-sin-degrees M))
420 (* (+ 0.019993 (* -0.000101 T)) (solar-sin-degrees (* 2 M)))
421 (* 0.000290 (solar-sin-degrees (* 3 M)))))
422 (true-longitude (+ Lo C))
423 (omega (+ 125.04 (* -1934.136 T)))
424 (apparent-longitude (mod
427 (* -0.00478 (solar-sin-degrees omega)))
431 (defun solar-ephemeris-correction (year)
432 "Ephemeris time minus Universal Time at astronomical (Julian) day D.
433 Result is in days For the years 1800-1987, the maximum error is 1.9 seconds.
434 For the other years, the maximum error is about 30 seconds."
435 (cond ((and (<= 1988 year) (< year 2020))
436 (/ (+ year -2000 67.0) 60.0 60.0 24.0))
437 ((and (<= 1900 year) (< year 1988))
438 (let* ((theta (/ (- (calendar-astro-from-absolute
439 (calendar-absolute-from-gregorian
441 (calendar-astro-from-absolute
442 (calendar-absolute-from-gregorian
445 (theta2 (* theta theta))
446 (theta3 (* theta2 theta))
447 (theta4 (* theta2 theta2))
448 (theta5 (* theta3 theta2)))
455 (* 0.677066 theta3 theta3)
456 (* -0.212591 theta4 theta3))))
457 ((and (<= 1800 year) (< year 1900))
458 (let* ((theta (/ (- (calendar-astro-from-absolute
459 (calendar-absolute-from-gregorian
461 (calendar-astro-from-absolute
462 (calendar-absolute-from-gregorian
465 (theta2 (* theta theta))
466 (theta3 (* theta2 theta))
467 (theta4 (* theta2 theta2))
468 (theta5 (* theta3 theta2)))
475 (* 31.332267 theta3 theta3)
476 (* 38.291999 theta4 theta3)
477 (* 28.316289 theta4 theta4)
478 (* 11.636204 theta4 theta5)
479 (* 2.043794 theta5 theta5))))
480 ((and (<= 1620 year) (< year 1800))
481 (let ((x (/ (- year 1600) 10.0)))
482 (/ (+ (* 2.19167 x x) (* -40.675 x) 196.58333) 60.0 60.0 24.0)))
483 (t (let* ((tmp (- (calendar-astro-from-absolute
484 (calendar-absolute-from-gregorian
487 (second (- (/ (* tmp tmp) 41048480.0) 15)))
488 (/ second 60.0 60.0 24.0)))))
491 (defun sunrise-sunset (&optional arg)
492 "Local time of sunrise and sunset for today. Accurate to +/- 2 minutes.
493 If called with an optional prefix argument, prompt for date.
495 If called with an optional double prefix argument, prompt for longitude,
496 latitude, time zone, and date, and always use standard time.
498 This function is suitable for execution in a .emacs file."
500 (or arg (setq arg 1))
502 (not (and calendar-latitude calendar-longitude calendar-time-zone)))
504 (let* ((calendar-longitude
505 (if (< arg 16) calendar-longitude
507 "Enter longitude (decimal fraction; + east, - west): ")))
509 (if (< arg 16) calendar-latitude
511 "Enter latitude (decimal fraction; + north, - south): ")))
513 (if (< arg 16) calendar-time-zone
515 "Enter difference from Coordinated Universal Time (in minutes): ")))
516 (calendar-location-name
517 (if (< arg 16) calendar-location-name
518 (let ((float-output-format "%.1f"))
520 (if (numberp calendar-latitude)
521 (abs calendar-latitude)
522 (+ (aref calendar-latitude 0)
523 (/ (aref calendar-latitude 1) 60.0)))
524 (if (numberp calendar-latitude)
525 (if (> calendar-latitude 0) "N" "S")
526 (if (equal (aref calendar-latitude 2) 'north) "N" "S"))
527 (if (numberp calendar-longitude)
528 (abs calendar-longitude)
529 (+ (aref calendar-longitude 0)
530 (/ (aref calendar-longitude 1) 60.0)))
531 (if (numberp calendar-longitude)
532 (if (> calendar-longitude 0) "E" "W")
533 (if (equal (aref calendar-longitude 2) 'east)
535 (calendar-standard-time-zone-name
536 (if (< arg 16) calendar-standard-time-zone-name
537 (cond ((= calendar-time-zone 0) "UTC")
538 ((< calendar-time-zone 0)
539 (format "UTC%dmin" calendar-time-zone))
540 (t (format "UTC+%dmin" calendar-time-zone)))))
541 (calendar-daylight-savings-starts
542 (if (< arg 16) calendar-daylight-savings-starts))
543 (calendar-daylight-savings-ends
544 (if (< arg 16) calendar-daylight-savings-ends))
545 (date (if (< arg 4) (calendar-current-date) (calendar-read-date)))
546 (date-string (calendar-date-string date t))
547 (time-string (solar-sunrise-sunset date))
548 (msg (format "%s: %s" date-string time-string))
549 (one-window (one-window-p t)))
550 (if (<= (length msg) (frame-width))
552 (with-output-to-temp-buffer "*temp*"
553 (princ (concat date-string "\n" time-string)))
554 (message (substitute-command-keys
557 "Type \\[delete-other-windows] to remove temp window."
558 "Type \\[switch-to-buffer] RET to remove temp window.")
559 "Type \\[switch-to-buffer-other-window] RET to restore old contents of temp window."))))))
561 (defun calendar-sunrise-sunset ()
562 "Local time of sunrise and sunset for date under cursor.
563 Accurate to +/- 2 minutes."
565 (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
567 (let ((date (calendar-cursor-to-date t)))
569 (calendar-date-string date t t)
570 (solar-sunrise-sunset date))))
572 (defun diary-sunrise-sunset ()
573 "Local time of sunrise and sunset as a diary entry.
574 Accurate to +/- 2 minutes."
575 (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
577 (solar-sunrise-sunset date))
579 (defun diary-sabbath-candles ()
580 "Local time of candle lighting diary entry--applies if date is a Friday.
581 No diary entry if there is no sunset on that date."
582 (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
584 (if (= (% (calendar-absolute-from-gregorian date) 7) 5);; Friday
585 (let* ((sunset (solar-sunset date))
589 (- sunset (/ 18.0 60.0))))))
590 (if (and light (calendar-date-equal date (car light)))
591 (format "%s Sabbath candle lighting"
592 (apply 'solar-time-string (cdr light)))))))
595 (defun solar-equinoxes-solstices ()
596 "Date and time of equinoxes and solstices, if visible in the calendar window.
597 Requires floating point."
598 (let ((m displayed-month)
600 (increment-calendar-month m y (cond ((= 1 (% m 3)) -1)
603 (let* ((calendar-standard-time-zone-name
604 (if calendar-time-zone calendar-standard-time-zone-name "UTC"))
605 (calendar-daylight-savings-starts
606 (if calendar-time-zone calendar-daylight-savings-starts))
607 (calendar-daylight-savings-ends
608 (if calendar-time-zone calendar-daylight-savings-ends))
609 (calendar-time-zone (if calendar-time-zone calendar-time-zone 0))
611 (d (solar-date-next-longitude
612 (calendar-astro-from-absolute
613 (calendar-absolute-from-gregorian
614 (list (+ 3 (* k 3)) 15 y)))
616 (abs-day (calendar-absolute-from-astro d)))
618 (list (calendar-gregorian-from-absolute (floor abs-day))
620 (nth k (if (and calendar-latitude
621 (< (calendar-latitude) 0))
623 solar-n-hemi-seasons))
625 (* 24 (- abs-day (floor abs-day)))
626 (if (dst-in-effect abs-day)
627 calendar-daylight-time-zone-name
628 calendar-standard-time-zone-name))))))))
633 ;;; solar.el ends here