1;;; solar.el --- calendar functions for solar events 2 3;; Copyright (C) 1992, 1993, 1995, 1997, 2001, 2002, 2003, 2004, 2005, 4;; 2006, 2007 Free Software Foundation, Inc. 5 6;; Author: Edward M. Reingold <reingold@cs.uiuc.edu> 7;; Denis B. Roegel <Denis.Roegel@loria.fr> 8;; Maintainer: Glenn Morris <rgm@gnu.org> 9;; Keywords: calendar 10;; Human-Keywords: sunrise, sunset, equinox, solstice, calendar, diary, 11;; holidays 12 13;; This file is part of GNU Emacs. 14 15;; GNU Emacs is free software; you can redistribute it and/or modify 16;; it under the terms of the GNU General Public License as published by 17;; the Free Software Foundation; either version 2, or (at your option) 18;; any later version. 19 20;; GNU Emacs is distributed in the hope that it will be useful, 21;; but WITHOUT ANY WARRANTY; without even the implied warranty of 22;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 23;; GNU General Public License for more details. 24 25;; You should have received a copy of the GNU General Public License 26;; along with GNU Emacs; see the file COPYING. If not, write to the 27;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 28;; Boston, MA 02110-1301, USA. 29 30;;; Commentary: 31 32;; This collection of functions implements the features of calendar.el, 33;; diary.el, and holiday.el that deal with times of day, sunrise/sunset, and 34;; equinoxes/solstices. 35 36;; Based on the ``Almanac for Computers 1984,'' prepared by the Nautical 37;; Almanac Office, United States Naval Observatory, Washington, 1984, on 38;; ``Astronomical Formulae for Calculators,'' 3rd ed., by Jean Meeus, 39;; Willmann-Bell, Inc., 1985, on ``Astronomical Algorithms'' by Jean Meeus, 40;; Willmann-Bell, Inc., 1991, and on ``Planetary Programs and Tables from 41;; -4000 to +2800'' by Pierre Bretagnon and Jean-Louis Simon, Willmann-Bell, 42;; Inc., 1986. 43 44;; 45;; Accuracy: 46;; 1. Sunrise/sunset times will be accurate to the minute for years 47;; 1951--2050. For other years the times will be within +/- 2 minutes. 48;; 49;; 2. Equinox/solstice times will be accurate to the minute for years 50;; 1951--2050. For other years the times will be within +/- 1 minute. 51 52;; Technical details of all the calendrical calculations can be found in 53;; ``Calendrical Calculations: The Millennium Edition'' by Edward M. Reingold 54;; and Nachum Dershowitz, Cambridge University Press (2001). 55 56;;; Code: 57 58(defvar date) 59(defvar displayed-month) 60(defvar displayed-year) 61 62(if (fboundp 'atan) 63 (require 'lisp-float-type) 64 (error "Solar/lunar calculations impossible since floating point is unavailable")) 65 66(require 'cal-dst) 67(require 'cal-julian) 68 69;;;###autoload 70(defcustom calendar-time-display-form 71 '(12-hours ":" minutes am-pm 72 (if time-zone " (") time-zone (if time-zone ")")) 73 "*The pseudo-pattern that governs the way a time of day is formatted. 74 75A pseudo-pattern is a list of expressions that can involve the keywords 76`12-hours', `24-hours', and `minutes', all numbers in string form, 77and `am-pm' and `time-zone', both alphabetic strings. 78 79For example, the form 80 81 '(24-hours \":\" minutes 82 (if time-zone \" (\") time-zone (if time-zone \")\")) 83 84would give military-style times like `21:07 (UTC)'." 85 :type 'sexp 86 :group 'calendar) 87 88;;;###autoload 89(defcustom calendar-latitude nil 90 "*Latitude of `calendar-location-name' in degrees. 91 92The value can be either a decimal fraction (one place of accuracy is 93sufficient), + north, - south, such as 40.7 for New York City, or the value 94can be a vector [degrees minutes north/south] such as [40 50 north] for New 95York City. 96 97This variable should be set in `site-start'.el." 98 :type '(choice (const nil) 99 (number :tag "Exact") 100 (vector :value [0 0 north] 101 (integer :tag "Degrees") 102 (integer :tag "Minutes") 103 (choice :tag "Position" 104 (const north) 105 (const south)))) 106 :group 'calendar) 107 108;;;###autoload 109(defcustom calendar-longitude nil 110 "*Longitude of `calendar-location-name' in degrees. 111 112The value can be either a decimal fraction (one place of accuracy is 113sufficient), + east, - west, such as -73.9 for New York City, or the value 114can be a vector [degrees minutes east/west] such as [73 55 west] for New 115York City. 116 117This variable should be set in `site-start'.el." 118 :type '(choice (const nil) 119 (number :tag "Exact") 120 (vector :value [0 0 west] 121 (integer :tag "Degrees") 122 (integer :tag "Minutes") 123 (choice :tag "Position" 124 (const east) 125 (const west)))) 126 :group 'calendar) 127 128(defsubst calendar-latitude () 129 "Convert calendar-latitude to a signed decimal fraction, if needed." 130 (if (numberp calendar-latitude) 131 calendar-latitude 132 (let ((lat (+ (aref calendar-latitude 0) 133 (/ (aref calendar-latitude 1) 60.0)))) 134 (if (equal (aref calendar-latitude 2) 'north) 135 lat 136 (- lat))))) 137 138(defsubst calendar-longitude () 139 "Convert calendar-longitude to a signed decimal fraction, if needed." 140 (if (numberp calendar-longitude) 141 calendar-longitude 142 (let ((long (+ (aref calendar-longitude 0) 143 (/ (aref calendar-longitude 1) 60.0)))) 144 (if (equal (aref calendar-longitude 2) 'east) 145 long 146 (- long))))) 147 148;;;###autoload 149(defcustom calendar-location-name 150 '(let ((float-output-format "%.1f")) 151 (format "%s%s, %s%s" 152 (if (numberp calendar-latitude) 153 (abs calendar-latitude) 154 (+ (aref calendar-latitude 0) 155 (/ (aref calendar-latitude 1) 60.0))) 156 (if (numberp calendar-latitude) 157 (if (> calendar-latitude 0) "N" "S") 158 (if (equal (aref calendar-latitude 2) 'north) "N" "S")) 159 (if (numberp calendar-longitude) 160 (abs calendar-longitude) 161 (+ (aref calendar-longitude 0) 162 (/ (aref calendar-longitude 1) 60.0))) 163 (if (numberp calendar-longitude) 164 (if (> calendar-longitude 0) "E" "W") 165 (if (equal (aref calendar-longitude 2) 'east) "E" "W")))) 166 "*Expression evaluating to name of `calendar-longitude', `calendar-latitude'. 167For example, \"New York City\". Default value is just the latitude, longitude 168pair. 169 170This variable should be set in `site-start'.el." 171 :type 'sexp 172 :group 'calendar) 173 174(defcustom solar-error 0.5 175"*Tolerance (in minutes) for sunrise/sunset calculations. 176 177A larger value makes the calculations for sunrise/sunset faster, but less 178accurate. The default is half a minute (30 seconds), so that sunrise/sunset 179times will be correct to the minute. 180 181It is useless to set the value smaller than 4*delta, where delta is the 182accuracy in the longitude of the sun (given by the function 183`solar-ecliptic-coordinates') in degrees since (delta/360) x (86400/60) = 4 x 184delta. At present, delta = 0.01 degrees, so the value of the variable 185`solar-error' should be at least 0.04 minutes (about 2.5 seconds)." 186 :type 'number 187 :group 'calendar) 188 189(defvar solar-n-hemi-seasons 190 '("Vernal Equinox" "Summer Solstice" "Autumnal Equinox" "Winter Solstice") 191 "List of season changes for the northern hemisphere.") 192 193(defvar solar-s-hemi-seasons 194 '("Autumnal Equinox" "Winter Solstice" "Vernal Equinox" "Summer Solstice") 195 "List of season changes for the southern hemisphere.") 196 197(defvar solar-sidereal-time-greenwich-midnight 198 nil 199 "Sidereal time at Greenwich at midnight (universal time).") 200 201(defvar solar-northern-spring-or-summer-season nil 202 "Non-nil if northern spring or summer and nil otherwise. 203Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.") 204 205(defun solar-setup () 206 "Prompt user for latitude, longitude, and time zone." 207 (beep) 208 (if (not calendar-longitude) 209 (setq calendar-longitude 210 (solar-get-number 211 "Enter longitude (decimal fraction; + east, - west): "))) 212 (if (not calendar-latitude) 213 (setq calendar-latitude 214 (solar-get-number 215 "Enter latitude (decimal fraction; + north, - south): "))) 216 (if (not calendar-time-zone) 217 (setq calendar-time-zone 218 (solar-get-number 219 "Enter difference from Coordinated Universal Time (in minutes): ")))) 220 221(defun solar-get-number (prompt) 222 "Return a number from the minibuffer, prompting with PROMPT. 223Returns nil if nothing was entered." 224 (let ((x (read-string prompt ""))) 225 (if (not (string-equal x "")) 226 (string-to-number x)))) 227 228;; The condition-case stuff is needed to catch bogus arithmetic 229;; exceptions that occur on some machines (like Sparcs) 230(defun solar-sin-degrees (x) 231 (condition-case nil 232 (sin (degrees-to-radians (mod x 360.0))) 233 (solar-sin-degrees x))) 234(defun solar-cosine-degrees (x) 235 (condition-case nil 236 (cos (degrees-to-radians (mod x 360.0))) 237 (solar-cosine-degrees x))) 238(defun solar-tangent-degrees (x) 239 (condition-case nil 240 (tan (degrees-to-radians (mod x 360.0))) 241 (solar-tangent-degrees x))) 242 243(defun solar-xy-to-quadrant (x y) 244 "Determines the quadrant of the point X, Y." 245 (if (> x 0) 246 (if (> y 0) 1 4) 247 (if (> y 0) 2 3))) 248 249(defun solar-degrees-to-quadrant (angle) 250 "Determines the quadrant of ANGLE." 251 (1+ (floor (mod angle 360) 90))) 252 253(defun solar-arctan (x quad) 254 "Arctangent of X in quadrant QUAD." 255 (let ((deg (radians-to-degrees (atan x)))) 256 (cond ((equal quad 2) (+ deg 180)) 257 ((equal quad 3) (+ deg 180)) 258 ((equal quad 4) (+ deg 360)) 259 (t deg)))) 260 261(defun solar-atn2 (x y) 262 "Arctan of point X, Y." 263 (if (= x 0) 264 (if (> y 0) 90 270) 265 (solar-arctan (/ y x) (solar-xy-to-quadrant x y)))) 266 267(defun solar-arccos (x) 268 "Arcos of X." 269 (let ((y (sqrt (- 1 (* x x))))) 270 (solar-atn2 x y))) 271 272(defun solar-arcsin (y) 273 "Arcsin of Y." 274 (let ((x (sqrt (- 1 (* y y))))) 275 (solar-atn2 x y) 276 )) 277 278(defsubst solar-degrees-to-hours (degrees) 279 "Convert DEGREES to hours." 280 (/ degrees 15.0)) 281 282(defsubst solar-hours-to-days (hour) 283 "Convert HOUR to decimal fraction of a day." 284 (/ hour 24.0)) 285 286(defun solar-right-ascension (longitude obliquity) 287 "Right ascension of the sun, in hours, given LONGITUDE and OBLIQUITY. 288Both arguments are in degrees." 289 (solar-degrees-to-hours 290 (solar-arctan 291 (* (solar-cosine-degrees obliquity) (solar-tangent-degrees longitude)) 292 (solar-degrees-to-quadrant longitude)))) 293 294(defun solar-declination (longitude obliquity) 295 "Declination of the sun, in degrees, given LONGITUDE and OBLIQUITY. 296Both arguments are in degrees." 297 (solar-arcsin 298 (* (solar-sin-degrees obliquity) 299 (solar-sin-degrees longitude)))) 300 301(defun solar-sunrise-and-sunset (time latitude longitude height) 302 "Sunrise, sunset and length of day. 303Parameters are the midday TIME and the LATITUDE, LONGITUDE of the location. 304 305TIME is a pair with the first component being the number of Julian centuries 306elapsed at 0 Universal Time, and the second component being the universal 307time. For instance, the pair corresponding to November 28, 1995 at 16 UT is 308\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between 309Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. 310 311HEIGHT is the angle the center of the sun has over the horizon for the contact 312we are trying to find. For sunrise and sunset, it is usually -0.61 degrees, 313accounting for the edge of the sun being on the horizon. 314 315Coordinates are included because this function is called with latitude=1 316degrees to find out if polar regions have 24 hours of sun or only night." 317 (let* ((rise-time (solar-moment -1 latitude longitude time height)) 318 (set-time (solar-moment 1 latitude longitude time height)) 319 (day-length)) 320 (if (not (and rise-time set-time)) 321 (if (or (and (> latitude 0) 322 solar-northern-spring-or-summer-season) 323 (and (< latitude 0) 324 (not solar-northern-spring-or-summer-season))) 325 (setq day-length 24) 326 (setq day-length 0)) 327 (setq day-length (- set-time rise-time))) 328 (list (if rise-time (+ rise-time (/ calendar-time-zone 60.0)) nil) 329 (if set-time (+ set-time (/ calendar-time-zone 60.0)) nil) 330 day-length))) 331 332(defun solar-moment (direction latitude longitude time height) 333 "Sunrise/sunset at location. 334Sunrise if DIRECTION =-1 or sunset if =1 at LATITUDE, LONGITUDE, with midday 335being TIME. 336 337TIME is a pair with the first component being the number of Julian centuries 338elapsed at 0 Universal Time, and the second component being the universal 339time. For instance, the pair corresponding to November 28, 1995 at 16 UT is 340\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between 341Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. 342 343HEIGHT is the angle the center of the sun has over the horizon for the contact 344we are trying to find. For sunrise and sunset, it is usually -0.61 degrees, 345accounting for the edge of the sun being on the horizon. 346 347Uses binary search." 348 (let* ((ut (car (cdr time))) 349 (possible t) ; we assume that rise or set are possible 350 (utmin (+ ut (* direction 12.0))) 351 (utmax ut) ; the time searched is between utmin and utmax 352 ; utmin and utmax are in hours 353 (utmoment-old 0.0) ; rise or set approximation 354 (utmoment 1.0) ; rise or set approximation 355 (hut 0) ; sun height at utmoment 356 (t0 (car time)) 357 (hmin (car (cdr 358 (solar-horizontal-coordinates (list t0 utmin) 359 latitude longitude t)))) 360 (hmax (car (cdr 361 (solar-horizontal-coordinates (list t0 utmax) 362 latitude longitude t))))) 363 ; -0.61 degrees is the height of the middle of the sun, when it rises 364 ; or sets. 365 (if (< hmin height) 366 (if (> hmax height) 367 (while ;(< i 20) ; we perform a simple dichotomy 368 ; (> (abs (- hut height)) epsilon) 369 (>= (abs (- utmoment utmoment-old)) 370 (/ solar-error 60)) 371 (setq utmoment-old utmoment) 372 (setq utmoment (/ (+ utmin utmax) 2)) 373 (setq hut (car (cdr 374 (solar-horizontal-coordinates 375 (list t0 utmoment) latitude longitude t)))) 376 (if (< hut height) (setq utmin utmoment)) 377 (if (> hut height) (setq utmax utmoment)) 378 ) 379 (setq possible nil)) ; the sun never rises 380 (setq possible nil)) ; the sun never sets 381 (if (not possible) nil utmoment))) 382 383(defun solar-time-string (time time-zone) 384 "Printable form for decimal fraction TIME in TIME-ZONE. 385Format used is given by `calendar-time-display-form'." 386 (let* ((time (round (* 60 time))) 387 (24-hours (/ time 60)) 388 (minutes (format "%02d" (% time 60))) 389 (12-hours (format "%d" (1+ (% (+ 24-hours 11) 12)))) 390 (am-pm (if (>= 24-hours 12) "pm" "am")) 391 (24-hours (format "%02d" 24-hours))) 392 (mapconcat 'eval calendar-time-display-form ""))) 393 394 395(defun solar-daylight (time) 396 "Printable form for time expressed in hours." 397 (format "%d:%02d" 398 (floor time) 399 (floor (* 60 (- time (floor time)))))) 400 401(defun solar-exact-local-noon (date) 402 "Date and Universal Time of local noon at *local date* date. 403 404The date may be different from the one asked for, but it will be the right 405local date. The second component of date should be an integer." 406 (let* ((nd date) 407 (ut (- 12.0 (/ (calendar-longitude) 15))) 408 (te (solar-time-equation date ut))) 409 (setq ut (- ut te)) 410 (if (>= ut 24) 411 (progn 412 (setq nd (list (car date) (+ 1 (car (cdr date))) 413 (car (cdr (cdr date))))) 414 (setq ut (- ut 24)))) 415 (if (< ut 0) 416 (progn 417 (setq nd (list (car date) (- (car (cdr date)) 1) 418 (car (cdr (cdr date))))) 419 (setq ut (+ ut 24)))) 420 (setq nd (calendar-gregorian-from-absolute 421 (calendar-absolute-from-gregorian nd))) 422 ; date standardization 423 (list nd ut))) 424 425(defun solar-sunrise-sunset (date) 426 "List of *local* times of sunrise, sunset, and daylight on Gregorian DATE. 427 428Corresponding value is nil if there is no sunrise/sunset." 429 (let* (; first, get the exact moment of local noon. 430 (exact-local-noon (solar-exact-local-noon date)) 431 ; get the time from the 2000 epoch. 432 (t0 (solar-julian-ut-centuries (car exact-local-noon))) 433 ; store the sidereal time at Greenwich at midnight of UT time. 434 ; find if summer or winter slightly above the equator 435 (equator-rise-set 436 (progn (setq solar-sidereal-time-greenwich-midnight 437 (solar-sidereal-time t0)) 438 (solar-sunrise-and-sunset 439 (list t0 (car (cdr exact-local-noon))) 440 1.0 441 (calendar-longitude) 0))) 442 ; store the spring/summer information, 443 ; compute sunrise and sunset (two first components of rise-set). 444 ; length of day is the third component (it is only the difference 445 ; between sunset and sunrise when there is a sunset and a sunrise) 446 (rise-set 447 (progn 448 (setq solar-northern-spring-or-summer-season 449 (if (> (car (cdr (cdr equator-rise-set))) 12) t nil)) 450 (solar-sunrise-and-sunset 451 (list t0 (car (cdr exact-local-noon))) 452 (calendar-latitude) 453 (calendar-longitude) -0.61))) 454 (rise (car rise-set)) 455 (adj-rise (if rise (dst-adjust-time date rise) nil)) 456 (set (car (cdr rise-set))) 457 (adj-set (if set (dst-adjust-time date set) nil)) 458 (length (car (cdr (cdr rise-set)))) ) 459 (list 460 (and rise (calendar-date-equal date (car adj-rise)) (cdr adj-rise)) 461 (and set (calendar-date-equal date (car adj-set)) (cdr adj-set)) 462 (solar-daylight length)))) 463 464(defun solar-sunrise-sunset-string (date) 465 "String of *local* times of sunrise, sunset, and daylight on Gregorian DATE." 466 (let ((l (solar-sunrise-sunset date))) 467 (format 468 "%s, %s at %s (%s hours daylight)" 469 (if (car l) 470 (concat "Sunrise " (apply 'solar-time-string (car l))) 471 "No sunrise") 472 (if (car (cdr l)) 473 (concat "sunset " (apply 'solar-time-string (car (cdr l)))) 474 "no sunset") 475 (eval calendar-location-name) 476 (car (cdr (cdr l)))))) 477 478(defun solar-julian-ut-centuries (date) 479 "Number of Julian centuries elapsed since 1 Jan, 2000 at noon U.T. for Gregorian DATE." 480 (/ (- (calendar-absolute-from-gregorian date) 481 (calendar-absolute-from-gregorian '(1 1.5 2000))) 482 36525.0)) 483 484(defun solar-ephemeris-time(time) 485 "Ephemeris Time at moment TIME. 486 487TIME is a pair with the first component being the number of Julian centuries 488elapsed at 0 Universal Time, and the second component being the universal 489time. For instance, the pair corresponding to November 28, 1995 at 16 UT is 490\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between 491Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. 492 493Result is in julian centuries of ephemeris time." 494 (let* ((t0 (car time)) 495 (ut (car (cdr time))) 496 (t1 (+ t0 (/ (/ ut 24.0) 36525))) 497 (y (+ 2000 (* 100 t1))) 498 (dt (* 86400 (solar-ephemeris-correction (floor y))))) 499 (+ t1 (/ (/ dt 86400) 36525)))) 500 501(defun solar-date-next-longitude (d l) 502 "First moment on or after Julian day number D when sun's longitude is a 503multiple of L degrees at calendar-location-name with that location's 504local time (including any daylight saving rules). 505 506L must be an integer divisor of 360. 507 508Result is in local time expressed astronomical (Julian) day numbers. 509 510The values of calendar-daylight-savings-starts, 511calendar-daylight-savings-starts-time, calendar-daylight-savings-ends, 512calendar-daylight-savings-ends-time, calendar-daylight-time-offset, and 513calendar-time-zone are used to interpret local time." 514 (let* ((long) 515 (start d) 516 (start-long (solar-longitude d)) 517 (next (mod (* l (1+ (floor (/ start-long l)))) 360)) 518 (end (+ d (* (/ l 360.0) 400))) 519 (end-long (solar-longitude end))) 520 (while ;; bisection search for nearest minute 521 (< 0.00001 (- end start)) 522 ;; start <= d < end 523 ;; start-long <= next < end-long when next != 0 524 ;; when next = 0, we look for the discontinuity (start-long is near 360 525 ;; and end-long is small (less than l). 526 (setq d (/ (+ start end) 2.0)) 527 (setq long (solar-longitude d)) 528 (if (or (and (/= next 0) (< long next)) 529 (and (= next 0) (< l long))) 530 (progn 531 (setq start d) 532 (setq start-long long)) 533 (setq end d) 534 (setq end-long long))) 535 (/ (+ start end) 2.0))) 536 537(defun solar-horizontal-coordinates 538 (time latitude longitude for-sunrise-sunset) 539 "Azimuth and height of the sun at TIME, LATITUDE, and LONGITUDE. 540 541TIME is a pair with the first component being the number of Julian centuries 542elapsed at 0 Universal Time, and the second component being the universal 543time. For instance, the pair corresponding to November 28, 1995 at 16 UT is 544\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between 545Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. 546 547The azimuth is given in degrees as well as the height (between -180 and 180)." 548 (let* ((ut (car (cdr time))) 549 (ec (solar-equatorial-coordinates time for-sunrise-sunset)) 550 (st (+ solar-sidereal-time-greenwich-midnight 551 (* ut 1.00273790935))) 552 (ah (- (* st 15) (* 15 (car ec)) (* -1 (calendar-longitude)))) 553 ; hour angle (in degrees) 554 (de (car (cdr ec))) 555 (azimuth (solar-atn2 (- (* (solar-cosine-degrees ah) 556 (solar-sin-degrees latitude)) 557 (* (solar-tangent-degrees de) 558 (solar-cosine-degrees latitude))) 559 (solar-sin-degrees ah))) 560 (height (solar-arcsin 561 (+ (* (solar-sin-degrees latitude) (solar-sin-degrees de)) 562 (* (solar-cosine-degrees latitude) 563 (solar-cosine-degrees de) 564 (solar-cosine-degrees ah)))))) 565 (if (> height 180) (setq height (- height 360))) 566 (list azimuth height))) 567 568(defun solar-equatorial-coordinates (time for-sunrise-sunset) 569 "Right ascension (in hours) and declination (in degrees) of the sun at TIME. 570 571TIME is a pair with the first component being the number of Julian centuries 572elapsed at 0 Universal Time, and the second component being the universal 573time. For instance, the pair corresponding to November 28, 1995 at 16 UT is 574\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between 575Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT." 576 (let* ((tm (solar-ephemeris-time time)) 577 (ec (solar-ecliptic-coordinates tm for-sunrise-sunset))) 578 (list (solar-right-ascension (car ec) (car (cdr ec))) 579 (solar-declination (car ec) (car (cdr ec)))))) 580 581(defun solar-ecliptic-coordinates (time for-sunrise-sunset) 582 "Apparent longitude of the sun, ecliptic inclination, (both in degrees) 583equation of time (in hours) and nutation in longitude (in seconds) 584at moment `time', expressed in julian centuries of Ephemeris Time 585since January 1st, 2000, at 12 ET." 586 (let* ((l (+ 280.46645 587 (* 36000.76983 time) 588 (* 0.0003032 time time))) ; sun mean longitude 589 (ml (+ 218.3165 590 (* 481267.8813 time))) ; moon mean longitude 591 (m (+ 357.52910 592 (* 35999.05030 time) 593 (* -0.0001559 time time) 594 (* -0.00000048 time time time))) ; sun mean anomaly 595 (i (+ 23.43929111 (* -0.013004167 time) 596 (* -0.00000016389 time time) 597 (* 0.0000005036 time time time))); mean inclination 598 (c (+ (* (+ 1.914600 599 (* -0.004817 time) 600 (* -0.000014 time time)) 601 (solar-sin-degrees m)) 602 (* (+ 0.019993 (* -0.000101 time)) 603 (solar-sin-degrees (* 2 m))) 604 (* 0.000290 605 (solar-sin-degrees (* 3 m))))) ; center equation 606 (L (+ l c)) ; total longitude 607 (omega (+ 125.04 608 (* -1934.136 time))) ; longitude of moon's ascending node 609 ; on the ecliptic 610 (nut (if (not for-sunrise-sunset) 611 (+ (* -17.20 (solar-sin-degrees omega)) 612 (* -1.32 (solar-sin-degrees (* 2 l))) 613 (* -0.23 (solar-sin-degrees (* 2 ml))) 614 (* 0.21 (solar-sin-degrees (* 2 omega)))) 615 nil)) 616 ; nut = nutation in longitude, measured in seconds of angle. 617 (ecc (if (not for-sunrise-sunset) 618 (+ 0.016708617 619 (* -0.000042037 time) 620 (* -0.0000001236 time time)) ; eccentricity of earth's orbit 621 nil)) 622 (app (+ L 623 -0.00569 624 (* -0.00478 625 (solar-sin-degrees omega)))) ; apparent longitude of sun 626 (y (if (not for-sunrise-sunset) 627 (* (solar-tangent-degrees (/ i 2)) 628 (solar-tangent-degrees (/ i 2))) 629 nil)) 630 (time-eq (if (not for-sunrise-sunset) 631 (/ (* 12 (+ (* y (solar-sin-degrees (* 2 l))) 632 (* -2 ecc (solar-sin-degrees m)) 633 (* 4 ecc y (solar-sin-degrees m) 634 (solar-cosine-degrees (* 2 l))) 635 (* -0.5 y y (solar-sin-degrees (* 4 l))) 636 (* -1.25 ecc ecc (solar-sin-degrees (* 2 m))))) 637 3.1415926535) 638 nil))) 639 ; equation of time, in hours 640 (list app i time-eq nut))) 641 642(defconst solar-data-list 643 '((403406 4.721964 1.621043) 644 (195207 5.937458 62830.348067) 645 (119433 1.115589 62830.821524) 646 (112392 5.781616 62829.634302) 647 (3891 5.5474 125660.5691) 648 (2819 1.5120 125660.984) 649 (1721 4.1897 62832.4766) 650 (0 1.163 0.813) 651 (660 5.415 125659.31) 652 (350 4.315 57533.85) 653 (334 4.553 -33.931) 654 (314 5.198 777137.715) 655 (268 5.989 78604.191) 656 (242 2.911 5.412) 657 (234 1.423 39302.098) 658 (158 0.061 -34.861) 659 (132 2.317 115067.698) 660 (129 3.193 15774.337) 661 (114 2.828 5296.670) 662 (99 0.52 58849.27) 663 (93 4.65 5296.11) 664 (86 4.35 -3980.70) 665 (78 2.75 52237.69) 666 (72 4.50 55076.47) 667 (68 3.23 261.08) 668 (64 1.22 15773.85) 669 (46 0.14 188491.03) 670 (38 3.44 -7756.55) 671 (37 4.37 264.89) 672 (32 1.14 117906.27) 673 (29 2.84 55075.75) 674 (28 5.96 -7961.39) 675 (27 5.09 188489.81) 676 (27 1.72 2132.19) 677 (25 2.56 109771.03) 678 (24 1.92 54868.56) 679 (21 0.09 25443.93) 680 (21 5.98 -55731.43) 681 (20 4.03 60697.74) 682 (18 4.47 2132.79) 683 (17 0.79 109771.63) 684 (14 4.24 -7752.82) 685 (13 2.01 188491.91) 686 (13 2.65 207.81) 687 (13 4.98 29424.63) 688 (12 0.93 -7.99) 689 (10 2.21 46941.14) 690 (10 3.59 -68.29) 691 (10 1.50 21463.25) 692 (10 2.55 157208.40))) 693 694(defun solar-longitude (d) 695 "Longitude of sun on astronomical (Julian) day number D. 696Accurary is about 0.0006 degree (about 365.25*24*60*0.0006/360 = 1 minutes). 697 698The values of calendar-daylight-savings-starts, 699calendar-daylight-savings-starts-time, calendar-daylight-savings-ends, 700calendar-daylight-savings-ends-time, calendar-daylight-time-offset, and 701calendar-time-zone are used to interpret local time." 702 (let* ((a-d (calendar-absolute-from-astro d)) 703 ;; get Universal Time 704 (date (calendar-astro-from-absolute 705 (- a-d 706 (if (dst-in-effect a-d) 707 (/ calendar-daylight-time-offset 24.0 60.0) 0) 708 (/ calendar-time-zone 60.0 24.0)))) 709 ;; get Ephemeris Time 710 (date (+ date (solar-ephemeris-correction 711 (extract-calendar-year 712 (calendar-gregorian-from-absolute 713 (floor 714 (calendar-absolute-from-astro 715 date))))))) 716 (U (/ (- date 2451545) 3652500)) 717 (longitude 718 (+ 4.9353929 719 (* 62833.1961680 U) 720 (* 0.0000001 721 (apply '+ 722 (mapcar '(lambda (x) 723 (* (car x) 724 (sin (mod 725 (+ (car (cdr x)) 726 (* (car (cdr (cdr x))) U)) 727 (* 2 pi))))) 728 solar-data-list))))) 729 (aberration 730 (* 0.0000001 (- (* 17 (cos (+ 3.10 (* 62830.14 U)))) 973))) 731 (A1 (mod (+ 2.18 (* U (+ -3375.70 (* 0.36 U)))) (* 2 pi))) 732 (A2 (mod (+ 3.51 (* U (+ 125666.39 (* 0.10 U)))) (* 2 pi))) 733 (nutation (* -0.0000001 (+ (* 834 (sin A1)) (* 64 (sin A2)))))) 734 (mod (radians-to-degrees (+ longitude aberration nutation)) 360.0))) 735 736(defun solar-ephemeris-correction (year) 737 "Ephemeris time minus Universal Time during Gregorian year. 738Result is in days. 739 740For the years 1800-1987, the maximum error is 1.9 seconds. 741For the other years, the maximum error is about 30 seconds." 742 (cond ((and (<= 1988 year) (< year 2020)) 743 (/ (+ year -2000 67.0) 60.0 60.0 24.0)) 744 ((and (<= 1900 year) (< year 1988)) 745 (let* ((theta (/ (- (calendar-astro-from-absolute 746 (calendar-absolute-from-gregorian 747 (list 7 1 year))) 748 (calendar-astro-from-absolute 749 (calendar-absolute-from-gregorian 750 '(1 1 1900)))) 751 36525.0)) 752 (theta2 (* theta theta)) 753 (theta3 (* theta2 theta)) 754 (theta4 (* theta2 theta2)) 755 (theta5 (* theta3 theta2))) 756 (+ -0.00002 757 (* 0.000297 theta) 758 (* 0.025184 theta2) 759 (* -0.181133 theta3) 760 (* 0.553040 theta4) 761 (* -0.861938 theta5) 762 (* 0.677066 theta3 theta3) 763 (* -0.212591 theta4 theta3)))) 764 ((and (<= 1800 year) (< year 1900)) 765 (let* ((theta (/ (- (calendar-astro-from-absolute 766 (calendar-absolute-from-gregorian 767 (list 7 1 year))) 768 (calendar-astro-from-absolute 769 (calendar-absolute-from-gregorian 770 '(1 1 1900)))) 771 36525.0)) 772 (theta2 (* theta theta)) 773 (theta3 (* theta2 theta)) 774 (theta4 (* theta2 theta2)) 775 (theta5 (* theta3 theta2))) 776 (+ -0.000009 777 (* 0.003844 theta) 778 (* 0.083563 theta2) 779 (* 0.865736 theta3) 780 (* 4.867575 theta4) 781 (* 15.845535 theta5) 782 (* 31.332267 theta3 theta3) 783 (* 38.291999 theta4 theta3) 784 (* 28.316289 theta4 theta4) 785 (* 11.636204 theta4 theta5) 786 (* 2.043794 theta5 theta5)))) 787 ((and (<= 1620 year) (< year 1800)) 788 (let ((x (/ (- year 1600) 10.0))) 789 (/ (+ (* 2.19167 x x) (* -40.675 x) 196.58333) 60.0 60.0 24.0))) 790 (t (let* ((tmp (- (calendar-astro-from-absolute 791 (calendar-absolute-from-gregorian 792 (list 1 1 year))) 793 2382148)) 794 (second (- (/ (* tmp tmp) 41048480.0) 15))) 795 (/ second 60.0 60.0 24.0))))) 796 797(defun solar-sidereal-time (t0) 798 "Sidereal time (in hours) in Greenwich. 799 800At T0=Julian centuries of universal time. 801T0 must correspond to 0 hours UT." 802 (let* ((mean-sid-time (+ 6.6973746 803 (* 2400.051337 t0) 804 (* 0.0000258622 t0 t0) 805 (* -0.0000000017222 t0 t0 t0))) 806 (et (solar-ephemeris-time (list t0 0.0))) 807 (nut-i (solar-ecliptic-coordinates et nil)) 808 (nut (car (cdr (cdr (cdr nut-i))))) ; nutation 809 (i (car (cdr nut-i)))) ; inclination 810 (mod (+ (mod (+ mean-sid-time 811 (/ (/ (* nut (solar-cosine-degrees i)) 15) 3600)) 24.0) 812 24.0) 813 24.0))) 814 815(defun solar-time-equation (date ut) 816 "Equation of time expressed in hours at Gregorian DATE at Universal time UT." 817 (let* ((et (solar-date-to-et date ut)) 818 (ec (solar-ecliptic-coordinates et nil))) 819 (car (cdr (cdr ec))))) 820 821(defun solar-date-to-et (date ut) 822 "Ephemeris Time at Gregorian DATE at Universal Time UT (in hours). 823Expressed in julian centuries of Ephemeris Time." 824 (let ((t0 (solar-julian-ut-centuries date))) 825 (solar-ephemeris-time (list t0 ut)))) 826 827;;;###autoload 828(defun sunrise-sunset (&optional arg) 829 "Local time of sunrise and sunset for today. Accurate to a few seconds. 830If called with an optional prefix argument, prompt for date. 831 832If called with an optional double prefix argument, prompt for longitude, 833latitude, time zone, and date, and always use standard time. 834 835This function is suitable for execution in a .emacs file." 836 (interactive "p") 837 (or arg (setq arg 1)) 838 (if (and (< arg 16) 839 (not (and calendar-latitude calendar-longitude calendar-time-zone))) 840 (solar-setup)) 841 (let* ((calendar-longitude 842 (if (< arg 16) calendar-longitude 843 (solar-get-number 844 "Enter longitude (decimal fraction; + east, - west): "))) 845 (calendar-latitude 846 (if (< arg 16) calendar-latitude 847 (solar-get-number 848 "Enter latitude (decimal fraction; + north, - south): "))) 849 (calendar-time-zone 850 (if (< arg 16) calendar-time-zone 851 (solar-get-number 852 "Enter difference from Coordinated Universal Time (in minutes): "))) 853 (calendar-location-name 854 (if (< arg 16) calendar-location-name 855 (let ((float-output-format "%.1f")) 856 (format "%s%s, %s%s" 857 (if (numberp calendar-latitude) 858 (abs calendar-latitude) 859 (+ (aref calendar-latitude 0) 860 (/ (aref calendar-latitude 1) 60.0))) 861 (if (numberp calendar-latitude) 862 (if (> calendar-latitude 0) "N" "S") 863 (if (equal (aref calendar-latitude 2) 'north) "N" "S")) 864 (if (numberp calendar-longitude) 865 (abs calendar-longitude) 866 (+ (aref calendar-longitude 0) 867 (/ (aref calendar-longitude 1) 60.0))) 868 (if (numberp calendar-longitude) 869 (if (> calendar-longitude 0) "E" "W") 870 (if (equal (aref calendar-longitude 2) 'east) 871 "E" "W")))))) 872 (calendar-standard-time-zone-name 873 (if (< arg 16) calendar-standard-time-zone-name 874 (cond ((= calendar-time-zone 0) "UTC") 875 ((< calendar-time-zone 0) 876 (format "UTC%dmin" calendar-time-zone)) 877 (t (format "UTC+%dmin" calendar-time-zone))))) 878 (calendar-daylight-savings-starts 879 (if (< arg 16) calendar-daylight-savings-starts)) 880 (calendar-daylight-savings-ends 881 (if (< arg 16) calendar-daylight-savings-ends)) 882 (date (if (< arg 4) (calendar-current-date) (calendar-read-date))) 883 (date-string (calendar-date-string date t)) 884 (time-string (solar-sunrise-sunset-string date)) 885 (msg (format "%s: %s" date-string time-string)) 886 (one-window (one-window-p t))) 887 (if (<= (length msg) (frame-width)) 888 (message "%s" msg) 889 (with-output-to-temp-buffer "*temp*" 890 (princ (concat date-string "\n" time-string))) 891 (message "%s" 892 (substitute-command-keys 893 (if one-window 894 (if pop-up-windows 895 "Type \\[delete-other-windows] to remove temp window." 896 "Type \\[switch-to-buffer] RET to remove temp window.") 897 "Type \\[switch-to-buffer-other-window] RET to restore old contents of temp window.")))))) 898 899(defun calendar-sunrise-sunset () 900 "Local time of sunrise and sunset for date under cursor. 901Accurate to a few seconds." 902 (interactive) 903 (if (not (and calendar-latitude calendar-longitude calendar-time-zone)) 904 (solar-setup)) 905 (let ((date (calendar-cursor-to-date t))) 906 (message "%s: %s" 907 (calendar-date-string date t t) 908 (solar-sunrise-sunset-string date)))) 909 910(defun diary-sunrise-sunset () 911 "Local time of sunrise and sunset as a diary entry. 912Accurate to a few seconds." 913 (if (not (and calendar-latitude calendar-longitude calendar-time-zone)) 914 (solar-setup)) 915 (solar-sunrise-sunset-string date)) 916 917(defcustom diary-sabbath-candles-minutes 18 918 "*Number of minutes before sunset for sabbath candle lighting." 919 :group 'diary 920 :type 'integer 921 :version "21.1") 922 923(defun diary-sabbath-candles (&optional mark) 924 "Local time of candle lighting diary entry--applies if date is a Friday. 925No diary entry if there is no sunset on that date. 926 927An optional parameter MARK specifies a face or single-character string to 928use when highlighting the day in the calendar." 929 (if (not (and calendar-latitude calendar-longitude calendar-time-zone)) 930 (solar-setup)) 931 (if (= (% (calendar-absolute-from-gregorian date) 7) 5);; Friday 932 (let* ((sunset (car (cdr (solar-sunrise-sunset date)))) 933 (light (if sunset 934 (cons (- (car sunset) 935 (/ diary-sabbath-candles-minutes 60.0)) 936 (cdr sunset))))) 937 (if sunset 938 (cons mark 939 (format "%s Sabbath candle lighting" 940 (apply 'solar-time-string light))))))) 941 942; from Meeus, 1991, page 167 943(defconst solar-seasons-data 944 '((485 324.96 1934.136) 945 (203 337.23 32964.467) 946 (199 342.08 20.186) 947 (182 27.85 445267.112) 948 (156 73.14 45036.886) 949 (136 171.52 22518.443) 950 (77 222.54 65928.934) 951 (74 296.72 3034.906) 952 (70 243.58 9037.513) 953 (58 119.81 33718.147) 954 (52 297.17 150.678) 955 (50 21.02 2281.226) 956 (45 247.54 29929.562) 957 (44 325.15 31555.956) 958 (29 60.93 4443.417) 959 (18 155.12 67555.328) 960 (17 288.79 4562.452) 961 (16 198.04 62894.029) 962 (14 199.76 31436.921) 963 (12 95.39 14577.848) 964 (12 287.11 31931.756) 965 (12 320.81 34777.259) 966 (9 227.73 1222.114) 967 (8 15.45 16859.074))) 968 969(defun solar-equinoxes/solstices (k year) 970 "Date of equinox/solstice K for YEAR. 971K=0, spring equinox; K=1, summer solstice; K=2, fall equinox; 972K=3, winter solstice. 973RESULT is a gregorian local date. 974 975Accurate to less than a minute between 1951 and 2050." 976 (let* ((JDE0 (solar-mean-equinoxes/solstices k year)) 977 (T (/ (- JDE0 2451545.0) 36525)) 978 (W (- (* 35999.373 T) 2.47)) 979 (Delta-lambda (+ 1 (* 0.0334 (solar-cosine-degrees W)) 980 (* 0.0007 (solar-cosine-degrees (* 2 W))))) 981 (S (apply '+ (mapcar '(lambda(x) 982 (* (car x) (solar-cosine-degrees 983 (+ (* (car (cdr (cdr x))) T) 984 (car (cdr x)))))) 985 solar-seasons-data))) 986 (JDE (+ JDE0 (/ (* 0.00001 S) Delta-lambda))) 987 (correction (+ 102.3 (* 123.5 T) (* 32.5 T T))) 988 ; ephemeris time correction 989 (JD (- JDE (/ correction 86400))) 990 (date (calendar-gregorian-from-absolute (floor (- JD 1721424.5)))) 991 (time (- (- JD 0.5) (floor (- JD 0.5)))) 992 ) 993 (list (car date) (+ (car (cdr date)) time 994 (/ (/ calendar-time-zone 60.0) 24.0)) 995 (car (cdr (cdr date)))))) 996 997; from Meeus, 1991, page 166 998(defun solar-mean-equinoxes/solstices (k year) 999 "Julian day of mean equinox/solstice K for YEAR. 1000K=0, spring equinox; K=1, summer solstice; K=2, fall equinox; K=3, winter 1001solstice. These formulas are only to be used between 1000 BC and 3000 AD." 1002 (let ((y (/ year 1000.0)) 1003 (z (/ (- year 2000) 1000.0))) 1004 (if (< year 1000) ; actually between -1000 and 1000 1005 (cond ((equal k 0) (+ 1721139.29189 1006 (* 365242.13740 y) 1007 (* 0.06134 y y) 1008 (* 0.00111 y y y) 1009 (* -0.00071 y y y y))) 1010 ((equal k 1) (+ 1721233.25401 1011 (* 365241.72562 y) 1012 (* -0.05323 y y) 1013 (* 0.00907 y y y) 1014 (* 0.00025 y y y y))) 1015 ((equal k 2) (+ 1721325.70455 1016 (* 365242.49558 y) 1017 (* -0.11677 y y) 1018 (* -0.00297 y y y) 1019 (* 0.00074 y y y y))) 1020 ((equal k 3) (+ 1721414.39987 1021 (* 365242.88257 y) 1022 (* -0.00769 y y) 1023 (* -0.00933 y y y) 1024 (* -0.00006 y y y y)))) 1025 ; actually between 1000 and 3000 1026 (cond ((equal k 0) (+ 2451623.80984 1027 (* 365242.37404 z) 1028 (* 0.05169 z z) 1029 (* -0.00411 z z z) 1030 (* -0.00057 z z z z))) 1031 ((equal k 1) (+ 2451716.56767 1032 (* 365241.62603 z) 1033 (* 0.00325 z z) 1034 (* 0.00888 z z z) 1035 (* -0.00030 z z z z))) 1036 ((equal k 2) (+ 2451810.21715 1037 (* 365242.01767 z) 1038 (* -0.11575 z z) 1039 (* 0.00337 z z z) 1040 (* 0.00078 z z z z))) 1041 ((equal k 3) (+ 2451900.05952 1042 (* 365242.74049 z) 1043 (* -0.06223 z z) 1044 (* -0.00823 z z z) 1045 (* 0.00032 z z z z))))))) 1046 1047;;;###autoload 1048(defun solar-equinoxes-solstices () 1049 "*local* date and time of equinoxes and solstices, if visible in the calendar window. 1050Requires floating point." 1051 (let ((m displayed-month) 1052 (y displayed-year)) 1053 (increment-calendar-month m y (cond ((= 1 (% m 3)) -1) 1054 ((= 2 (% m 3)) 1) 1055 (t 0))) 1056 (let* ((calendar-standard-time-zone-name 1057 (if calendar-time-zone calendar-standard-time-zone-name "UTC")) 1058 (calendar-daylight-savings-starts 1059 (if calendar-time-zone calendar-daylight-savings-starts)) 1060 (calendar-daylight-savings-ends 1061 (if calendar-time-zone calendar-daylight-savings-ends)) 1062 (calendar-time-zone (if calendar-time-zone calendar-time-zone 0)) 1063 (k (1- (/ m 3))) 1064 (d0 (solar-equinoxes/solstices k y)) 1065 (d1 (list (car d0) (floor (car (cdr d0))) (car (cdr (cdr d0))))) 1066 (h0 (* 24 (- (car (cdr d0)) (floor (car (cdr d0)))))) 1067 (adj (dst-adjust-time d1 h0)) 1068 (d (list (car (car adj)) 1069 (+ (car (cdr (car adj)) ) 1070 (/ (car (cdr adj)) 24.0)) 1071 (car (cdr (cdr (car adj)))))) 1072 ; The following is nearly as accurate, but not quite: 1073 ;(d0 (solar-date-next-longitude 1074 ; (calendar-astro-from-absolute 1075 ; (calendar-absolute-from-gregorian 1076 ; (list (+ 3 (* k 3)) 15 y))) 1077 ; 90)) 1078 ;(abs-day (calendar-absolute-from-astro d))) 1079 (abs-day (calendar-absolute-from-gregorian d))) 1080 (list 1081 (list (calendar-gregorian-from-absolute (floor abs-day)) 1082 (format "%s %s" 1083 (nth k (if (and calendar-latitude 1084 (< (calendar-latitude) 0)) 1085 solar-s-hemi-seasons 1086 solar-n-hemi-seasons)) 1087 (solar-time-string 1088 (* 24 (- abs-day (floor abs-day))) 1089 (if (dst-in-effect abs-day) 1090 calendar-daylight-time-zone-name 1091 calendar-standard-time-zone-name)))))))) 1092 1093 1094(provide 'solar) 1095 1096;;; arch-tag: bc0ff693-df58-4666-bde4-2a7837ccb8fe 1097;;; solar.el ends here 1098