solar.el 44.3 KB
Newer Older
1
;;; solar.el --- calendar functions for solar events
Jim Blandy's avatar
Jim Blandy committed
2

3
;; Copyright (C) 1992, 1993, 1995, 1997, 2003 Free Software Foundation, Inc.
Jim Blandy's avatar
Jim Blandy committed
4 5

;; Author: Edward M. Reingold <reingold@cs.uiuc.edu>
6
;;	Denis B. Roegel <Denis.Roegel@loria.fr>
Eric S. Raymond's avatar
Eric S. Raymond committed
7 8 9
;; Keywords: calendar
;; Human-Keywords: sunrise, sunset, equinox, solstice, calendar, diary,
;;	holidays
Jim Blandy's avatar
Jim Blandy committed
10 11 12

;; This file is part of GNU Emacs.

13 14 15 16 17
;; GNU Emacs is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2, or (at your option)
;; any later version.

Jim Blandy's avatar
Jim Blandy committed
18
;; GNU Emacs is distributed in the hope that it will be useful,
19 20 21 22 23
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.

;; You should have received a copy of the GNU General Public License
Erik Naggum's avatar
Erik Naggum committed
24 25 26
;; along with GNU Emacs; see the file COPYING.  If not, write to the
;; Free Software Foundation, Inc., 59 Temple Place - Suite 330,
;; Boston, MA 02111-1307, USA.
Jim Blandy's avatar
Jim Blandy committed
27 28 29

;;; Commentary:

30 31
;; This collection of functions implements the features of calendar.el,
;; diary.el, and holiday.el that deal with times of day, sunrise/sunset, and
Paul Eggert's avatar
Paul Eggert committed
32
;; equinoxes/solstices.
Jim Blandy's avatar
Jim Blandy committed
33 34

;; Based on the ``Almanac for Computers 1984,'' prepared by the Nautical
35 36
;; Almanac Office, United States Naval Observatory, Washington, 1984, on
;; ``Astronomical Formulae for Calculators,'' 3rd ed., by Jean Meeus,
37 38 39 40
;; 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.
41

Jim Blandy's avatar
Jim Blandy committed
42
;;
43 44 45
;; 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.
Jim Blandy's avatar
Jim Blandy committed
46
;;
47 48
;;    2. Equinox/solstice times will be accurate to the minute for years
;;       1951--2050.  For other years the times will be within +/- 1 minute.
Jim Blandy's avatar
Jim Blandy committed
49

50
;; Technical details of all the calendrical calculations can be found in
51 52
;; ``Calendrical Calculations: The Millennium Edition'' by Edward M. Reingold
;; and Nachum Dershowitz, Cambridge University Press (2001).
53

Jim Blandy's avatar
Jim Blandy committed
54 55 56 57 58 59 60 61
;; Comments, corrections, and improvements should be sent to
;;  Edward M. Reingold               Department of Computer Science
;;  (217) 333-6733                   University of Illinois at Urbana-Champaign
;;  reingold@cs.uiuc.edu             1304 West Springfield Avenue
;;                                   Urbana, Illinois 61801

;;; Code:

62 63 64
(defvar displayed-month)
(defvar displayed-year)

