From e7148377c1a45b7b2cfbe9c9fcb553db131a8bde Mon Sep 17 00:00:00 2001 From: Glenn Morris Date: Fri, 14 Mar 2008 07:24:41 +0000 Subject: [PATCH] Reorder so that functions are defined before use. (displayed-month, displayed-year): Move declarations where needed. (solar-get-number): Move definition before use. Use unless. (solar-equatorial-coordinates): Simplify. (solar-sunrise-and-sunset): Use let rather than let*. (solar-longitude, solar-equinoxes-solstices): Use cadr, nth --- lisp/ChangeLog | 7 + lisp/calendar/solar.el | 549 ++++++++++++++++++++--------------------- 2 files changed, 281 insertions(+), 275 deletions(-) diff --git a/lisp/ChangeLog b/lisp/ChangeLog index 75125030b4b..ced8820c45f 100644 --- a/lisp/ChangeLog +++ b/lisp/ChangeLog @@ -1,5 +1,12 @@ 2008-03-14 Glenn Morris + * calendar/solar.el: Reorder so that functions are defined before use. + (displayed-month, displayed-year): Move declarations where needed. + (solar-get-number): Move definition before use. Use unless. + (solar-equatorial-coordinates): Simplify. + (solar-sunrise-and-sunset): Use let rather than let*. + (solar-longitude, solar-equinoxes-solstices): Use cadr, nth + * startup.el (command-line-1): Rename -internal-script back to -scriptload (reverts previous change). diff --git a/lisp/calendar/solar.el b/lisp/calendar/solar.el index d0b02b4b111..78071bb5db8 100644 --- a/lisp/calendar/solar.el +++ b/lisp/calendar/solar.el @@ -54,9 +54,6 @@ ;;; Code: -(defvar displayed-month) -(defvar displayed-year) - (if (fboundp 'atan) (require 'lisp-float-type) (error "Solar calculations impossible since floating point is unavailable")) @@ -206,6 +203,13 @@ Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.") long (- long))))) +(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 ""))) + (unless (string-equal x "") + (string-to-number x)))) + (defun solar-setup () "Prompt for `calendar-longitude', `calendar-latitude', `calendar-time-zone'." (beep) @@ -223,13 +227,6 @@ Needed for polar areas, in order to know whether the day lasts 0 or 24 hours.") "Enter difference from Coordinated Universal Time (in minutes): ") ))) -(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-number x)))) - (defun solar-sin-degrees (x) "Return sin of X degrees." (sin (degrees-to-radians (mod x 360.0)))) @@ -299,36 +296,182 @@ Both arguments are in degrees." (* (solar-sin-degrees obliquity) (solar-sin-degrees longitude)))) -(defun solar-sunrise-and-sunset (time latitude longitude height) - "Sunrise, sunset and length of day. -Parameters are the midday TIME and the LATITUDE, LONGITUDE of the location. +(defun solar-ecliptic-coordinates (time sunrise-flag) + "Return solar longitude, ecliptic inclination, equation of time, nutation. +Values are for TIME in Julian centuries of Ephemeris Time since +January 1st, 2000, at 12 ET. Longitude and inclination are in +degrees, equation of time in hours, and nutation in seconds of longitude. +If SUNRISE-FLAG is non-nil, only calculate longitude and inclination." + (let* ((l (+ 280.46645 + (* 36000.76983 time) + (* 0.0003032 time time))) ; sun mean longitude + (ml (+ 218.3165 + (* 481267.8813 time))) ; moon mean longitude + (m (+ 357.52910 + (* 35999.05030 time) + (* -0.0001559 time time) + (* -0.00000048 time time time))) ; sun mean anomaly + (i (+ 23.43929111 (* -0.013004167 time) + (* -0.00000016389 time time) + (* 0.0000005036 time time time))) ; mean inclination + (c (+ (* (+ 1.914600 + (* -0.004817 time) + (* -0.000014 time time)) + (solar-sin-degrees m)) + (* (+ 0.019993 (* -0.000101 time)) + (solar-sin-degrees (* 2 m))) + (* 0.000290 + (solar-sin-degrees (* 3 m))))) ; center equation + (L (+ l c)) ; total longitude + ;; Longitude of moon's ascending node on the ecliptic. + (omega (+ 125.04 + (* -1934.136 time))) + ;; nut = nutation in longitude, measured in seconds of angle. + (nut (unless sunrise-flag + (+ (* -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)))))) + (ecc (unless sunrise-flag ; eccentricity of earth's orbit + (+ 0.016708617 + (* -0.000042037 time) + (* -0.0000001236 time time)))) + (app (+ L ; apparent longitude of sun + -0.00569 + (* -0.00478 + (solar-sin-degrees omega)))) + (y (unless sunrise-flag + (* (solar-tangent-degrees (/ i 2)) + (solar-tangent-degrees (/ i 2))))) + ;; Equation of time, in hours. + (time-eq (unless sunrise-flag + (/ (* 12 (+ (* y (solar-sin-degrees (* 2 l))) + (* -2 ecc (solar-sin-degrees m)) + (* 4 ecc y (solar-sin-degrees m) + (solar-cosine-degrees (* 2 l))) + (* -0.5 y y (solar-sin-degrees (* 4 l))) + (* -1.25 ecc ecc (solar-sin-degrees (* 2 m))))) + 3.1415926535)))) + (list app i time-eq nut))) + +(defun solar-ephemeris-correction (year) + "Ephemeris time minus Universal Time during Gregorian YEAR. +Result is in days. For the years 1800-1987, the maximum error is +1.9 seconds. For the other years, the maximum error is about 30 seconds." + (cond ((and (<= 1988 year) (< year 2020)) + (/ (+ year -2000 67.0) 60.0 60.0 24.0)) + ((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))))) +(defun solar-ephemeris-time (time) + "Ephemeris Time at moment TIME. TIME is a pair with the first component being the number of Julian centuries elapsed at 0 Universal Time, and the second component being the universal time. For instance, the pair corresponding to November 28, 1995 at 16 UT is \(-0.040945 16), -0.040945 being the number of Julian centuries elapsed between Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. -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. +Result is in Julian centuries of ephemeris time." + (let* ((t0 (car time)) + (ut (cadr time)) + (t1 (+ t0 (/ (/ ut 24.0) 36525))) + (y (+ 2000 (* 100 t1))) + (dt (* 86400 (solar-ephemeris-correction (floor y))))) + (+ t1 (/ (/ dt 86400) 36525)))) -Coordinates are included because this function is called with latitude=1 -degrees to find out if polar regions have 24 hours of sun or only night." - (let* ((rise-time (solar-moment -1 latitude longitude time height)) - (set-time (solar-moment 1 latitude longitude time height)) - (day-length)) - (if (not (and rise-time set-time)) - (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))) - (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))) +(defun solar-equatorial-coordinates (time sunrise-flag) + "Right ascension (in hours) and declination (in degrees) of the sun at TIME. +TIME is a pair with the first component being the number of +Julian centuries elapsed at 0 Universal Time, and the second +component being the universal time. For instance, the pair +corresponding to November 28, 1995 at 16 UT is (-0.040945 16), +-0.040945 being the number of Julian centuries elapsed between +Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. SUNRISE-FLAG is passed +to `solar-ecliptic-coordinates'." + (let ((ec (solar-ecliptic-coordinates (solar-ephemeris-time time) + sunrise-flag))) + (list (solar-right-ascension (car ec) (cadr ec)) + (solar-declination (car ec) (cadr ec))))) + +(defun solar-horizontal-coordinates (time latitude longitude sunrise-flag) + "Azimuth and height of the sun at TIME, LATITUDE, and LONGITUDE. +TIME is a pair with the first component being the number of +Julian centuries elapsed at 0 Universal Time, and the second +component being the universal time. For instance, the pair +corresponding to November 28, 1995 at 16 UT is (-0.040945 16), +-0.040945 being the number of Julian centuries elapsed between +Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. SUNRISE-FLAG +is passed to `solar-ecliptic-coordinates'. Azimuth and +height (between -180 and 180) are both in degrees." + (let* ((ut (cadr time)) + (ec (solar-equatorial-coordinates time sunrise-flag)) + (st (+ solar-sidereal-time-greenwich-midnight + (* ut 1.00273790935))) + ;; Hour angle (in degrees). + (ah (- (* st 15) (* 15 (car ec)) (* -1 (calendar-longitude)))) + (de (cadr ec)) + (azimuth (solar-atn2 (- (* (solar-cosine-degrees ah) + (solar-sin-degrees latitude)) + (* (solar-tangent-degrees de) + (solar-cosine-degrees latitude))) + (solar-sin-degrees ah))) + (height (solar-arcsin + (+ (* (solar-sin-degrees latitude) (solar-sin-degrees de)) + (* (solar-cosine-degrees latitude) + (solar-cosine-degrees de) + (solar-cosine-degrees ah)))))) + (if (> height 180) (setq height (- height 360))) + (list azimuth height))) (defun solar-moment (direction latitude longitude time height) "Sunrise/sunset at location. @@ -377,6 +520,37 @@ Uses binary search." (setq possible nil)) ; the sun never sets (if possible utmoment))) +(defun solar-sunrise-and-sunset (time latitude longitude height) + "Sunrise, sunset and length of day. +Parameters are the midday TIME and the LATITUDE, LONGITUDE of the location. + +TIME is a pair with the first component being the number of Julian centuries +elapsed at 0 Universal Time, and the second component being the universal +time. For instance, the pair corresponding to November 28, 1995 at 16 UT is +\(-0.040945 16), -0.040945 being the number of Julian centuries elapsed between +Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. + +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 +degrees to find out if polar regions have 24 hours of sun or only night." + (let ((rise-time (solar-moment -1 latitude longitude time height)) + (set-time (solar-moment 1 latitude longitude time height)) + day-length) + (if (not (and rise-time set-time)) + (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))) + (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))) + (defun solar-time-string (time time-zone) "Printable form for decimal fraction TIME in TIME-ZONE. Format used is given by `calendar-time-display-form'." @@ -388,13 +562,27 @@ Format used is given by `calendar-time-display-form'." (24-hours (format "%02d" 24-hours))) (mapconcat 'eval calendar-time-display-form ""))) - (defun solar-daylight (time) "Printable form for TIME expressed in hours." (format "%d:%02d" (floor time) (floor (* 60 (- time (floor time)))))) +(defun solar-julian-ut-centuries (date) + "Number of Julian centuries since 1 Jan, 2000 at noon UT for Gregorian DATE." + (/ (- (calendar-absolute-from-gregorian date) + (calendar-absolute-from-gregorian '(1 1.5 2000))) + 36525.0)) + +(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." + (solar-ephemeris-time (list (solar-julian-ut-centuries date) ut))) + +(defun solar-time-equation (date ut) + "Equation of time expressed in hours at Gregorian DATE at Universal time UT." + (nth 2 (solar-ecliptic-coordinates (solar-date-to-et date ut) nil))) + (defun solar-exact-local-noon (date) "Date and Universal Time of local noon at *local date* DATE. The date may be different from the one asked for, but it will be the right @@ -415,6 +603,22 @@ local date. The second component of date should be an integer." (calendar-absolute-from-gregorian nd))) (list nd ut))) +(defun solar-sidereal-time (t0) + "Sidereal time (in hours) in Greenwich at T0 Julian centuries. +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 (nth 3 nut-i)) ; nutation + (i (cadr nut-i))) ; inclination + (mod (+ (mod (+ mean-sid-time + (/ (/ (* nut (solar-cosine-degrees i)) 15) 3600)) 24.0) + 24.0) + 24.0))) + (defun solar-sunrise-sunset (date) "List of *local* times of sunrise, sunset, and daylight on Gregorian DATE. Corresponding value is nil if there is no sunrise/sunset." @@ -467,161 +671,6 @@ Corresponding value is nil if there is no sunrise/sunset." (eval calendar-location-name) (nth 2 l)))) -(defun solar-julian-ut-centuries (date) - "Number of Julian centuries since 1 Jan, 2000 at noon UT for Gregorian DATE." - (/ (- (calendar-absolute-from-gregorian date) - (calendar-absolute-from-gregorian '(1 1.5 2000))) - 36525.0)) - -(defun solar-ephemeris-time (time) - "Ephemeris Time at moment TIME. -TIME is a pair with the first component being the number of Julian centuries -elapsed at 0 Universal Time, and the second component being the universal -time. For instance, the pair corresponding to November 28, 1995 at 16 UT is -\(-0.040945 16), -0.040945 being the number of Julian centuries elapsed between -Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. - -Result is in Julian centuries of ephemeris time." - (let* ((t0 (car time)) - (ut (cadr time)) - (t1 (+ t0 (/ (/ ut 24.0) 36525))) - (y (+ 2000 (* 100 t1))) - (dt (* 86400 (solar-ephemeris-correction (floor y))))) - (+ t1 (/ (/ dt 86400) 36525)))) - -(defun solar-date-next-longitude (d l) - "First time after day D when solar longitude is a multiple of L degrees. -D is a Julian day number. L must be an integer divisor of 360. -The result is for `calendar-location-name', and is in local time -\(including any daylight saving rules) expressed in 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) - long (solar-longitude d)) - (if (or (and (not (zerop next)) (< long next)) - (and (zerop next) (< l long))) - (setq start d - start-long long) - (setq end d - end-long long))) - (/ (+ start end) 2.0))) - -(defun solar-horizontal-coordinates (time latitude longitude sunrise-flag) - "Azimuth and height of the sun at TIME, LATITUDE, and LONGITUDE. -TIME is a pair with the first component being the number of -Julian centuries elapsed at 0 Universal Time, and the second -component being the universal time. For instance, the pair -corresponding to November 28, 1995 at 16 UT is (-0.040945 16), --0.040945 being the number of Julian centuries elapsed between -Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. SUNRISE-FLAG -is passed to `solar-ecliptic-coordinates'. Azimuth and -height (between -180 and 180) are both in degrees." - (let* ((ut (cadr time)) - (ec (solar-equatorial-coordinates time sunrise-flag)) - (st (+ solar-sidereal-time-greenwich-midnight - (* ut 1.00273790935))) - ;; Hour angle (in degrees). - (ah (- (* st 15) (* 15 (car ec)) (* -1 (calendar-longitude)))) - (de (cadr ec)) - (azimuth (solar-atn2 (- (* (solar-cosine-degrees ah) - (solar-sin-degrees latitude)) - (* (solar-tangent-degrees de) - (solar-cosine-degrees latitude))) - (solar-sin-degrees ah))) - (height (solar-arcsin - (+ (* (solar-sin-degrees latitude) (solar-sin-degrees de)) - (* (solar-cosine-degrees latitude) - (solar-cosine-degrees de) - (solar-cosine-degrees ah)))))) - (if (> height 180) (setq height (- height 360))) - (list azimuth height))) - -(defun solar-equatorial-coordinates (time sunrise-flag) - "Right ascension (in hours) and declination (in degrees) of the sun at TIME. -TIME is a pair with the first component being the number of -Julian centuries elapsed at 0 Universal Time, and the second -component being the universal time. For instance, the pair -corresponding to November 28, 1995 at 16 UT is (-0.040945 16), --0.040945 being the number of Julian centuries elapsed between -Jan 1, 2000 at 12 UT and November 28, 1995 at 0 UT. SUNRISE-FLAG is passed -to `solar-ecliptic-coordinates'." - (let* ((tm (solar-ephemeris-time time)) - (ec (solar-ecliptic-coordinates tm sunrise-flag))) - (list (solar-right-ascension (car ec) (car (cdr ec))) - (solar-declination (car ec) (car (cdr ec)))))) - -(defun solar-ecliptic-coordinates (time sunrise-flag) - "Return solar longitude, ecliptic inclination, equation of time, nutation. -Values are for TIME in Julian centuries of Ephemeris Time since -January 1st, 2000, at 12 ET. Longitude and inclination are in -degrees, equation of time in hours, and nutation in seconds of longitude. -If SUNRISE-FLAG is non-nil, only calculate longitude and inclination." - (let* ((l (+ 280.46645 - (* 36000.76983 time) - (* 0.0003032 time time))) ; sun mean longitude - (ml (+ 218.3165 - (* 481267.8813 time))) ; moon mean longitude - (m (+ 357.52910 - (* 35999.05030 time) - (* -0.0001559 time time) - (* -0.00000048 time time time))) ; sun mean anomaly - (i (+ 23.43929111 (* -0.013004167 time) - (* -0.00000016389 time time) - (* 0.0000005036 time time time))) ; mean inclination - (c (+ (* (+ 1.914600 - (* -0.004817 time) - (* -0.000014 time time)) - (solar-sin-degrees m)) - (* (+ 0.019993 (* -0.000101 time)) - (solar-sin-degrees (* 2 m))) - (* 0.000290 - (solar-sin-degrees (* 3 m))))) ; center equation - (L (+ l c)) ; total longitude - ;; Longitude of moon's ascending node on the ecliptic. - (omega (+ 125.04 - (* -1934.136 time))) - ;; nut = nutation in longitude, measured in seconds of angle. - (nut (unless sunrise-flag - (+ (* -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)))))) - (ecc (unless sunrise-flag ; eccentricity of earth's orbit - (+ 0.016708617 - (* -0.000042037 time) - (* -0.0000001236 time time)))) - (app (+ L ; apparent longitude of sun - -0.00569 - (* -0.00478 - (solar-sin-degrees omega)))) - (y (unless sunrise-flag - (* (solar-tangent-degrees (/ i 2)) - (solar-tangent-degrees (/ i 2))))) - ;; Equation of time, in hours. - (time-eq (unless sunrise-flag - (/ (* 12 (+ (* y (solar-sin-degrees (* 2 l))) - (* -2 ecc (solar-sin-degrees m)) - (* 4 ecc y (solar-sin-degrees m) - (solar-cosine-degrees (* 2 l))) - (* -0.5 y y (solar-sin-degrees (* 4 l))) - (* -1.25 ecc ecc (solar-sin-degrees (* 2 m))))) - 3.1415926535)))) - (list app i time-eq nut))) - (defconst solar-data-list '((403406 4.721964 1.621043) (195207 5.937458 62830.348067) @@ -705,8 +754,8 @@ The values of `calendar-daylight-savings-starts', (mapcar (lambda (x) (* (car x) (sin (mod - (+ (car (cdr x)) - (* (car (cdr (cdr x))) U)) + (+ (cadr x) + (* (nth 2 x) U)) (* 2 pi))))) solar-data-list))))) (aberration @@ -716,89 +765,36 @@ The values of `calendar-daylight-savings-starts', (nutation (* -0.0000001 (+ (* 834 (sin A1)) (* 64 (sin A2)))))) (mod (radians-to-degrees (+ longitude aberration nutation)) 360.0))) -(defun solar-ephemeris-correction (year) - "Ephemeris time minus Universal Time during Gregorian YEAR. -Result is in days. For the years 1800-1987, the maximum error is -1.9 seconds. For the other years, the maximum error is about 30 seconds." - (cond ((and (<= 1988 year) (< year 2020)) - (/ (+ year -2000 67.0) 60.0 60.0 24.0)) - ((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))))) - -(defun solar-sidereal-time (t0) - "Sidereal time (in hours) in Greenwich at T0 Julian centuries. -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 (nth 3 nut-i)) ; nutation - (i (cadr nut-i))) ; inclination - (mod (+ (mod (+ mean-sid-time - (/ (/ (* nut (solar-cosine-degrees i)) 15) 3600)) 24.0) - 24.0) - 24.0))) - -(defun solar-time-equation (date ut) - "Equation of time expressed in hours at Gregorian DATE at Universal time UT." - (nth 2 (solar-ecliptic-coordinates (solar-date-to-et date ut) nil))) - -(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." - (solar-ephemeris-time (list (solar-julian-ut-centuries date) ut))) +(defun solar-date-next-longitude (d l) + "First time after day D when solar longitude is a multiple of L degrees. +D is a Julian day number. L must be an integer divisor of 360. +The result is for `calendar-location-name', and is in local time +\(including any daylight saving rules) expressed in 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) + long (solar-longitude d)) + (if (or (and (not (zerop next)) (< long next)) + (and (zerop next) (< l long))) + (setq start d + start-long long) + (setq end d + end-long long))) + (/ (+ start end) 2.0))) ;;;###autoload (defun sunrise-sunset (&optional arg) @@ -1018,6 +1014,9 @@ solstice. These formulas are only to be used between 1000 BC and 3000 AD." (* -0.00823 z z z) (* 0.00032 z z z z))))))) +(defvar displayed-month) ; from generate-calendar +(defvar displayed-year) + ;;;###holiday-autoload (defun solar-equinoxes-solstices () "Local date and time of equinoxes and solstices, if visible in the calendar. @@ -1036,8 +1035,8 @@ Requires floating point." (calendar-time-zone (if calendar-time-zone calendar-time-zone 0)) (k (1- (/ m 3))) (d0 (solar-equinoxes/solstices k y)) - (d1 (list (car d0) (floor (car (cdr d0))) (car (cdr (cdr d0))))) - (h0 (* 24 (- (car (cdr d0)) (floor (car (cdr d0)))))) + (d1 (list (car d0) (floor (cadr d0)) (nth 2 d0))) + (h0 (* 24 (- (cadr d0) (floor (cadr d0))))) (adj (dst-adjust-time d1 h0)) (d (list (caar adj) (+ (car (cdar adj)) -- 2.39.5