code.delx.au - gnu-emacs/blob - lisp/calendar/solar.el

   1 ;;; solar.el --- calendar functions for solar events.
   2
   3 ;; Copyright (C) 1992 Free Software Foundation, Inc.
   4
   5 ;; Author: Edward M. Reingold <reingold@cs.uiuc.edu>
   6 ;; Keywords: calendar
   7 ;; Human-Keywords: sunrise, sunset, equinox, solstice, calendar, diary,
   8 ;;      holidays
   9
  10 ;; This file is part of GNU Emacs.
  11
  12 ;; GNU Emacs is distributed in the hope that it will be useful,
  13 ;; but WITHOUT ANY WARRANTY.  No author or distributor
  14 ;; accepts responsibility to anyone for the consequences of using it
  15 ;; or for whether it serves any particular purpose or works at all,
  16 ;; unless he says so in writing.  Refer to the GNU Emacs General Public
  17 ;; License for full details.
  18
  19 ;; Everyone is granted permission to copy, modify and redistribute
  20 ;; GNU Emacs, but only under the conditions described in the
  21 ;; GNU Emacs General Public License.   A copy of this license is
  22 ;; supposed to have been given to you along with GNU Emacs so you
  23 ;; can know your rights and responsibilities.  It should be in a
  24 ;; file named COPYING.  Among other things, the copyright notice
  25 ;; and this notice must be preserved on all copies.
  26
  27 ;;; Commentary:
  28
  29 ;; This collection of functions implements the features of calendar.el and
  30 ;; diary.el that deal with sunrise/sunset and eqinoxes/solstices.
  31
  32 ;; Based on the ``Almanac for Computers 1984,'' prepared by the Nautical
  33 ;; Almanac Office, United States Naval Observatory, Washington, 1984 and
  34 ;; on ``Astronomical Formulae for Calculators,'' 3rd ed., by Jean Meeus,
  35 ;; Willmann-Bell, Inc., 1985.
  36 ;;
  37 ;; WARNINGS:
  38 ;;    1. SUNRISE/SUNSET calculations will be accurate only to +/- 2 minutes.
  39 ;;       Locations should be between +/- 65 degrees of latitude.
  40 ;;       Dates should be in the latter half of the 20th century.
  41 ;;
  42 ;;    2. Equinox/solstice times will be accurate only to +/- 15 minutes.
  43
  44 ;; The author would be delighted to have an astronomically more sophisticated
  45 ;; person rewrite the code for the solar calculations in this file!
  46
  47 ;; Comments, corrections, and improvements should be sent to
  48 ;;  Edward M. Reingold               Department of Computer Science
  49 ;;  (217) 333-6733                   University of Illinois at Urbana-Champaign
  50 ;;  reingold@cs.uiuc.edu             1304 West Springfield Avenue
  51 ;;                                   Urbana, Illinois 61801
  52
  53 ;;; Code:
  54
  55 (if (fboundp 'atan)
  56     (require 'lisp-float-type)
  57   (error "Solar calculations impossible since floating point is unavailable."))
  58
  59 (require 'calendar)
  60
  61 (defun solar-setup ()
  62   "Prompt user for latitude, longitude, and time zone."
  63   (beep)
  64   (if (not calendar-longitude)
  65       (setq calendar-longitude
  66             (solar-get-number
  67              "Enter longitude (decimal fraction; + east, - west): ")))
  68   (if (not calendar-latitude)
  69       (setq calendar-latitude
  70             (solar-get-number
  71              "Enter latitude (decimal fraction; + north, - south): ")))
  72   (if (not calendar-time-zone)
  73       (setq calendar-time-zone
  74             (solar-get-number
  75              "Enter difference from Universal Time (in minutes): "))))
  76
  77 (defun solar-get-number (prompt)
  78   "Return a number from the minibuffer, prompting with PROMPT.
  79 Returns nil if nothing was entered."
  80   (let ((x (read-string prompt "")))
  81     (if (not (string-equal x ""))
  82         (string-to-int x))))
  83
  84 (defun solar-sin-degrees (x)
  85   (sin (degrees-to-radians x)))
  86
  87 (defun solar-cosine-degrees (x)
  88   (cos (degrees-to-radians x)))
  89
  90 (defun solar-tangent-degrees (x)
  91   (tan (degrees-to-radians x)))
  92
  93 (defun solar-xy-to-quadrant (x y)
  94   "Determines the quadrant of the point X, Y."
  95   (if (> x 0)
  96       (if (> y 0) 1 4)
  97       (if (> y 0) 2 3)))
  98
  99 (defun solar-degrees-to-quadrant (angle)
 100   "Determines the quadrant of ANGLE."
 101   (1+ (truncate (/ (solar-mod angle 360.0) 90.0))))
 102
 103 (defun solar-arctan (x quad)
 104   "Arctangent of X in quadrant QUAD."
 105   (let ((deg (radians-to-degrees (atan x))))
 106     (cond ((equal quad 2)   (+ deg 180))
 107           ((equal quad 3)   (+ deg 180))
 108           ((equal quad 4)   (+ deg 360))
 109           (t                deg))))
 110
 111 (defun solar-arccos (x)
 112   (let ((y (sqrt (- 1 (* x x)))))
 113     (solar-arctan (/ y x) (solar-xy-to-quadrant x y))))
 114
 115 (defun solar-arcsin (y)
 116   (let ((x (sqrt (- 1 (* y y)))))
 117     (solar-arctan (/ y x) (solar-xy-to-quadrant x y))))
 118
 119 (defun solar-mod (x y)
 120   "Returns X mod Y; value is *always* non-negative."
 121   (let ((v (% x y)))
 122     (if (> 0 v)
 123         (+ v y)
 124       v)))
 125
 126 (defconst solar-earth-inclination 23.441884
 127   "Inclination of earth's equator to its solar orbit in degrees.")
 128
 129 (defconst solar-cos-inclination (solar-cosine-degrees solar-earth-inclination)
 130   "Cosine of earth's inclination.")
 131
 132 (defconst solar-sin-inclination (solar-sin-degrees solar-earth-inclination)
 133   "Sine of earth's inclination.")
 134
 135 (defconst solar-earth-orbit-eccentricity 0.016718
 136   "Eccentricity of orbit of the earth around the sun.")
 137
 138 (defmacro solar-degrees-to-hours (deg)
 139   (list '/ deg 15))
 140
 141 (defmacro solar-hours-to-days (hour)
 142   (list '/ hour 24))
 143
 144 (defun solar-longitude-of-sun (day)
 145   "Longitude of the sun at DAY in the year."
 146   (let ((mean-anomaly (- (* 0.9856 day) 3.289)))
 147     (solar-mod (+ mean-anomaly
 148                   (* 1.916 (solar-sin-degrees mean-anomaly))
 149                   (* 0.020 (solar-sin-degrees (* 2 mean-anomaly)))
 150                   282.634)
 151                360)))
 152
 153 (defun solar-right-ascension (longitude)
 154   "Right ascension of the sun, given its LONGITUDE."
 155   (solar-degrees-to-hours
 156    (solar-arctan
 157     (* solar-cos-inclination (solar-tangent-degrees longitude))
 158     (solar-degrees-to-quadrant longitude))))
 159
 160 (defun solar-declination (longitude)
 161   "Declination of the sun, given its LONGITUDE."
 162   (solar-arcsin
 163    (* solar-sin-inclination
 164       (solar-sin-degrees longitude))))
 165
 166 (defun solar-sunrise (date)
 167   "Calculates the *standard* time of sunrise for Gregorian DATE for location
 168 given by `calendar-latitude' and `calendar-longitude'.  Returns a decimal fraction
 169 of hours.  Returns nil if the sun does not rise at that location on that day."
 170   (let* ((day-of-year (calendar-day-number date))
 171          (approx-sunrise
 172           (+ day-of-year
 173              (solar-hours-to-days
 174               (-  6 (solar-degrees-to-hours calendar-longitude)))))
 175          (solar-longitude-of-sun-at-sunrise
 176           (solar-longitude-of-sun approx-sunrise))
 177          (solar-right-ascension-at-sunrise
 178           (solar-right-ascension solar-longitude-of-sun-at-sunrise))
 179          (solar-declination-at-sunrise
 180           (solar-declination solar-longitude-of-sun-at-sunrise))
 181          (cos-local-sunrise
 182           (/ (- (solar-cosine-degrees (+ 90 (/ 50.0 60.0)))
 183                 (* (solar-sin-degrees solar-declination-at-sunrise)
 184                    (solar-sin-degrees calendar-latitude)))
 185              (* (solar-cosine-degrees solar-declination-at-sunrise)
 186                 (solar-cosine-degrees calendar-latitude)))))
 187     (if (<= (abs cos-local-sunrise) 1);; otherwise, no sunrise that day
 188       (let* ((local-sunrise (solar-degrees-to-hours
 189                              (- 360 (solar-arccos cos-local-sunrise))))
 190              (local-mean-sunrise
 191               (solar-mod (- (+ local-sunrise solar-right-ascension-at-sunrise)
 192                             (+ (* 0.065710 approx-sunrise)
 193                                6.622))
 194                          24)))
 195         (+ (- local-mean-sunrise (solar-degrees-to-hours calendar-longitude))
 196            (/ calendar-time-zone 60.0))))))
 197
 198 (defun solar-sunset (date)
 199   "Calculates the *standard* time of sunset for Gregorian DATE for location
 200 given by `calendar-latitude' and `calendar-longitude'.  Returns a decimal fractions
 201 of hours.  Returns nil if the sun does not set at that location on that day."
 202   (let* ((day-of-year (calendar-day-number date))
 203          (approx-sunset
 204           (+ day-of-year
 205              (solar-hours-to-days
 206               (- 18 (solar-degrees-to-hours calendar-longitude)))))
 207          (solar-longitude-of-sun-at-sunset
 208           (solar-longitude-of-sun approx-sunset))
 209          (solar-right-ascension-at-sunset
 210           (solar-right-ascension solar-longitude-of-sun-at-sunset))
 211          (solar-declination-at-sunset
 212           (solar-declination solar-longitude-of-sun-at-sunset))
 213          (cos-local-sunset
 214           (/ (- (solar-cosine-degrees (+ 90 (/ 50.0 60.0)))
 215                 (* (solar-sin-degrees solar-declination-at-sunset)
 216                    (solar-sin-degrees calendar-latitude)))
 217              (* (solar-cosine-degrees solar-declination-at-sunset)
 218                 (solar-cosine-degrees calendar-latitude)))))
 219     (if (<= (abs cos-local-sunset) 1);; otherwise, no sunset that day
 220       (let* ((local-sunset (solar-degrees-to-hours
 221                             (solar-arccos cos-local-sunset)))
 222              (local-mean-sunset
 223               (solar-mod (- (+ local-sunset solar-right-ascension-at-sunset)
 224                             (+ (* 0.065710 approx-sunset) 6.622))
 225                          24)))
 226         (+ (- local-mean-sunset (solar-degrees-to-hours calendar-longitude))
 227            (/ calendar-time-zone 60.0))))))
 228
 229 (defun solar-time-string (time date)
 230   "Printable form for decimal fraction standard TIME on DATE.
 231 Format used is given by `calendar-time-display-form'.  Converted to daylight
 232 savings time according to `calendar-daylight-savings-starts' and
 233 `calendar-daylight-savings-ends', if those variables are not nil."
 234   (let* ((year (extract-calendar-year date))
 235          (abs-date (calendar-absolute-from-gregorian date))
 236          (dst (and calendar-daylight-savings-starts
 237                    calendar-daylight-savings-ends
 238                    (<= (calendar-absolute-from-gregorian
 239                         (eval calendar-daylight-savings-starts))
 240                        abs-date)
 241                    (< abs-date
 242                       (calendar-absolute-from-gregorian
 243                        (eval calendar-daylight-savings-ends)))))
 244          (time (if dst (1+ time) time))
 245          (time-zone (if dst
 246                         calendar-daylight-time-zone-name
 247                         calendar-standard-time-zone-name))
 248          (24-hours (truncate time))
 249          (minutes (round (* 60 (- time 24-hours))))
 250          (24-hours (if (= minutes 60) (1+ 24-hours) 24-hours))
 251          (minutes (if (= minutes 60) 0 minutes))
 252          (minutes (format "%02d" minutes))
 253          (12-hours (format "%d" (if (> 24-hours 12)
 254                                     (- 24-hours 12)
 255                                   (if (= 24-hours 0) 12 24-hours))))
 256          (am-pm (if (>= 24-hours 12) "pm" "am"))
 257          (24-hours (format "%02d" 24-hours)))
 258     (mapconcat 'eval calendar-time-display-form "")))
 259
 260 (defun solar-sunrise-sunset (date)
 261   "String giving local times of sunrise and sunset on Gregorian DATE."
 262   (let ((rise (solar-sunrise date))
 263         (set (solar-sunset date)))
 264     (format "%s, %s at %s"
 265             (if rise
 266                 (concat "Sunrise " (solar-time-string rise date))
 267               "No sunrise")
 268             (if set
 269                 (concat "sunset " (solar-time-string set date))
 270               "no sunset")
 271             (eval calendar-location-name))))
 272
 273 (defun solar-apparent-longitude-of-sun (date)
 274   "Apparent longitude of the sun on Gregorian DATE."
 275   (let* ((time (/ (- (calendar-absolute-from-gregorian date)
 276                      (calendar-absolute-from-gregorian '(1 0.5 1900)))
 277                   36525))
 278          (l (+ 279.69668
 279                (* 36000.76892 time)
 280                (* 0.0003025 time time)))
 281          (m (+ 358.47583
 282                (* 35999.04975 time)
 283                (* -0.000150 time time)
 284                (* -0.0000033 time time time)))
 285          (c (+ (* (+ 1.919460
 286                      (* -0.004789 time)
 287                      (* -0.000014 time time))
 288                   (solar-sin-degrees m))
 289                (* (+ 0.020094
 290                      (* -0.000100 time))
 291                   (solar-sin-degrees (* 2 m)))
 292                (* 0.000293
 293                   (solar-sin-degrees (* 3 m)))))
 294          (L (+ l c))
 295          (omega (+ 259.18
 296                    (* -1934.142 time)))
 297          (app (+ L
 298                  -0.00569
 299                  (* -0.00479
 300                     (solar-sin-degrees omega)))))
 301     app))
 302
 303 (defun solar-ephemeris-correction (year)
 304   "Difference in minutes between Ephemeris time an Universal time in YEAR.
 305 Value is only an approximation."
 306   (let ((T (/ (- year 1900) 100.0)))
 307     (+ 0.41 (* 1.2053 T) (* 0.4992 T T))))
 308
 309 (defun solar-equinoxes/solstices (k year)
 310   "Date of equinox/solstice K for YEAR.  K=0, spring equinox; K=1, summer
 311 solstice; K=2, fall equinox; K=3, winter solstice.  Accurate to within
 312 several minutes."
 313   (let ((date (list (+ 3 (* k 3)) 21 year))
 314         (correction 1000))
 315     (while (> correction 0.00001)
 316       (setq app (solar-mod (solar-apparent-longitude-of-sun date) 360.0))
 317       (setq correction (* 58 (solar-sin-degrees (- (* k 90) app))))
 318       (setq date (list (extract-calendar-month date)
 319                        (+ (extract-calendar-day date) correction)
 320                        year)))
 321     (list (extract-calendar-month date)
 322           (+ (extract-calendar-day date) (/ calendar-time-zone 60.0 24.0)
 323              (- (/ (solar-ephemeris-correction year) 60.0 24.0)))
 324           year)))
 325
 326 ;;;###autoload
 327 (defun sunrise-sunset (&optional arg)
 328   "Local time of sunrise and sunset for today.  Accurate to +/- 2 minutes.
 329 If called with an optional prefix argument, prompts for date.
 330
 331 If called with an optional double prefix argument, prompts for longitude,
 332 latitude, time zone, and date.
 333
 334 This function is suitable for execution in a .emacs file."
 335  (interactive "p")
 336  (if (and (< arg 16)
 337           (not (and calendar-latitude calendar-longitude calendar-time-zone)))
 338      (solar-setup))
 339  (let* ((calendar-longitude
 340          (if (< arg 16) calendar-longitude
 341            (solar-get-number
 342             "Enter longitude (decimal fraction; + east, - west): ")))
 343         (calendar-latitude
 344          (if (< arg 16) calendar-latitude
 345            (solar-get-number
 346             "Enter latitude (decimal fraction; + north, - south): ")))
 347         (calendar-time-zone
 348          (if (< arg 16) calendar-time-zone
 349            (solar-get-number
 350             "Enter difference from Universal Time (in minutes): ")))
 351         (calendar-location-name
 352          (if (< arg 16) calendar-location-name
 353            (let ((float-output-format "%.1f"))
 354              (format "%s%s, %s%s"
 355                      (abs calendar-latitude)
 356                      (if (> calendar-latitude 0) "N" "S")
 357                      (abs calendar-longitude)
 358                      (if (> calendar-longitude 0) "E" "W")))))
 359         (calendar-standard-time-zone-name
 360          (if (< arg 16) calendar-standard-time-zone-name
 361            (cond ((= calendar-time-zone 0) "UT")
 362                  ((< calendar-time-zone 0)
 363                      (format "UT%dmin" calendar-time-zone))
 364                  (t  (format "UT+%dmin" calendar-time-zone)))))
 365         (calendar-daylight-savings-starts
 366          (if (< arg 16) calendar-daylight-savings-starts))
 367         (calendar-daylight-savings-ends
 368          (if (< arg 16) calendar-daylight-savings-ends))
 369         (date (if (< arg 4) (calendar-current-date) (calendar-read-date)))
 370         (date-string (calendar-date-string date t))
 371         (time-string (solar-sunrise-sunset date))
 372         (msg (format "%s: %s" date-string time-string))
 373         (one-window (one-window-p t)))
 374    (if (<= (length msg) (screen-width))
 375        (message msg)
 376      (with-output-to-temp-buffer "*temp*"
 377        (princ (concat date-string "\n" time-string)))
 378      (message (substitute-command-keys
 379                (if one-window
 380                    (if pop-up-windows
 381                        "Type \\[delete-other-windows] to remove temp window."
 382                      "Type \\[switch-to-buffer] RET to remove temp window.")
 383                  "Type \\[switch-to-buffer-other-window] RET to restore old contents of temp window."))))))
 384
 385 (defun calendar-sunrise-sunset ()
 386   "Local time of sunrise and sunset for date under cursor.
 387 Accurate to +/- 2 minutes."
 388   (interactive)
 389   (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
 390       (solar-setup))
 391   (message
 392    (solar-sunrise-sunset
 393     (or (calendar-cursor-to-date)
 394         (error "Cursor is not on a date!")))))
 395
 396 (defun diary-sunrise-sunset ()
 397   "Local time of sunrise and sunset as a diary entry.
 398 Accurate to +/- 2 minutes."
 399   (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
 400       (solar-setup))
 401   (solar-sunrise-sunset date))
 402
 403 (defun diary-sabbath-candles ()
 404   "Local time of candle lighting diary entry--applies if date is a Friday.
 405 No diary entry if there is no sunset on that date."
 406   (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
 407       (solar-setup))
 408   (if (= (% (calendar-absolute-from-gregorian date) 7) 5);;  Friday
 409       (let* ((sunset (solar-sunset date))
 410              (light (if sunset (- sunset (/ 18.0 60.0)))))
 411         (if light (format "%s Sabbath candle lighting"
 412                           (solar-time-string light date))))))
 413
 414 (defun calendar-holiday-function-solar-equinoxes-solstices ()
 415   "Date and time of equinoxes and solstices, if visible in the calendar window.
 416 Requires floating point."
 417   (let* ((m displayed-month)
 418          (y displayed-year))
 419     (increment-calendar-month m y (cond ((= 1 (% m 3)) -1)
 420                                         ((= 2 (% m 3))  1)
 421                                         (t              0)))
 422     (let* ((calendar-standard-time-zone-name
 423             (if calendar-time-zone calendar-standard-time-zone-name "UT"))
 424            (calendar-daylight-savings-starts
 425             (if calendar-time-zone calendar-daylight-savings-starts))
 426            (calendar-daylight-savings-ends
 427             (if calendar-time-zone calendar-daylight-savings-ends))
 428            (calendar-time-zone (if calendar-time-zone calendar-time-zone 0))
 429            (k (1- (/ m 3)))
 430            (date (solar-equinoxes/solstices k y))
 431            (day (extract-calendar-day date))
 432            (time (* 24 (- day (truncate day))))
 433            ;; Time zone/DST can't move the date out of range,
 434            ;; so let solar-time-string do the conversion.
 435            (date (list (extract-calendar-month date)
 436                        (truncate day)
 437                        (extract-calendar-year date))))
 438       (list (list date
 439                   (format "%s %s"
 440                           (cond ((= k 0)  "Vernal Equinox")
 441                                 ((= k 1)  "Summer Solstice")
 442                                 ((= k 2)  "Fall Equinox")
 443                                 ((= k 3)  "Winter Solstice"))
 444                           (solar-time-string time date)))))))
 445
 446 (provide 'solar)
 447
 448 ;;; solar.el ends here