Jim Blandy's avatar
Jim Blandy committed
65 66
(if (fboundp 'atan)
    (require 'lisp-float-type)
67
  (error "Solar/lunar calculations impossible since floating point is unavailable"))
Jim Blandy's avatar
Jim Blandy committed
68

69
(require 'cal-dst)
70
(require 'cal-julian)
71 72

;;;###autoload
Richard M. Stallman's avatar
Richard M. Stallman committed
73
(defcustom calendar-time-display-form
74 75 76 77 78
  '(12-hours ":" minutes am-pm
    (if time-zone " (") time-zone (if time-zone ")"))
  "*The pseudo-pattern that governs the way a time of day is formatted.

A pseudo-pattern is a list of expressions that can involve the keywords
79 80
`12-hours', `24-hours', and `minutes', all numbers in string form,
and `am-pm' and `time-zone', both alphabetic strings.
81 82 83 84 85 86

For example, the form

  '(24-hours \":\" minutes
    (if time-zone \" (\") time-zone (if time-zone \")\"))

Richard M. Stallman's avatar
Richard M. Stallman committed
87 88 89
would give military-style times like `21:07 (UTC)'."
  :type 'sexp
  :group 'calendar)
90 91

;;;###autoload
Richard M. Stallman's avatar
Richard M. Stallman committed
92
(defcustom calendar-latitude nil
93 94 95 96 97 98 99
  "*Latitude of `calendar-location-name' in degrees.

The value can be either a decimal fraction (one place of accuracy is
sufficient), + north, - south, such as 40.7 for New York City, or the value
can be a vector [degrees minutes north/south] such as [40 50 north] for New
York City.

Richard M. Stallman's avatar
Richard M. Stallman committed
100 101 102 103 104 105 106 107 108 109
This variable should be set in `site-start'.el."
  :type '(choice (const nil)
		 (number :tag "Exact")
		 (vector :value [0 0 north]
			 (integer :tag "Degrees")
			 (integer :tag "Minutes")
			 (choice :tag "Position"
				 (const north)
				 (const south))))
  :group 'calendar)
110 111

;;;###autoload
Richard M. Stallman's avatar
Richard M. Stallman committed
112
(defcustom calendar-longitude nil
113 114 115 116 117 118 119
  "*Longitude of `calendar-location-name' in degrees.

The value can be either a decimal fraction (one place of accuracy is
sufficient), + east, - west, such as -73.9 for New York City, or the value
can be a vector [degrees minutes east/west] such as [73 55 west] for New
York City.

Richard M. Stallman's avatar
Richard M. Stallman committed
120 121 122 123 124 125 126 127 128 129
This variable should be set in `site-start'.el."
  :type '(choice (const nil)
		 (number :tag "Exact")
		 (vector :value [0 0 west]
			 (integer :tag "Degrees")
			 (integer :tag "Minutes")
			 (choice :tag "Position"
				 (const east)
				 (const west))))
  :group 'calendar)
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149

(defsubst calendar-latitude ()
  "Convert calendar-latitude to a signed decimal fraction, if needed."
  (if (numberp calendar-latitude)
      calendar-latitude
    (let ((lat (+ (aref calendar-latitude 0)
                  (/ (aref calendar-latitude 1) 60.0))))
      (if (equal (aref calendar-latitude 2) 'north)
          lat
        (- lat)))))

(defsubst calendar-longitude ()
  "Convert calendar-longitude to a signed decimal fraction, if needed."
  (if (numberp calendar-longitude)
      calendar-longitude
    (let ((long (+ (aref calendar-longitude 0)
                  (/ (aref calendar-longitude 1) 60.0))))
      (if (equal (aref calendar-longitude 2) 'east)
          long
        (- long)))))
150 151

;;;###autoload
Richard M. Stallman's avatar
Richard M. Stallman committed
152
(defcustom calendar-location-name
153 154
  '(let ((float-output-format "%.1f"))
     (format "%s%s, %s%s"
155 156 157 158 159 160 161 162 163 164 165 166 167
             (if (numberp calendar-latitude)
                 (abs calendar-latitude)
               (+ (aref calendar-latitude 0)
                  (/ (aref calendar-latitude 1) 60.0)))
             (if (numberp calendar-latitude)
                 (if (> calendar-latitude 0) "N" "S")
               (if (equal (aref calendar-latitude 2) 'north) "N" "S"))
             (if (numberp calendar-longitude)
                 (abs calendar-longitude)
               (+ (aref calendar-longitude 0)
                  (/ (aref calendar-longitude 1) 60.0)))
             (if (numberp calendar-longitude)
                 (if (> calendar-longitude 0) "E" "W")
168
               (if (equal (aref calendar-longitude 2) 'east) "E" "W"))))
169
  "*Expression evaluating to name of `calendar-longitude', `calendar-latitude'.
170 171 172
For example, \"New York City\".  Default value is just the latitude, longitude
pair.

Richard M. Stallman's avatar
Richard M. Stallman committed
173 174 175
This variable should be set in `site-start'.el."
  :type 'sexp
  :group 'calendar)
176

Richard M. Stallman's avatar
Richard M. Stallman committed
177
(defcustom solar-error 0.5
178 179 180 181 182 183 184 185 186 187
"*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
Richard M. Stallman's avatar
Richard M. Stallman committed
188 189 190
`solar-error' should be at least 0.04 minutes (about 2.5 seconds)."
  :type 'number
  :group 'calendar)
Jim Blandy's avatar
Jim Blandy committed
191

192
(defvar solar-n-hemi-seasons
193 194 195
  '("Vernal Equinox" "Summer Solstice" "Autumnal Equinox" "Winter Solstice")
  "List of season changes for the northern hemisphere.")

196
(defvar solar-s-hemi-seasons
197 198 199
  '("Autumnal Equinox" "Winter Solstice" "Vernal Equinox" "Summer Solstice")
  "List of season changes for the southern hemisphere.")

200 201
(defvar solar-sidereal-time-greenwich-midnight
   nil
202 203
   "Sidereal time at Greenwich at midnight (universal time).")

204
(defvar solar-northern-spring-or-summer-season nil
205
  "Non-nil if northern spring or summer and nil otherwise.
206 207
Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.")

Jim Blandy's avatar
Jim Blandy committed
208 209 210 211 212 213 214 215 216 217 218 219 220 221
(defun solar-setup ()
  "Prompt user for latitude, longitude, and time zone."
  (beep)
  (if (not calendar-longitude)
      (setq calendar-longitude
            (solar-get-number
             "Enter longitude (decimal fraction; + east, - west): ")))
  (if (not calendar-latitude)
      (setq calendar-latitude
            (solar-get-number
             "Enter latitude (decimal fraction; + north, - south): ")))
  (if (not calendar-time-zone)
      (setq calendar-time-zone
            (solar-get-number
222
             "Enter difference from Coordinated Universal Time (in minutes): "))))
Jim Blandy's avatar
Jim Blandy committed
223 224 225 226 227 228 229 230

(defun solar-get-number (prompt)
  "Return a number from the minibuffer, prompting with PROMPT.
Returns nil if nothing was entered."
  (let ((x (read-string prompt "")))
    (if (not (string-equal x ""))
        (string-to-int x))))

231 232 233 234 235 236 237 238 239 240 241 242 243 244
;; 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)))
245

Jim Blandy's avatar
Jim Blandy committed
246 247 248 249 250 251 252 253
(defun solar-xy-to-quadrant (x y)
  "Determines the quadrant of the point X, Y."
  (if (> x 0)
      (if (> y 0) 1 4)
      (if (> y 0) 2 3)))

(defun solar-degrees-to-quadrant (angle)
  "Determines the quadrant of ANGLE."
254
  (1+ (floor (mod angle 360) 90)))
Jim Blandy's avatar
Jim Blandy committed
255 256 257 258 259 260 261 262 263

(defun solar-arctan (x quad)
  "Arctangent of X in quadrant QUAD."
  (let ((deg (radians-to-degrees (atan x))))
    (cond ((equal quad 2)   (+ deg 180))
	  ((equal quad 3)   (+ deg 180))
	  ((equal quad 4)   (+ deg 360))
	  (t                deg))))

264 265
(defun solar-atn2 (x y)
   "Arctan of point X, Y."
266 267
   (if (= x 0)
       (if (> y 0) 90 270)
268
     (solar-arctan (/ y x) (solar-xy-to-quadrant x y))))
269

Jim Blandy's avatar
Jim Blandy committed
270
(defun solar-arccos (x)
271 272 273
     "Arcos of X."
     (let ((y (sqrt (- 1 (* x x)))))
       (solar-atn2 x y)))
Jim Blandy's avatar
Jim Blandy committed
274 275

(defun solar-arcsin (y)
276 277 278 279
     "Arcsin of Y."
     (let ((x (sqrt (- 1 (* y y)))))
       (solar-atn2 x y)
       ))
Jim Blandy's avatar
Jim Blandy committed
280

281 282 283
(defsubst solar-degrees-to-hours (degrees)
  "Convert DEGREES to hours."
  (/ degrees 15.0))
Jim Blandy's avatar
Jim Blandy committed
284

285
(defsubst solar-hours-to-days (hour)
286
  "Convert HOUR to decimal fraction of a day."
287
  (/ hour 24.0))
Jim Blandy's avatar
Jim Blandy committed
288

289 290 291
(defun solar-right-ascension (longitude obliquity)
  "Right ascension of the sun, in hours, given LONGITUDE and OBLIQUITY.
Both arguments are in degrees."
Jim Blandy's avatar
Jim Blandy committed
292 293
  (solar-degrees-to-hours
   (solar-arctan
294
    (* (solar-cosine-degrees obliquity) (solar-tangent-degrees longitude))
Jim Blandy's avatar
Jim Blandy committed
295 296
    (solar-degrees-to-quadrant longitude))))

297 298 299
(defun solar-declination (longitude obliquity)
  "Declination of the sun, in degrees, given LONGITUDE and OBLIQUITY.
Both arguments are in degrees."
Jim Blandy's avatar
Jim Blandy committed
300
  (solar-arcsin
301
   (* (solar-sin-degrees obliquity)
Jim Blandy's avatar
Jim Blandy committed
302 303
      (solar-sin-degrees longitude))))

304 305
(defun solar-sunrise-and-sunset (time latitude longitude height)
  "Sunrise, sunset and length of day.
306 307 308 309 310
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
311
\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
312 313
Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT.

314 315 316 317 318
HEIGHT is the angle the center of the sun has over the horizon for the contact
we are trying to find. For sunrise and sunset, it is usually -0.61 degrees,
accounting for the edge of the sun being on the horizon.

Coordinates are included because this function is called with latitude=1
319
degrees to find out if polar regions have 24 hours of sun or only night."
320 321
  (let* ((rise-time (solar-moment -1 latitude longitude time height))
         (set-time (solar-moment 1 latitude longitude time height))
322 323
         (day-length))
    (if (not (and rise-time set-time))
324 325 326 327 328 329 330
        (if (or (and (> latitude 0)
                     solar-northern-spring-or-summer-season)
                (and (< latitude 0)
                     (not solar-northern-spring-or-summer-season)))
            (setq day-length 24)
	  (setq day-length 0))
      (setq day-length (- set-time rise-time)))
331 332 333
    (list (if rise-time (+ rise-time (/ calendar-time-zone 60.0)) nil)
          (if set-time (+ set-time (/ calendar-time-zone 60.0)) nil)
          day-length)))
334

335
(defun solar-moment (direction latitude longitude time height)
336 337 338 339 340 341 342
  "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
343
\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
344 345
Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT.

346 347 348 349
HEIGHT is the angle the center of the sun has over the horizon for the contact
we are trying to find. For sunrise and sunset, it is usually -0.61 degrees,
accounting for the edge of the sun being on the horizon.

350 351
Uses binary search."
  (let* ((ut (car (cdr time)))
352 353
         (possible t) ; we assume that rise or set are possible
         (utmin (+ ut (* direction 12.0)))
354 355 356 357 358 359
         (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))
360 361
         (hmin (car (cdr
               (solar-horizontal-coordinates (list t0 utmin)
362
                                                latitude longitude t))))
363 364
         (hmax (car (cdr
               (solar-horizontal-coordinates (list t0 utmax)
365 366 367
                                                latitude longitude t)))))
       ; -0.61 degrees is the height of the middle of the sun, when it rises
       ;   or sets.
368 369
     (if (< hmin height)
              (if (> hmax height)
370
                  (while ;(< i 20) ; we perform a simple dichotomy
371
                         ; (> (abs (- hut height)) epsilon)
372 373 374 375
                         (>= (abs (- utmoment utmoment-old))
                             (/ solar-error 60))
                    (setq utmoment-old utmoment)
                    (setq utmoment (/ (+ utmin utmax) 2))
376 377
                    (setq hut (car (cdr
                                    (solar-horizontal-coordinates
378
                                   (list t0 utmoment) latitude longitude t))))
379 380
                    (if (< hut height) (setq utmin utmoment))
                    (if (> hut height) (setq utmax utmoment))
381
                   )
382 383 384
                (setq possible nil)) ; the sun never rises
                (setq possible nil)) ; the sun never sets
     (if (not possible) nil utmoment)))
Jim Blandy's avatar
Jim Blandy committed
385

386
(defun solar-time-string (time time-zone)
387
  "Printable form for decimal fraction TIME in TIME-ZONE.
388 389
Format used is given by `calendar-time-display-form'."
  (let* ((time (round (* 60 time)))
390 391 392
	 (24-hours (/ time 60))
	 (minutes (format "%02d" (% time 60)))
	 (12-hours (format "%d" (1+ (% (+ 24-hours 11) 12))))
Jim Blandy's avatar
Jim Blandy committed
393 394 395 396
	 (am-pm (if (>= 24-hours 12) "pm" "am"))
	 (24-hours (format "%02d" 24-hours)))
    (mapconcat 'eval calendar-time-display-form "")))

397 398 399 400 401 402 403 404

(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)
405
  "Date and Universal Time of local noon at *local date* date.
406 407 408 409

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)
410
         (ut (- 12.0 (/ (calendar-longitude) 15)))
411 412 413
         (te (solar-time-equation date ut)))
    (setq ut (- ut te))
    (if (>= ut 24)
414
        (progn
415 416 417 418
          (setq nd (list (car date) (+ 1 (car (cdr date)))
                         (car (cdr (cdr date)))))
          (setq ut (- ut 24))))
    (if (< ut 0)
419
        (progn
420 421 422 423 424 425 426 427
          (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)))

Jim Blandy's avatar
Jim Blandy committed
428
(defun solar-sunrise-sunset (date)
429 430 431 432 433
  "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))
Pavel Janík's avatar
Pavel Janík committed
434
         ; get the time from the 2000 epoch.
435 436 437 438
         (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
439
          (progn (setq solar-sidereal-time-greenwich-midnight
440
                       (solar-sidereal-time t0))
441
                 (solar-sunrise-and-sunset
442
                  (list t0 (car (cdr exact-local-noon)))
443 444
                  1.0
                  (calendar-longitude) 0)))
445 446 447 448 449 450
         ; 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
451 452 453
            (setq solar-northern-spring-or-summer-season
                  (if (> (car (cdr (cdr equator-rise-set))) 12) t nil))
            (solar-sunrise-and-sunset
454
             (list t0 (car (cdr exact-local-noon)))
455
             (calendar-latitude)
456
             (calendar-longitude) -0.61)))
457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482
         (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."
483
  (/ (- (calendar-absolute-from-gregorian date)
484 485
        (calendar-absolute-from-gregorian '(1 1.5 2000)))
     36525.0))
486

487 488 489 490 491 492
(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
493
\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
494 495 496 497 498 499 500 501 502
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))))
Jim Blandy's avatar
Jim Blandy committed
503

504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539
(defun solar-date-next-longitude (d l)
  "First moment on or after Julian day number D when sun's longitude is a
multiple of L degrees at calendar-location-name with that location's
local time (including any daylight savings rules).

L must be an integer divisor of 360.

Result is in local time expressed astronomical (Julian) day numbers.

The values of calendar-daylight-savings-starts,
calendar-daylight-savings-starts-time, calendar-daylight-savings-ends,
calendar-daylight-savings-ends-time, calendar-daylight-time-offset, and
calendar-time-zone are used to interpret local time."
  (let* ((long)
         (start d)
         (start-long (solar-longitude d))
         (next (mod (* l (1+ (floor (/ start-long l)))) 360))
         (end (+ d (* (/ l 360.0) 400)))
         (end-long (solar-longitude end)))
    (while                 ;; bisection search for nearest minute
        (< 0.00001 (- end start))
      ;; start   <= d    < end
      ;; start-long <= next < end-long when next != 0
      ;; when next = 0, we look for the discontinuity (start-long is near 360
      ;;                and end-long is small (less than l).
      (setq d (/ (+ start end) 2.0))
      (setq long (solar-longitude d))
      (if (or (and (/= next 0) (< long next))
              (and (= next 0) (< l long)))
          (progn
            (setq start d)
            (setq start-long long))
        (setq end d)
        (setq end-long long)))
    (/ (+ start end) 2.0)))

540
(defun solar-horizontal-coordinates
541 542 543 544 545 546
          (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
547
\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
548 549 550 551 552 553 554
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)))
555
         (ah (- (* st 15) (* 15 (car ec)) (* -1 (calendar-longitude))))
556 557
                       ; hour angle (in degrees)
         (de (car (cdr ec)))
558
         (azimuth (solar-atn2 (- (* (solar-cosine-degrees ah)
559 560
                                   (solar-sin-degrees latitude))
                                (* (solar-tangent-degrees de)
561 562
                                   (solar-cosine-degrees latitude)))
                              (solar-sin-degrees ah)))
563
         (height (solar-arcsin
564 565 566 567 568 569 570 571 572 573 574 575 576
                  (+ (* (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
577
\(-0.040945 16), -0.040945 being the number of julian centuries elapsed between
578
Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT."
579
   (let* ((tm (solar-ephemeris-time time))
580 581 582 583 584 585 586 587 588 589 590
          (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)
591
               (* 0.0003032 time time))) ; sun mean longitude
592
         (ml (+ 218.3165
593
                (* 481267.8813 time))) ; moon mean longitude
594 595 596
         (m (+ 357.52910
               (* 35999.05030 time)
               (* -0.0001559 time time)
597
               (* -0.00000048 time time time))) ; sun mean anomaly
598 599
         (i (+ 23.43929111 (* -0.013004167 time)
               (* -0.00000016389 time time)
600
               (* 0.0000005036 time time time))); mean inclination
601 602 603 604 605 606 607
         (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
608 609
                  (solar-sin-degrees (* 3 m))))) ; center equation
         (L (+ l c)) ; total longitude
610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629
         (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)
630
                 (* (solar-tangent-degrees (/ i 2))
631 632 633 634 635
                  (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))
636
                     (* 4 ecc y (solar-sin-degrees m)
637 638 639 640 641 642 643 644
                                (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)))

645 646
(defun solar-longitude (d)
  "Longitude of sun on astronomical (Julian) day number D.
647
Accurary is about 0.0006 degree (about 365.25*24*60*0.0006/360 = 1 minutes).
648 649 650 651 652 653 654

The values of calendar-daylight-savings-starts,
calendar-daylight-savings-starts-time, calendar-daylight-savings-ends,
calendar-daylight-savings-ends-time, calendar-daylight-time-offset, and
calendar-time-zone are used to interpret local time."
  (let* ((a-d (calendar-absolute-from-astro d))
         ;; get Universal Time
655 656 657 658 659
         (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))))
660 661 662 663 664 665 666
         ;; get Ephemeris Time
         (date (+ date (solar-ephemeris-correction
                        (extract-calendar-year
                         (calendar-gregorian-from-absolute
                          (floor
                           (calendar-absolute-from-astro
                            date)))))))
667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737
         (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)))
Jim Blandy's avatar
Jim Blandy committed
738 739

(defun solar-ephemeris-correction (year)
740 741 742 743
  "Ephemeris time minus Universal Time during Gregorian year.
Result is in days.

For the years 1800-1987, the maximum error is 1.9 seconds.
744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798
For the other years, the maximum error is about 30 seconds."
  (cond ((and (<= 1988 year) (< year 2020))
         (/ (+ year -2000 67.0) 60.0 60.0 24.0))
        ((and (<= 1900 year) (< year 1988))
         (let* ((theta (/ (- (calendar-astro-from-absolute
                              (calendar-absolute-from-gregorian
                               (list 7 1 year)))
                             (calendar-astro-from-absolute
                              (calendar-absolute-from-gregorian
                               '(1 1 1900))))
                          36525.0))
                (theta2 (* theta theta))
                (theta3 (* theta2 theta))
                (theta4 (* theta2 theta2))
                (theta5 (* theta3 theta2)))
           (+ -0.00002
              (* 0.000297 theta)
              (* 0.025184 theta2)
              (* -0.181133 theta3)
              (* 0.553040 theta4)
              (* -0.861938 theta5)
              (* 0.677066 theta3 theta3)
              (* -0.212591 theta4 theta3))))
        ((and (<= 1800 year) (< year 1900))
         (let* ((theta (/ (- (calendar-astro-from-absolute
                              (calendar-absolute-from-gregorian
                               (list 7 1 year)))
                             (calendar-astro-from-absolute
                              (calendar-absolute-from-gregorian
                               '(1 1 1900))))
                          36525.0))
                (theta2 (* theta theta))
                (theta3 (* theta2 theta))
                (theta4 (* theta2 theta2))
                (theta5 (* theta3 theta2)))
           (+ -0.000009
              (* 0.003844 theta)
              (* 0.083563 theta2)
              (* 0.865736 theta3)
              (* 4.867575 theta4)
              (* 15.845535 theta5)
              (* 31.332267 theta3 theta3)
              (* 38.291999 theta4 theta3)
              (* 28.316289 theta4 theta4)
              (* 11.636204 theta4 theta5)
              (* 2.043794 theta5 theta5))))
        ((and (<= 1620 year) (< year 1800))
         (let ((x (/ (- year 1600) 10.0)))
           (/ (+ (* 2.19167 x x) (* -40.675 x) 196.58333) 60.0 60.0 24.0)))
        (t (let* ((tmp (- (calendar-astro-from-absolute
                           (calendar-absolute-from-gregorian
                            (list 1 1 year)))
                          2382148))
                  (second (- (/ (* tmp tmp) 41048480.0) 15)))
             (/ second 60.0 60.0 24.0)))))
Jim Blandy's avatar
Jim Blandy committed
799

800 801 802 803 804 805 806 807 808 809 810 811 812
(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
813
       (mod (+ (mod (+ mean-sid-time
814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829
                    (/ (/ (* 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))))

Jim Blandy's avatar
Jim Blandy committed
830 831
;;;###autoload
(defun sunrise-sunset (&optional arg)
832
  "Local time of sunrise and sunset for today.  Accurate to a few seconds.
833
If called with an optional prefix argument, prompt for date.
Jim Blandy's avatar
Jim Blandy committed
834

835 836
If called with an optional double prefix argument, prompt for longitude,
latitude, time zone, and date, and always use standard time.
Jim Blandy's avatar
Jim Blandy committed
837 838 839

This function is suitable for execution in a .emacs file."
 (interactive "p")
840
 (or arg (setq arg 1))
Jim Blandy's avatar
Jim Blandy committed
841 842 843
 (if (and (< arg 16)
          (not (and calendar-latitude calendar-longitude calendar-time-zone)))
     (solar-setup))
Jim Blandy's avatar
Jim Blandy committed
844
 (let* ((calendar-longitude
Jim Blandy's avatar
Jim Blandy committed
845
         (if (< arg 16) calendar-longitude
Jim Blandy's avatar
Jim Blandy committed
846 847 848
           (solar-get-number
            "Enter longitude (decimal fraction; + east, - west): ")))
        (calendar-latitude
Jim Blandy's avatar
Jim Blandy committed
849
         (if (< arg 16) calendar-latitude
Jim Blandy's avatar
Jim Blandy committed
850 851 852
           (solar-get-number
            "Enter latitude (decimal fraction; + north, - south): ")))
        (calendar-time-zone
Jim Blandy's avatar
Jim Blandy committed
853
         (if (< arg 16) calendar-time-zone
Jim Blandy's avatar
Jim Blandy committed
854
           (solar-get-number
855
            "Enter difference from Coordinated Universal Time (in minutes): ")))
Jim Blandy's avatar
Jim Blandy committed
856
        (calendar-location-name
Jim Blandy's avatar
Jim Blandy committed
857 858 859
         (if (< arg 16) calendar-location-name
           (let ((float-output-format "%.1f"))
             (format "%s%s, %s%s"
860 861 862 863 864 865 866 867 868 869 870 871 872
                     (if (numberp calendar-latitude)
                         (abs calendar-latitude)
                       (+ (aref calendar-latitude 0)
                          (/ (aref calendar-latitude 1) 60.0)))
                     (if (numberp calendar-latitude)
                         (if (> calendar-latitude 0) "N" "S")
                       (if (equal (aref calendar-latitude 2) 'north) "N" "S"))
                     (if (numberp calendar-longitude)
                         (abs calendar-longitude)
                       (+ (aref calendar-longitude 0)
                          (/ (aref calendar-longitude 1) 60.0)))
                     (if (numberp calendar-longitude)
                         (if (> calendar-longitude 0) "E" "W")
873
                       (if (equal (aref calendar-longitude 2) 'east)
874
                           "E" "W"))))))
Jim Blandy's avatar
Jim Blandy committed
875
        (calendar-standard-time-zone-name
Jim Blandy's avatar
Jim Blandy committed
876
         (if (< arg 16) calendar-standard-time-zone-name
877
           (cond ((= calendar-time-zone 0) "UTC")
Jim Blandy's avatar
Jim Blandy committed
878
                 ((< calendar-time-zone 0)
879 880
                     (format "UTC%dmin" calendar-time-zone))
                 (t  (format "UTC+%dmin" calendar-time-zone)))))
881 882 883 884
        (calendar-daylight-savings-starts
         (if (< arg 16) calendar-daylight-savings-starts))
        (calendar-daylight-savings-ends
         (if (< arg 16) calendar-daylight-savings-ends))
Jim Blandy's avatar
Jim Blandy committed
885
        (date (if (< arg 4) (calendar-current-date) (calendar-read-date)))
Jim Blandy's avatar
Jim Blandy committed
886
        (date-string (calendar-date-string date t))
887
        (time-string (solar-sunrise-sunset-string date))
Jim Blandy's avatar
Jim Blandy committed
888 889
        (msg (format "%s: %s" date-string time-string))
        (one-window (one-window-p t)))
890
   (if (<= (length msg) (frame-width))
891
       (message "%s" msg)
Jim Blandy's avatar
Jim Blandy committed
892 893
     (with-output-to-temp-buffer "*temp*"
       (princ (concat date-string "\n" time-string)))
894 895
     (message "%s"
	      (substitute-command-keys
Jim Blandy's avatar
Jim Blandy committed
896 897 898 899 900
               (if one-window
                   (if pop-up-windows
                       "Type \\[delete-other-windows] to remove temp window."
                     "Type \\[switch-to-buffer] RET to remove temp window.")
                 "Type \\[switch-to-buffer-other-window] RET to restore old contents of temp window."))))))
901

Jim Blandy's avatar
Jim Blandy committed
902 903
(defun calendar-sunrise-sunset ()
  "Local time of sunrise and sunset for date under cursor.
904
Accurate to a few seconds."
Jim Blandy's avatar
Jim Blandy committed
905 906 907
  (interactive)
  (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
      (solar-setup))
908
  (let ((date (calendar-cursor-to-date t)))
909 910
    (message "%s: %s"
             (calendar-date-string date t t)
911
             (solar-sunrise-sunset-string date))))
Jim Blandy's avatar
Jim Blandy committed
912 913 914

(defun diary-sunrise-sunset ()
  "Local time of sunrise and sunset as a diary entry.
915
Accurate to a few seconds."
Jim Blandy's avatar
Jim Blandy committed
916 917
  (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
      (solar-setup))
918
  (solar-sunrise-sunset-string date))
Jim Blandy's avatar
Jim Blandy committed
919

920 921 922 923 924 925
(defcustom diary-sabbath-candles-minutes 18
  "*Number of minutes before sunset for sabbath candle lighting."
  :group 'diary
  :type 'integer
  :version "21.1")

926
(defun diary-sabbath-candles (&optional mark)
Jim Blandy's avatar
Jim Blandy committed
927
  "Local time of candle lighting diary entry--applies if date is a Friday.
928 929
No diary entry if there is no sunset on that date.

930
An optional parameter MARK specifies a face or single-character string to
931
use when highlighting the day in the calendar."
Jim Blandy's avatar
Jim Blandy committed
932 933 934
  (if (not (and calendar-latitude calendar-longitude calendar-time-zone))
      (solar-setup))
  (if (= (% (calendar-absolute-from-gregorian date) 7) 5);;  Friday
935
      (let* ((sunset (car (cdr (solar-sunrise-sunset date))))
936 937 938 939
             (light (if sunset
                        (cons (- (car sunset)
                                 (/ diary-sabbath-candles-minutes 60.0))
                              (cdr sunset)))))
940
        (if sunset
941 942 943
            (cons mark
		  (format "%s Sabbath candle lighting"
                    (apply 'solar-time-string light)))))))
Jim Blandy's avatar
Jim Blandy committed
944

945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971
; 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)))

972 973 974
(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;
975
K=3, winter solstice.
976 977 978 979 980 981 982 983
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)))))
984 985
         (S (apply '+ (mapcar '(lambda(x)
                                 (* (car x) (solar-cosine-degrees
986
                                             (+ (* (car (cdr (cdr x))) T)
987
                                                  (car (cdr x))))))
988 989
                              solar-seasons-data)))
         (JDE (+ JDE0 (/ (* 0.00001 S) Delta-lambda)))
990
         (correction (+ 102.3 (* 123.5 T) (* 32.5 T T)))
991 992 993 994 995 996 997 998 999 1000 1001
             ; 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)
1002
  "Julian day of mean equinox/solstice K for YEAR.
1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049
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)))))))

1050
;;;###autoload
1051
(defun solar-equinoxes-solstices ()
1052
  "*local* date and time of equinoxes and solstices, if visible in the calendar window.
Jim Blandy's avatar
Jim Blandy committed
1053
Requires floating point."
1054 1055
  (let ((m displayed-month)
        (y displayed-year))
Jim Blandy's avatar
Jim Blandy committed
1056 1057 1058 1059
    (increment-calendar-month m y (cond ((= 1 (% m 3)) -1)
					((= 2 (% m 3))  1)
					(t              0)))
    (let* ((calendar-standard-time-zone-name
1060
            (if calendar-time-zone calendar-standard-time-zone-name "UTC"))
Jim Blandy's avatar
Jim Blandy committed
1061 1062 1063 1064 1065 1066
           (calendar-daylight-savings-starts
            (if calendar-time-zone calendar-daylight-savings-starts))
           (calendar-daylight-savings-ends
            (if calendar-time-zone calendar-daylight-savings-ends))
           (calendar-time-zone (if calendar-time-zone calendar-time-zone 0))
           (k (1- (/ m 3)))
1067
           (d0 (solar-equinoxes/solstices k y))
1068 1069 1070
           (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))
1071 1072 1073 1074
           (d (list (car (car adj))
                    (+ (car (cdr (car adj))  )
                       (/ (car (cdr adj)) 24.0))
                    (car (cdr (cdr (car adj))))))
1075 1076 1077 1078 1079 1080 1081 1082
           ; 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)))
1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095
      (list
       (list (calendar-gregorian-from-absolute (floor abs-day))
             (format "%s %s"
                     (nth k (if (and calendar-latitude
                                     (< (calendar-latitude) 0))
                                solar-s-hemi-seasons
                              solar-n-hemi-seasons))
                     (solar-time-string
                      (* 24 (- abs-day (floor abs-day)))
                      (if (dst-in-effect abs-day)
                          calendar-daylight-time-zone-name
                        calendar-standard-time-zone-name))))))))

Jim Blandy's avatar
Jim Blandy committed
1096 1097 1098

(provide 'solar)