(* Helmer Aslaksen, 2011 May 31. Full became a built-in Mathematica symbol in V6, so I changed Full[] to FullMoon[]. DateDifference became a built-in Mathematica symbol in V6, so I changed DateDifference[] to GregorianDateDifference[]. Added Off[Remove::"rmnsm"] and On[Remove::"rmnsm"]. *) (* Calendrica 1.0b1 *) (** ;;;; The Functions (code, comments, and definitions) contained in this ;;;; file (the "Program") were written by Nachum Dershowitz and Edward ;;;; M. Reingold (the "Authors"), who retain all rights to them except ;;;; as granted in the License and subject to the warranty and ;;;; liability limitations below. These Functions are explained in ;;;; the Authors' book, "Calendrical Calculations" (Cambridge ;;;; University Press, 1997), and are subject to an international ;;;; copyright and a Patent Pending on certain functions. ;;;; ;;;; The Authors' public service intent is more liberal than suggested ;;;; by the License below, as are their licensing policies for ;;;; otherwise nonallowed uses such as--without limitation--those in ;;;; commercial, web-site, and large-scale academic contexts. Please ;;;; see the web-site ;;;; ;;;; http://emr.cs.uiuc.edu/~reingold/calendar-book/index.html ;;;; ;;;; for all uses not authorized below; in case there is cause ;;;; for doubt about whether a use you contemplate is authorized, ;;;; please contact the Authors (e-mail: reingold@cs.uiuc.edu). ;;;; For commercial licensing information, contact the authors at ;;;; the Department of Computer Science, 1304 West Springfield ;;;; Avenue, University of Illinois at Urbana-Champaign, Urbana, ;;;; IL 61801-2987, USA. ;;;; ;;;; 1. LICENSE. The Authors grant you a license for personal ;;;; use only. This means that for strictly personal use you may ;;;; copy and use the code, and keep a backup or archival copy ;;;; also. Any other uses, including without limitation, ;;;; allowing the code or its output to be accessed, used, or ;;;; available to others, is not permitted. ;;;; ;;;; 2. WARRANTY. ;;;; ;;;; (a) THE AUTHORS PROVIDE NO WARRANTIES OF ANY KIND, EITHER ;;;; EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITING THE ;;;; GENERALITY OF THE FOREGOING, ANY IMPLIED WARRANTY OF ;;;; MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. ;;;; ;;;; (b) THE AUTHORS SHALL NOT BE LIABLE TO YOU OR ANY THIRD PARTIES ;;;; FOR DAMAGES OF ANY KIND, INCLUDING WITHOUT LIMITATION, ANY LOST ;;;; PROFITS, LOST SAVINGS, OR ANY OTHER INCIDENTAL OR CONSEQUENTIAL ;;;; DAMAGES ARISING OUT OF OR RELATED TO THE USE, INABILITY TO USE, ;;;; OR INACCURACY OF CALCULATIONS, OF THE CODE AND FUNCTIONS ;;;; CONTAINED HEREIN, OR THE BREACH OF ANY EXPRESS OR IMPLIED ;;;; WARRANTY, EVEN IF THE AUTHORS OR PUBLISHER HAVE BEEN ADVISED OF ;;;; THE POSSIBILITY OF THOSE DAMAGES. ;;;; ;;;; (c) The foregoing warranty may give you specific legal ;;;; rights which may vary from state to state in the U.S.A. ;;;; ;;;; 3. LIMITATION OF LICENSEE REMEDIES. You acknowledge and ;;;; agree that your exclusive remedy (in law or in equity), and ;;;; Authors' entire liability with respect to ;;;; the material herein, for any breach of representation or for ;;;; any inaccuracy shall be a refund of the license fee or service ;;;; and handling charge which you paid the Authors, if any. ;;;; ;;;; SOME STATES IN THE U.S.A. DO NOT ALLOW THE EXCLUSION OR ;;;; LIMITATION OF LIABILITY FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES, ;;;; SO THE ABOVE EXCLUSIONS OR LIMITATION MAY NOT APPLY TO YOU. ;;;; 4. DISCLAIMER. Except as expressly set forth above, the Authors: ;;;; ;;;; (a) make no other warranties with respect to the material in ;;;; the Program and expressly disclaim any others; ;;;; ;;;; (b) do not warrant that the material contained in the Program ;;;; will meet your requirements or that their operation shall be ;;;; uninterrupted or error-free; ;;;; ;;;; (c) license this material on an "as is" basis, and the entire ;;;; risk as to the quality, accuracy, and performance of the ;;;; Program is yours, should the code prove defective (except as ;;;; expressly warranted herein). You alone assume the entire cost ;;;; of all necessary corrections. ;;;; ;;;; Sample values for the functions (useful for debugging) are ;;;; given in Appendix C of the book. These sample values are ;;;; not available electronically. ;;;; Last modified 27 May 1997. **) Off[Remove::"rmnsm"] Remove["CalendricaV1`*"] Remove["CalendricaV1`Private`*"] On[Remove::"rmnsm"] BeginPackage["CalendricaV1`"] (*----------------------------------------------------------------------------*) (*=============================================== Public Usage information ===*) (*----------------------------------------------------------------------------*) (* General Information *) Calendrica::usage = "Calendrica is a package for generalized manipulation and conversion of dates through various modern and historical calendars. It is based on the work of Nachum Dershowitz and Edward M. Reingold in their book _Calendrical Calculations_, published by Cambridge University Press, ISBN 0-521-56474-3. This _Mathematica_ translation was performed by Robert C. McNally ." (* Special Names *) Bogus::usage = "Bogus represents nonexistent dates." Calendars::usage = "Calendars[] returns the set of all supported calendars" (* Arithmetical Calendars *) Map[(Evaluate[#[[1]]]::"usage" = StringJoin[ToString[#[[1]]], " represents the ", #[[2]], " calendar. ", ToString[#[[1]]], "[] returns a list of components in ", #[[3]], " calendar dates. ", ToString[#[[1]]], "[fixedDate] returns the ", #[[3]], " calendar date corresponding to the given fixed date."])&, { {Gregorian, "present civil (Gregorian)", "Gregorian"}, {ISO, "ISO commercial", "ISO"}, {Julian, "old civil (Julian)", "Julian"}, {Coptic, "Coptic", "Coptic"}, {Ethiopic, "Ethiopic", "Ethiopic"}, {Islamic, "Islamic (Moslem)", "Islamic"}, {Persian, "Persian (solar)", "Persian"}, {Bahai, "Baha'i", "Baha'i"}, {Hebrew, "Hebrew (Jewish)", "Hebrew"}, {MayanLongCount, "Mayan long count", "long count"}, {MayanHaab, "Mayan haab", "haab"}, {MayanTzolkin, "Mayan tzolkin", "tzolkin"}, {OldHinduSolar, "old (mean) Hindu (Indian)", "old Hindu solar"}, {OldHinduLunar, "old (mean) Hindu (Indian)", "old Hindu lunar"}, {French, "original form of the French Revolutionary", "French Revolutionary"}, {ModifiedFrench, "modified form of the French Revolutionary", "modified French Revolutionary"}, {Chinese, "Chinese", "Chinese"}, {HinduSolar, "new (true) Hindu (Indian) solar", "new Hindu solar"}, {HinduLunar, "new (true) Hindu (Indian) lunar", "new Hindu lunar"} }] (* Names of standard week days *) DayOfWeekNames[] = {Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday} Map[(Evaluate[#]::"usage" = StringJoin[ToString[#], " represents a day of the week. ", ToString[#], "[] returns a corresponding numerical constant."])&, DayOfWeekNames[]] (* Names of standard months *) MonthNames[] = {January, February, March, April, May, June, July, August, September, October, November, December} Map[(Evaluate[#]::"usage" = StringJoin[ToString[#], " represents a month on the Julian and Gregorian calendars. ", ToString[#], "[] returns a corresponding numerical constant."])&, MonthNames[]] (* Conversions from fixed dates to week days *) DayOfWeekFromFixed::usage = "DayOfWeekFromFixed[fixedDate] returns the week day on which the given fixed date falls." (* Conversions between standard week days and standard week day names *) NameFromDayOfWeek::usage = "NameFromDayOfWeek[dayOfWeek] returns the name of the given week day." DayOfWeekNames::usage = "DayOfWeekNames[] returns a list of the names of the standard week days." (* Conversion between standard months and standard month names *) NameFromMonth::usage = "NameFromMonth[month] returns the name of the given standard month." MonthNames::usage = "MonthNames[] returns a list of the names of the standard months." (* Inspectors for components of dates and times *) Map[(Evaluate[#[[1]]]::"usage" = StringJoin[ToString[#[[1]]], "[date] returns the ", #[[2]], " component of the given calendar date. ", ToString[#[[1]]], "[] returns the set of calendars whose dates have a ", #[[2]], " component. "])&, { {CMonth, "month"}, {CDay, "day"}, {CYear, "year"}, {CWeek, "week"}, {CBaktun, "baktun"}, {CKatun, "katun"}, {CTun, "tun"}, {CUinal, "uinal"}, {CKin, "kin"}, {CNumber, "number"}, {CName, "name"}, {CMajor, "major"}, {CCycle, "cycle"}, {CLeap, "leap"}, {CLeapMonth, "leap month"}, {CLeapDay, "leap day"} }] (* Time of Day *) TimeOfDay::usage = "TimeOfDay represents a time of day. TimeOfDay[moment] returns the time of day from the given moment. TimeOfDay[] returns a list of components in a time of day." CHour::usage = "CHour[time] returns the hour component of the given time of day." CMinute::usage = "CMinute[time] returns the minute component of the given time of day." CSecond::usage = "CSecond[time] returns the second component of the given time of day." ToMoment::usage = "ToMoment[timeOfDay] returns the moment corresponding to the given time of day." (* General Calculations *) LeapYearQ::usage = "LeapYearQ[calendar, year] returns True if the given year is a leap year on the given calendar, and False otherwise. LeapYearQ[] returns the set of calendars which support LeapYearQ[calendar, year]" (* Conversion between fixed dates and Julian days *) JDStart::usage = "JDStart[] returns the fixed time of start of the Julian day number." MomentFromJD::usage = "MomentFromJD[jDay] returns the fixed moment corresponding to the given Julian day number." JDFromMoment::usage = "JDFromMoment[moment] returns the Julian day number corresponding to the given fixed moment." FixedFromJD::usage = "FixedFromJD[jDay] returns the fixed date (R.D.) corresponding to the given Julian day number." JDFromFixed::usage = "JDFromFixed[fixedDate] returns the Julian day number corresponding to the given fixed date." (* Gregorian Calendar *) IndependenceDay::usage = "IndependenceDay[gYear] returns the fixed date of American Independence Day in the given Gregorian year." KDayOnOrBefore::usage = "KDayOnOrBefore[fixedDate, k] returns the fixed date of the k-day on or before the given fixed date. A k-day of 0 means Sunday, 1 means Monday, etc." KDayOnOrAfter::usage = "KDayOnOrAfter[fixedDate, k] returns the fixed date of the k-day on or after the given fixed date. A k-day of 0 means Sunday, 1 means Monday, etc." KDayNearest::usage = "KDayNearest[fixedDate, k] returns the fixed date of the k-day nearest the given fixed date. A k-day of 0 means Sunday, 1 means Monday, etc." KDayAfter::usage = "KDayAfter[fixedDate, k] returns the fixed date of the k-day after the given fixed date. A k-day of 0 means Sunday, 1 means Monday, etc." KDayBefore::usage = "KDayBefore[fixedDate, k] returns the fixed date of the k-day before the given fixed date. A k-day of 0 means Sunday, 1 means Monday, etc." NthKDay::usage = "NthKDay[n, k, date] returns the fixed date of the n-th k-day after the given calendar date. If n > 0 the n-th k-day on or after the given date is returned. If n < 0, the n-th k-day on or before the given date is returned. A k-day of 0 means Sunday, 1 means Monday, etc." OFirst::usage = "OFirst[] returns an index for selecting a k-day." OLast::usage = "OLast[] returns an index for selecting a k-day." LaborDay::usage = "LaborDay[gYear] returns the fixed date of America Labor Day in the given Gregorian year-- the first Monday in September." MemorialDay::usage = "MemorialDay[gYear] returns the fixed date of American Memorial Day in the given Gregorian year-- the last Monday in May." ElectionDay::usage = "ElectionDay[gYear] returns the fixed date of American Election Day in the given Gregorian year-- the Tuesday after the first Monday in November." DaylightSavingsStart::usage = "DaylightSavingsStart[gYear] returns the fixed date of the start of American daylight savings time in the given Gregorian year-- the first Sunday in April." DaylightSavingsEnd::usage = "DaylightSavingsEnd[gYear] returns the fixed date of the end of American daylight savings time in the given Gregorian year-- the last Sunday in October." Christmas::usage = "Christmas[gYear] returns the fixed date of Christmas in the given Gregorian year." Advent::usage = "Advent[gYear] returns the fixed date of Advent in the given Gregorian year." Epiphany::usage = "Epiphany[gYear] returns the fixed date of Epiphany in the given Gregorian year." (* Julian Calendar *) NicaeanRuleEaster::usage = "NicaeanRuleEaster[jYear] returns the fixed date of Easter in the given positive Julian year, according to the rule of the Council of Nicaea." Easter::usage = "Easter[gYear] returns the fixed date of Easter in the given Gregorian year." JulianInGregorian::usage = "JulianInGregorian[jMonth, jDay, gYear] returns a list of the fixed dates of the given Julian month and day that occur in the given Gregorian year." EasternOrthodoxChristmas::usage = "EasternOrthodoxChristmas[gYear] returns a list of zero or one fixed dates of Eastern Orthodox Christmas in the given Gregorian year." Pentecost::usage = "Pentecost[gYear] returns the fixed date of Pentecost in the given Gregorian year." (* Coptic Calendar *) CopticInGregorian::usage = "CopticInGregorian[cMonth, cDay, gYear] returns a list of the fixed dates of the given Coptic month and day that occur in the given Gregorian year." CopticChristmas::usage = "CopticChristmas[gYear] returns a list of zero or one fixed dates of Coptic Christmas in the given Gregorian year." (* Islamic Calendar *) IslamicInGregorian::usage = "IslamicInGregorian[iMonth, iDay, gYear] returns a list of the fixed dates of the given Islamic month and day that occur in the given Gregorian year." MawlidAnNabi::usage = "MawlidAnNabi[gYear] returns a list of fixed dates of Mawlid-an-Nabi occurring in the given Gregorian year." (* Baha'i Calendar *) BahaiNewYear::usage = "BahaiNewYear[gYear] returns the fixed date of the Baha'i New Year in the given Gregorian year." FeastOfRidvan::usage = "FeastOfRidvan[gYear] returns the fixed date of the Baha'i Feast of Ridvan in the given Gregorian year." (* Persian Calendar *) NawRuz::usage = "NawRuz[gYear] returns the fixed date of the Persian New Year (Naw-Ruz) in the given Gregorian year." (* Hebrew Calendar *) YomKippur::usage = "YomKippur[gYear] returns the fixed date of Yom Kippur occurring in the given Gregorian year." Passover::usage = "Passover[gYear] returns the fixed date of Passover occurring in the given Gregorian year." Omer::usage = "Omer[fixedDate] returns a list of the form {weeks, days} containing the number of elapsed weeks and days in the omer at the given date. Returns Bogus if that date does not fall during the omer." Purim::usage = "Purim[gYear] returns the fixed date of Purim occurring in the given Gregorian year." TaAnithEsther::usage = "TaAnithEsther[gYear] returns the fixed date of Ta'anith Esther occurring in the given Gregorian year." TishaBAv::usage = "TishaBAv[gYear] returns the fixed date of Tisha B'Av occurring in the given Gregorian year." BirkathHaHama::usage = "BirkathHaHama[gYear] returns a list of the fixed dates of Birkath HaHama occurring in the given Gregorian year, if it occurs." ShEla::usage = "ShEla[gYear] returns the fixed date of Sh'ela occurring in the given Gregorian year." YomHaZikaron::usage = "YomHaZikarom[gYear] returns the fixed date of Yom HaZikaron occurring in the given Gregorian year." HebrewBirthday::usage = "HebrewBirthday[hBirthDate, hYear] returns the fixed date of the anniversary of the given Hebrew calendar date occurring in the given Hebrew year." Yahrzeit::usage = "Yahrzeit[hDeathDate, hYear] returns the fixed date of the anniversary of the given Hebrew death-date occurring in the given Hebrew year." (* Mayan Calendars *) MayanHaabDifference::usage = "MayanHaabDifference[hDate1, hDate2] returns the number of days from the first given Mayan haab calendar date hDate1 to the next occurrence of the second given Mayan haab calendar date hDate2." MayanHaabOnOrBefore::usage = "MayanHaabOnOrBefore[hDate, date] returns the fixed date of the latest date on or before the given fixed date that is the given Mayan haab calendar date hDate." MayanTzolkinDifference::usage = "MayanTzolkinDifference[tDate1, tDate2] returns the number of days from the first given Mayan Tzolkin calendar date tDate1 to the next occurrence of the second given Mayan Tzolkin calendar date tDate2." MayanTzolkinOnOrBefore::usage = "MayanTzolkinOnOrBefore[tDate, fixedDate] returns the fixed date of the latest date on or before the given fixed date that is Mayan tzolkin date tDate." MayanHaabTzolkinOnOrBefore::usage = "MayanHaabTzolkinOnOrBefore[hDate, tDate, fixedDate] returns the fixed date of the latest date on or before the given fixed date that is Mayan haab date hDate and tzolkin date tDate. Returns Bogus if such a haab-tzolkin combination is impossible." (* Old Hindu Calendars *) JovianYear::usage = "JovianYear[fixedDate] returns the year of the Jupiter cycle at the given fixed date." HinduDayCount::usage = "HinduDayCount[date] returns the number of elapsed days (Ahargana) to date since the Hindu epoch (K.Y.) from the given fixed date." (* Modern Hindu Calendars *) SunriseAtUjjain::usage = "SunriseAtUjjain[kyTime] returns the hindu-moment of astronomical sunrise at Ujjain (latitude 1389 minutes) at the give kyTime, in whole days since the Hindu epoch." LunarMansion::usage = "LunarMansion[date] returns the Hindu lunar mansion (nakshatra) at sunrise on the given fixed date." Samkranti::usage = "Samkranti[gYear, m] returns the fixed moment of the start of the m-th solar month of the Hindu year beginning in the given Gregorian year." MeshaSamkranti::usage = "MeshaSamkranti[gYear] returns the fixed moment of Mesha samkranti (Vernal equinox) in the given Gregorian year." HinduLunarNewYear::usage = "HinduLunarNewYear[gYear] returns the fixed date of the Hindu luni-solar new year in the given Gregorian year." Karana::usage = "Karana[n] returns the number (0-10) of the name of the n-th (1-60) Hindu karana." Yoga::usage = "Yoga[kyTime] returns the Hindu yoga at ky-time days since the Hindu epoch." SacredWednesdaysInGregorian::usage = "SacredWednesdaysInGregorian[gYear] returns a list of Wednesdays in the given Gregorian year that are day 8 of Hindu lunar months." SacredWednesdays::usage = "SacredWednesdays[start, end] returns a list of Wednesdays between the given fixed dates (inclusive) that are day 8 of Hindu lunar months." (* Time and Astronomy *) Sunrise::usage = "Sunrise[fixedDate, latitude, longitude] returns the local time (fraction of day) of sunrise at the given latitude (+ is N) and longitude (+ is E) (in nonpolar regions) for the given fixed date." Sunset::usage = "Sunset[fixedDate, latitude, longitude] returns the local time (fraction of day) of sunset at the given latitude (+ is N) and longitude (+ is E) (in nonpolar regions) for the given fixed date." LocalFromUniversal::usage = "LocalFromUniversal[uTime, zone] returns the local time (fraction of day) in the given time zone (in minutes of time after Greenwich) for the given universal time." UniversalFromLocal::usage = "UniversalFromLocal[lTime, zone] returns the universal time (fraction of day) for the given local time in the given time zone (in minutes of time after Greenwich)." LocalFromStandard::usage = "LocalFromStandard[sTime, offset] returns the local time (moment) from the given standard time (moment) at distance offset (in minutes) from the local time zone." StandardFromLocal::usage = "StandardFromLocal[lTime, offset] returns the standard time (moment) from the given local time (moment) at distance offset (in minutes) from the local time zone." LocalFromApparent::usage = "LocalFromApparent[moment] returns the local time from the given sundial time." ApparentFromLocal::usage = "ApparentFromLocal[moment] returns the sundial time at the given local time." SolarLongitude::usage = "SolarLongitude[jd] returns the longitude of the sun on the given astronomical (julian) day." DateNextSolarLongitude::usage = "DateNextSolarLongitude[jd, l] returns the Julian day number of the first date at or after julian day number jd (in Greenwich) when the solar longitude will be a multiple of l degrees; l must be a proper divisor of 360." LunarLongitude::usage = "LunarLongitude[uTime] returns the Longitude of the moon (in degrees) at uTime (Universal time)." NewMoonAtOrAfter::usage = "NewMoonAtOrAfter[jd] returns the astronomical (julian) day number of the first new moon at or after astronomical (julian) day number jd (in Greenwich). The fractional part is the time of day." NewMoonBefore::usage = "NewMoonBefore[jd] returns the astronomical (julian) day number of the last new moon before astronomical (julian) day number jd (in Greenwich). The fractional part is the time of day." New::usage = "New represents the new moon phase. New[] returns a corresponding excess of lunar longitude over solar longitude at the new moon." FirstQuarter::usage = "FirstQuarter represents the first quarter moon phase. FirstQuarter[] returns a corresponding excess of lunar longitude over solar longitude at the first quarter." FullMoon::usage = "FullMoon represents the new moon phase. FullMoon[] returns a corresponding excess of lunar longitude over solar longitude at the full moon." LastQuarter::usage = "LastQuarter represents the last quarter moon phase. LastQuarter[] returns a corresponding excess of lunar longitude over solar longitude at the last quarter." LunarSolarAngle::usage = "LunarSolarAngle[jd] returns the lunar phase, as an angle in degrees, at astronomical (julian) day number jd. An angle of 0 means a new moon, 90 degrees means the first quarter, 180 means a full moon, and 270 degrees means the last quarter." LunarPhaseAtOrBefore::usage = "LunarPhaseAtOrBefore[phase, jd] returns the astronomical (julian) day number of last time moon was in phase phase (in degrees) at or before astronomical (julian) day number jd (in Greenwich). The fractional part is the time of day." NewMoonAtOrBefore::usage = "NewMoonAtOrBefore[jd] returns the astronomical (julian) day number of last new moon at or before astronomical (julian) day number jd (in Greenwich). The fractional part is the time of day." FullMoonAtOrBefore::usage = "FullMoonAtOrBefore[jd] returns the astronomical (julian) day number of last full moon at or before astronomical (julian) day number jd (in Greenwich). The fractional part is the time of day." FirstQuarterMoonAtOrBefore::usage = "FirstQuarterMoonAtOrBefore[jd] returns the astronomical (julian) day number of last first-quarter moon at or before astronomical (julian) day number jd (in Greenwich). The fractional part is the time of day." LastQuarterMoonAtOrBefore::usage = "LastQuarterMoonAtOrBefore[jd] returns the astronomical (julian) day number of last last-quarter moon at or before astronomical (julian) day number jd (in Greenwich). The fractional part is the time of day." SiderealFromJD::usage = "SiderealFromJD[jd] returns the siderial time from the given julian day number." (*----------------------------------------------------------------------------*) (* Conversion between fixed dates and calendar dates *) ToFixed::usage = "ToFixed[date] returns the fixed (R.D.) date corresponding to the given calendar date." (* Conversion between calendars *) ConvertDateTo::usage = "ConvertDateTo[date, calendar] converts the given calendar date to the specified calendar." (* Other manipulations on dates *) GregorianDateDifference::usage = "GregorianDateDifference[date1, date2] returns the number of days from calendar date 1 until calendar date 2, performing calendar conversion as necessary." DayNumber::usage = "DayNumber[date] returns the number of the day in the year of the given calendar date. DayNumber[] returns the set of calendars supported by DayNumber[date]." DaysRemaining::usage = "DaysRemaining[date] returns the number of days remaining in the year after the given calendar date. DaysRemaining[] returns the set of calendars supported by DaysRemaining[date]." (* Testing *) CalendricaTest::usage = "CalendricaTest[] runs a series of tests on the date conversion algorithms, comparing the results against test data transcribed from Appendix C of _Calendrical Calculations_ and elsewhere. As it runs, it prints a list of tests performed and either \"OK\" if the algorithms appear to be working correctly, or a lists containing any discrepancies found." (*----------------------------------------------------------------------------*) (*==================================================== Private Definitions ===*) (*----------------------------------------------------------------------------*) Begin["`Private`"] (** ;;;; Section: Basic Code **) (* ADDED *) Calendars[] = {Gregorian, ISO, Julian, Coptic, Ethiopic, Islamic, Persian, Bahai, Hebrew, MayanLongCount, MayanHaab, MayanTzolkin, OldHinduSolar, OldHinduLunar, French, ModifiedFrench, Chinese, HinduSolar, HinduLunar} LeapYearQ[] = {Gregorian, Julian, Coptic, Islamic, Persian, Hebrew, OldHinduLunar, ModifiedFrench} (** (defconstant false ;; TYPE boolean ;; Constant representing false. nil) **) (* False is a pre-defined Mathematica constant *) (** (defconstant bogus ;; TYPE string ;; Used to denote nonexistent dates. "bogus") **) Bogus (** (defun quotient (m n) ;; TYPE (real non-zero-real) -> integer ;; Whole part of m/n. (floor m n)) **) (* Quotient is a pre-defined Mathematica function defined exactly the same way as it is here *) (** (defun adjusted-mod (m n) ;; TYPE (integer positive-integer) -> positive-integer ;; Positive remainder of m/n with n instead of 0. (1+ (mod (1- m) n))) **) AdjustedMod[m_, n_] := Mod[m - 1, n] + 1 (** (defmacro sum (expression index initial condition) ;; TYPE ((integer->real) * integer (integer->boolean)) ;; TYPE -> real ;; Sum expression for index = initial and successive ;; integers, as long as condition holds. (let* ((temp (gensym))) `(do ((,temp 0 (+ ,temp ,expression)) (,index ,initial (1+ ,index))) ((not ,condition) ,temp)))) **) (* Since Sum is already a Mathematica function, we rename the sum macro MSum.*) MSum[expr_, initial_, condition_] := Module[{total, i}, For[total = 0; i = initial, condition[i], i++, total += expr[i]]; total ] (** (defmacro binary-search (l lo h hi x test end) ;; TYPE ( * real * real * (real->boolean) ;; TYPE ((real real)->boolean)) -> real ;; Bisection search for x in lo..hi such that end holds. ;; test determines when to go left. (let* ((left (gensym))) `(do* ((,x false (/ (+ ,h ,l) 2)) (,left false ,test) (,l ,lo (if ,left ,l ,x)) (,h ,hi (if ,left ,x ,h))) (,end (/ (+ ,h ,l) 2))))) **) BinarySearch[lo_, hi_, test_, end_] := Module[{l, h, x}, For[ l = lo; h = hi; x = (h + l) / 2, (* establish the low and high bounds, and the point between them *) !end[l, h], (* continue until the end condition is true *) Identity, (* no action in increment phase *) If[test[x], h = x, l = x]; (* body determines whether to go left or right, and resets the upper or lower bounds accordingly *) x = (h + l) / 2 (* then it resets the test point to lie between the bounds *) ]; x (* return the last test point visited *) ] (** (defmacro sigma (list body) ;; TYPE (list-of-pairs (list-of-reals->real)) ;; TYPE -> real ;; list is of the form ((i1 l1)..(in ln)). ;; Sum of body for indices i1..in ;; running simultaneously thru lists l1..ln. `(apply '+ (mapcar (function (lambda ,(mapcar 'car list) ,body)) ,@(mapcar 'cadr list)))) **) Sigma[list_, body_] := Apply[Plus, Apply[Function[Evaluate[Map[First, list]], body], Map[Part[#, 2]&, list] ]] (** (defun poly (x a) ;; TYPE (real list-of-reals) -> real ;; Sum powers of x with coefficients (from order 0 up) ;; in list a. (if (equal a nil) 0 (+ (first a) ( * x (poly x (cdr a)))))) **) Poly[x_, a_] := a[[1]] + Sum[a[[i]] x^(i - 1), {i, 2, Length[a]}] (** (defconstant sunday ;; TYPE day-of-week ;; Residue class for Sunday. 0) (defconstant monday ;; TYPE day-of-week ;; Residue class for Monday. (1+ sunday)) (defconstant tuesday ;; TYPE day-of-week ;; Residue class for Tuesday. (+ sunday 2)) (defconstant wednesday ;; TYPE day-of-week ;; Residue class for Wednesday. (+ sunday 3)) (defconstant thursday ;; TYPE day-of-week ;; Residue class for Thursday. (+ sunday 4)) (defconstant friday ;; TYPE day-of-week ;; Residue class for Friday. (+ sunday 5)) (defconstant saturday ;; TYPE day-of-week ;; Residue class for Saturday. (+ sunday 6)) **) Sunday[] = 0 Monday[] = Sunday[] + 1 Tuesday[] = Sunday[] + 2 Wednesday[] = Sunday[] + 3 Thursday[] = Sunday[] + 4 Friday[] = Sunday[] + 5 Saturday[] = Sunday[] + 6 (* ADDED *) NameFromNumber[number_Integer, nameList_List] := nameList[[ AdjustedMod[number, Length[nameList]] ]] (** (defun day-of-week-from-fixed (date) ;; TYPE fixed-date -> day-of-week ;; The residue class of the day of the week of date. (mod date 7)) **) DayOfWeekFromFixed[date_Integer] := Mod[date, 7] (* ADDED *) NameFromDayOfWeek[dayOfWeek_Integer] := NameFromNumber[dayOfWeek + 1, DayOfWeekNames[]] (** (defun standard-month (date) ;; TYPE standard-date -> standard-month ;; Month field of date = (month day year). (first date)) (defun standard-day (date) ;; TYPE standard-date -> standard-day ;; Day field of date = (month day year). (second date)) (defun standard-year (date) ;; TYPE standard-date -> standard-year ;; Year field of date = (month day year). (third date)) **) CMonth[ date_Gregorian | date_Julian | date_Islamic | date_Hebrew | date_French | date_ModifiedFrench | date_Coptic | date_Ethiopic | date_Persian | date_OldHinduSolar | date_HinduSolar] := date[[1]] CDay[ date_Gregorian | date_Julian | date_Islamic | date_Hebrew | date_French | date_ModifiedFrench | date_Coptic | date_Ethiopic | date_Persian | date_OldHinduSolar | date_HinduSolar] := date[[2]] CYear[ date_Gregorian | date_Julian | date_Islamic | date_Hebrew | date_French | date_ModifiedFrench | date_Coptic | date_Ethiopic | date_Persian | date_OldHinduSolar | date_HinduSolar] := date[[3]] (** (defun iso-week (date) ;; TYPE iso-date -> iso-week (first date)) (defun iso-day (date) ;; TYPE iso-date -> day-of-week (second date)) (defun iso-year (date) ;; TYPE iso-date -> iso-year (third date)) **) CWeek[date_ISO] := date[[1]] CDay[date_ISO] := date[[2]] CYear[date_ISO] := date[[3]] (** (defun mayan-baktun (date) ;; TYPE mayan-long-count-date -> mayan-baktun (first date)) (defun mayan-katun (date) ;; TYPE mayan-long-count-date -> mayan-katun (second date)) (defun mayan-tun (date) ;; TYPE mayan-long-count-date -> mayan-tun (third date)) (defun mayan-uinal (date) ;; TYPE mayan-long-count-date -> mayan-uinal (fourth date)) (defun mayan-kin (date) ;; TYPE mayan-long-count-date -> mayan-kin (fifth date)) **) CBaktun[date_MayanLongCount] := date[[1]] CKatun[date_MayanLongCount] := date[[2]] CTun[date_MayanLongCount] := date[[3]] CUinal[date_MayanLongCount] := date[[4]] CKin[date_MayanLongCount] := date[[5]] (** (defun mayan-haab-day (date) ;; TYPE mayan-haab-date -> mayan-haab-day (first date)) (defun mayan-haab-month (date) ;; TYPE mayan-haab-date -> mayan-haab-month (second date)) **) CDay[date_MayanHaab] := date[[1]] CMonth[date_MayanHaab] := date[[2]] (** (defun mayan-tzolkin-number (date) ;; TYPE mayan-tzolkin-date -> mayan-tzolkin-number (first date)) (defun mayan-tzolkin-name (date) ;; TYPE mayan-tzolkin-date -> mayan-tzolkin-name (second date)) **) CNumber[date_MayanTzolkin] := date[[1]] CName[date_MayanTzolkin] := date[[2]] (** (defun bahai-major (date) ;; TYPE bahai-date -> bahai-major (first date)) (defun bahai-cycle (date) ;; TYPE bahai-date -> bahai-cycle (second date)) (defun bahai-year (date) ;; TYPE bahai-date -> bahai-year (third date)) (defun bahai-month (date) ;; TYPE bahai-date -> bahai-month (fourth date)) (defun bahai-day (date) ;; TYPE bahai-date -> bahai-day (fifth date)) **) CMajor[date_Bahai] := date[[1]] CCycle[date_Bahai] := date[[2]] CYear[date_Bahai] := date[[3]] CMonth[date_Bahai] := date[[4]] CDay[date_Bahai] := date[[5]] (** (defun chinese-cycle (date) ;; TYPE chinese-date -> chinese-cycle (first date)) (defun chinese-year (date) ;; TYPE chinese-date -> chinese-year (second date)) (defun chinese-month (date) ;; TYPE chinese-date -> chinese-month (third date)) (defun chinese-leap (date) ;; TYPE chinese-date -> chinese-leap (fourth date)) (defun chinese-day (date) ;; TYPE chinese-date -> chinese-day (fifth date)) **) CCycle[date_Chinese] := date[[1]] CYear[date_Chinese] := date[[2]] CMonth[date_Chinese] := date[[3]] CLeap[date_Chinese] := date[[4]] CDay[date_Chinese] := date[[5]] (** (defun old-hindu-lunar-month (date) ;; TYPE old-hindu-lunar-date -> old-hindu-lunar-month (first date)) (defun old-hindu-lunar-leap (date) ;; TYPE old-hindu-lunar-date -> old-hindu-lunar-leap (second date)) (defun old-hindu-lunar-day (date) ;; TYPE old-hindu-lunar-date -> old-hindu-lunar-day (third date)) (defun old-hindu-lunar-year (date) ;; TYPE old-hindu-lunar-date -> old-hindu-lunar-year (fourth date)) **) CMonth[date_OldHinduLunar] := date[[1]] CLeap[date_OldHinduLunar] := date[[2]] CDay[date_OldHinduLunar] := date[[3]] CYear[date_OldHinduLunar] := date[[4]] (** (defun hindu-lunar-month (date) ;; TYPE hindu-lunar-date -> hindu-lunar-month (first date)) (defun hindu-lunar-leap-month (date) ;; TYPE hindu-lunar-date -> hindu-lunar-leap-month (second date)) (defun hindu-lunar-day (date) ;; TYPE hindu-lunar-date -> hindu-lunar-day (third date)) (defun hindu-lunar-leap-day (date) ;; TYPE hindu-lunar-date -> hindu-lunar-leap-day (fourth date)) (defun hindu-lunar-year (date) ;; TYPE hindu-lunar-date -> hindu-lunar-year (fifth date)) **) CMonth[date_HinduLunar] := date[[1]] CLeapMonth[date_HinduLunar] := date[[2]] CDay[date_HinduLunar] := date[[3]] CLeapDay[date_HinduLunar] := date[[4]] CYear[date_HinduLunar] := date[[5]] (* ADDED *) CHour[time_TimeOfDay] := time[[1]] CMinute[time_TimeOfDay] := time[[2]] CSecond[time_TimeOfDay] := time[[3]] (* ADDED *) CMonth[] = {Gregorian, Julian, Islamic, Hebrew, French, ModifiedFrench, Coptic, Ethiopic, Persian, OldHinduSolar, HinduSolar, Bahai, Chinese, OldHinduLunar, HinduLunar} CDay[] = {Gregorian, Julian, Islamic, Hebrew, French, ModifiedFrench, Coptic, Ethiopic, Persian, OldHinduSolar, HinduSolar, ISO, MayanHaab, Bahai, Chinese, OldHinduLunar, HinduLunar} CYear[] = {Gregorian, Julian, Islamic, Hebrew, French, ModifiedFrench, Coptic, Ethiopic, Persian, OldHinduSolar, HinduSolar, ISO, MayanHaab, Bahai, Chinese, OldHinduLunar, HinduLunar} CWeek[] = {ISO} CBaktun[] = {MayanLongCount} CKatun[] = {MayanLongCount} CTun[] = {MayanLongCount} CUinal[] = {MayanLongCount} CKin[] = {MayanLongCount} CNumber[] = {MayanTzolkin} CName[] = {MayanTzolkin} CMajor[] = {Bahai} CCycle[] = {Bahai, Chinese} CLeap[] = {Chinese, OldHinduLunar} CLeapMonth[] = {HinduLunar} CLeapDay[] = {HinduLunar} Gregorian[] = {CMonth, CDay, CYear} Julian[] = {CMonth, CDay, CYear} Islamic[] = {CMonth, CDay, CYear} Hebrew[] = {CMonth, CDay, CYear} French[] = {CMonth, CDay, CYear} ModifiedFrench[] = {CMonth, CDay, CYear} Coptic[] = {CMonth, CDay, CYear} Ethiopic[] = {CMonth, CDay, CYear} Persian[] = {CMonth, CDay, CYear} OldHinduSolar[] = {CMonth, CDay, CYear} HinduSolar[] = {CMonth, CDay, CYear} ISO[] = {CWeek, CDay, CYear} MayanLongCount[] = {CBaktun, CKatun, CTun, CUinal, CKin} MayanHaab[] = {CDay, CMonth} MayanTzolkin[] = {CNumber, CName} Chinese[] = {CCycle, CYear, CMonth, CLeap, CDay} Bahai[] = {CMajor, CCycle, CYear, CMonth, CDay} OldHinduLunar[] = {CMonth, CLeap, CDay, CYear} HinduLunar[] = {CMonth, CLeapMonth, CDay, CLeapDay, CYear} TimeOfDay[] = {CHour, CMinute, CSecond} (** (defun gregorian-date (month day year) ;; TYPE (gregorian-month gregorian-day gregorian-year) ;; TYPE -> gregorian-date (list month day year)) (defun julian-date (month day year) ;; TYPE (julian-month julian-day julian-year) ;; TYPE -> julian-date (list month day year)) (defun islamic-date (month day year) ;; TYPE (islamic-month islamic-day islamic-year) ;; TYPE -> islamic-date (list month day year)) (defun hebrew-date (month day year) ;; TYPE (hebrew-month hebrew-day hebrew-year) -> hebrew-date (list month day year)) (defun french-date (month day year) ;; TYPE (french-month french-day french-year) -> french-date (list month day year)) (defun coptic-date (month day year) ;; TYPE (coptic-month coptic-day coptic-year) -> coptic-date (list month day year)) (defun ethiopic-date (month day year) ;; TYPE (ethiopic-month ethiopic-day ethiopic-year) ;; TYPE -> ethiopic-date (list month day year)) (defun persian-date (month day year) ;; TYPE (persian-month persian-day persian-year) ;; TYPE -> persian-date (list month day year)) (defun hindu-solar-date (month day year) ;; TYPE (hindu-solar-month hindu-solar-day hindu-solar-year) ;; TYPE -> hindu-solar-date (list month day year)) (defun time-of-day (hour minute second) ;; TYPE (hour minute second) -> moment (list hour minute second)) (defun iso-date (week day year) ;; TYPE (iso-week iso-day iso-year) -> iso-date (list week day year)) (defun mayan-long-count-date (baktun katun tun uinal kin) ;; TYPE (mayan-baktun mayan-katun mayan-tun mayan-uinal ;; TYPE mayan-kin) -> mayan-long-count-date (list baktun katun tun uinal kin)) (defun mayan-haab-date (day month) ;; TYPE (mayan-haab-day mayan-haab-month) -> mayan-haab-date (list day month)) (defun mayan-tzolkin-date (number name) ;; TYPE (mayan-tzolkin-number mayan-tzolkin-name) ;; TYPE -> mayan-tzolkin-date (list number name)) (defun chinese-date (cycle year month leap day) ;; TYPE (chinese-cycle chinese-year chinese-month ;; TYPE chinese-leap chinese-day) -> chinese-date (list cycle year month leap day)) (defun bahai-date (major cycle year month day) ;; TYPE (bahai-major bahai-cycle bahai-year ;; TYPE bahai-month bahai-day) -> bahai-date (list major cycle year month day)) (defun old-hindu-lunar-date (month leap day year) ;; TYPE (old-hindu-lunar-month old-hindu-lunar-leap ;; TYPE old-hindu-lunar-day old-hindu-lunar-year) ;; TYPE -> old-hindu-lunar-date (list month leap day year)) (defun hindu-lunar-date (month leap-month day leap-day year) ;; TYPE (hindu-lunar-month hindu-lunar-leap-month ;; TYPE hindu-lunar-day hindu-lunar-leap-day ;; TYPE hindu-lunar-year) -> hindu-lunar-date (list month leap-month day leap-day year)) **) (* Not needed in Mathematica; simply write [, , ... ] or TimeOfDay[, , ] *) (** (defconstant jd-start ;; TYPE moment ;; Fixed time of start of the julian day number. -1721424.5) **) JDStart[] = -1721424.5 (** (defun moment-from-jd (jd) ;; TYPE julian-day-number -> moment ;; Fixed time of astronomical (julian) day number jd. (+ jd jd-start)) **) MomentFromJD[jd_] := jd + JDStart[] (** (defun fixed-from-jd (jd) ;; TYPE julian-day-number -> fixed-date ;; Fixed date of astronomical (julian) day number jd. (floor (moment-from-jd jd))) **) FixedFromJD[jd_] := Floor[MomentFromJD[jd]] (** (defun jd-from-moment (moment) ;; TYPE moment -> julian-day-number ;; Astronomical (julian) day number of fixed moment moment. (- moment jd-start)) **) JDFromMoment[moment_] := moment - JDStart[] (* ADDED *) JDFromFixed[date_] := JDFromMoment[date] (** ;;;; Section: Gregorian Calendar **) (** (defconstant gregorian-epoch ;; TYPE fixed-date ;; Fixed date of start of the (proleptic) Gregorian ;; calendar. 1) **) GregorianEpoch[] = 1 (* (defconstant january ;; TYPE standard-month ;; January on Julian/Gregorian calendar. 1) (defconstant february ;; TYPE standard-month ;; February on Julian/Gregorian calendar. (1+ january)) (defconstant march ;; TYPE standard-month ;; March on Julian/Gregorian calendar. (+ january 2)) (defconstant april ;; TYPE standard-month ;; April on Julian/Gregorian calendar. (+ january 3)) (defconstant may ;; TYPE standard-month ;; May on Julian/Gregorian calendar. (+ january 4)) (defconstant june ;; TYPE standard-month ;; June on Julian/Gregorian calendar. (+ january 5)) (defconstant july ;; TYPE standard-month ;; July on Julian/Gregorian calendar. (+ january 6)) (defconstant august ;; TYPE standard-month ;; August on Julian/Gregorian calendar. (+ january 7)) (defconstant september ;; TYPE standard-month ;; September on Julian/Gregorian calendar. (+ january 8)) (defconstant october ;; TYPE standard-month ;; October on Julian/Gregorian calendar. (+ january 9)) (defconstant november ;; TYPE standard-month ;; November on Julian/Gregorian calendar. (+ january 10)) (defconstant december ;; TYPE standard-month ;; December on Julian/Gregorian calendar. (+ january 11)) **) January[] = 1 February[] = January[] + 1 March[] = January[] + 2 April[] = January[] + 3 May[] = January[] + 4 June[] = January[] + 5 July[] = January[] + 6 August[] = January[] + 7 September[] = January[] + 8 October[] = January[] + 9 November[] = January[] + 10 December[] = January[] + 11 (* ADDED *) NameFromMonth[month_] := NameFromNumber[month, MonthNames[]] (** (defun gregorian-leap-year? (g-year) ;; TYPE gregorian-year -> boolean ;; True if year is a leap year on the Gregorian calendar. (and (= (mod g-year 4) 0) (not (member (mod g-year 400) (list 100 200 300))))) **) LeapYearQ[Gregorian, year_] := Mod[year, 4] == 0 && !MemberQ[{100, 200, 300}, Mod[year, 400]] (** (defun fixed-from-gregorian (g-date) ;; TYPE gregorian-date -> fixed-date ;; Fixed date equivalent to the Gregorian date. (let* ((month (standard-month g-date)) (day (standard-day g-date)) (year (standard-year g-date))) (+ (1- gregorian-epoch); Days before start of calendar ( * 365 (1- year)); Ordinary days since epoch (quotient (1- year) 4); Julian leap days since epoch... (- ; ...minus century years since epoch... (quotient (1- year) 100)) (quotient ; ...plus years since epoch divisible... (1- year) 400) ; ...by 400. (quotient ; Days in prior months this year... (- ( * 367 month) 362); ...assuming 30-day Feb 12) (if (<= month 2) ; Correct for 28- or 29-day Feb 0 (if (gregorian-leap-year? year) -1 -2)) day))) ; Days so far this month. **) ToFixed[date_Gregorian] := Module[{month = CMonth[date], year = CYear[date]}, (GregorianEpoch[] - 1) + 365 (year - 1) + Quotient[year - 1, 4] - Quotient[year - 1, 100] + Quotient[year - 1, 400] + Quotient[367 month - 362, 12] + If[month <= 2, 0, If[LeapYearQ[Gregorian, year], -1, -2]] + CDay[date] ] (** (defun gregorian-year-from-fixed (date) ;; TYPE fixed-date -> gregorian-year ;; Gregorian year corresponding to the fixed date. (let* ((d0 ; Prior days. (- date gregorian-epoch)) (n400 ; Completed 400-year cycles. (quotient d0 146097)) (d1 ; Prior days not in n400. (mod d0 146097)) (n100 ; 100-year cycles not in n400. (quotient d1 36524)) (d2 ; Prior days not in n400 or n100. (mod d1 36524)) (n4 ; 4-year cycles not in n400 or n100. (quotient d2 1461)) (d3 ; Prior days not in n400, n100, or n4. (mod d2 1461)) (n1 ; Years not in n400, n100, or n4. (quotient d3 365)) (d4 ; Prior days not in n400, n100, n4, or n1. (1+ (mod d3 365))) (year (+ ( * 400 n400) ( * 100 n100) ( * 4 n4) n1))) (if (or (= n100 4) (= n1 4)) year ; Date is December 31 in year. (1+ year)))); Date is ordinal day (1+ d4) in (1+ year). **) GregorianYearFromFixed[date_] := Module[{d0, n400, d1, n100, d2, n4, d3, n1, d4, year}, d0 = date - GregorianEpoch[]; n400 = Quotient[d0, 146097]; d1 = Mod[d0, 146097]; n100 = Quotient[d1, 36524]; d2 = Mod[d1, 36524]; n4 = Quotient[d2, 1461]; d3 = Mod[d2, 1461]; n1 = Quotient[d3, 365]; d4 = Mod[d3, 365] + 1; year = 400 n400 + 100 n100 + 4 n4 + n1; If[n100 == 4 || n1 == 4, year, year + 1] ] (** (defun gregorian-from-fixed (date) ;; TYPE fixed-date -> gregorian-date ;; Gregorian (month day year) corresponding to fixed date. (let* ((year (gregorian-year-from-fixed date)) (prior-days; This year (- date (fixed-from-gregorian (gregorian-date january 1 year)))) (correction; To simulate a 30-day Feb (if (< date (fixed-from-gregorian (gregorian-date march 1 year))) 0 (if (gregorian-leap-year? year) 1 2))) (month ; Assuming a 30-day Feb (quotient (+ ( * 12 (+ prior-days correction)) 373) 367)) (day ; Calculate the day by subtraction. (1+ (- date (fixed-from-gregorian (gregorian-date month 1 year)))))) (gregorian-date month day year))) **) Gregorian[date_Integer] := Module[{year, priorDays, correction, month, day}, year = GregorianYearFromFixed[date]; priorDays = date - ToFixed[ Gregorian[January[], 1, year] ]; correction = If[ date < ToFixed[ Gregorian[March[], 1, year] ], 0, If[ LeapYearQ[Gregorian, year], 1, 2] ]; month = Quotient[12 (priorDays + correction) + 373, 367]; day = date - ToFixed[ Gregorian[month, 1, year] ] + 1; Gregorian[month, day, year] ] (** (defun gregorian-date-difference (g-date1 g-date2) ;; TYPE (gregorian-date gregorian-date) -> integer ;; Number of days from Gregorian date g-date1 until g-date2. (- (fixed-from-gregorian g-date2) (fixed-from-gregorian g-date1))) **) GregorianDateDifference[date1_, date2_] := ToFixed[date2] - ToFixed[date1] (** (defun day-number (g-date) ;; TYPE gregorian-date -> non-negative-integer ;; Day number in year of Gregorian date g-date. (gregorian-date-difference (gregorian-date december 31 (1- (standard-year g-date))) g-date)) **) DayNumber[] = {Gregorian} DayNumber[date_Gregorian] := GregorianDateDifference[Gregorian[December[], 31, CYear[date] - 1], date] (** (defun days-remaining (g-date) ;; TYPE gregorian-date -> non-negative-integer ;; Days remaining in year after Gregorian date g-date. (gregorian-date-difference g-date (gregorian-date december 31 (standard-year g-date)))) **) DaysRemaining[] = {Gregorian} DaysRemaining[date_Gregorian] := GregorianDateDifference[date, Gregorian[December[], 31, CYear[date]]] (** (defun independence-day (year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of American Independence Day in ;; Gregorian year. (fixed-from-gregorian (gregorian-date july 4 year))) **) IndependenceDay[year_] := ToFixed[Gregorian[July[], 4, year]]; (** (defun kday-on-or-before (date k) ;; TYPE (fixed-date weekday) -> fixed-date ;; Fixed date of the k-day on or before fixed date. ;; k=0 means Sunday, k=1 means Monday, and so on. (- date (day-of-week-from-fixed (- date k)))) **) KDayOnOrBefore[date_, k_] := date - DayOfWeekFromFixed[date - k] (** (defun kday-on-or-after (date k) ;; TYPE (fixed-date weekday) -> fixed-date ;; Fixed date of the k-day on or after fixed date. ;; k=0 means Sunday, k=1 means Monday, and so on. (kday-on-or-before (+ date 6) k)) **) KDayOnOrAfter[date_, k_] := KDayOnOrBefore[date + 6, k] (** (defun kday-nearest (date k) ;; TYPE (fixed-date weekday) -> fixed-date ;; Fixed date of the k-day nearest fixed date. k=0 ;; means Sunday, k=1 means Monday, and so on. (kday-on-or-before (+ date 3) k)) **) KDayNearest[date_, k_] := KDayOnOrBefore[date + 3, k] (** (defun kday-after (date k) ;; TYPE (fixed-date weekday) -> fixed-date ;; Fixed date of the k-day after fixed date. k=0 ;; means Sunday, k=1 means Monday, and so on. (kday-on-or-before (+ date 7) k)) **) KDayAfter[date_, k_] := KDayOnOrBefore[date + 7, k] (** (defun kday-before (date k) ;; TYPE (fixed-date weekday) -> fixed-date ;; Fixed date of the k-day before fixed date. k=0 ;; means Sunday, k=1 means Monday, and so on. (kday-on-or-before (1- date) k)) **) KDayBefore[date_, k_] := KDayOnOrBefore[date - 1, k] (** (defun nth-kday (n k date) ;; TYPE (integer weekday gregorian-date) -> fixed-date ;; Fixed date of n-th k-day after Gregorian date. If ;; n>0, return the n-th k-day on or after date. ;; If n<0, return the n-th k-day on or before date. ;; A k-day of 0 means Sunday, 1 means Monday, and so on. (if (> n 0) (+ ( * 7 n) (kday-before (fixed-from-gregorian date) k)) (+ ( * 7 n) (kday-after (fixed-from-gregorian date) k)))) **) NthKDay[n_, k_, date_] := If[n > 0, KDayBefore[ToFixed[date], k] + 7 n, KDayAfter[ToFixed[date], k] + 7 n] (** (defconstant first ;; TYPE integer ;; Index for selecting a k-day. 1) (defconstant last ;; TYPE integer ;; Index for selecting a k-day. -1) **) OFirst[] = 1 OLast[] = -1 (** (defun labor-day (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of American Labor Day in Gregorian ;; year--the first Monday in September. (nth-kday first monday (gregorian-date september 1 g-year))) **) LaborDay[year_] := NthKDay[OFirst[], Monday[], Gregorian[September[], 1, year]] (** (defun memorial-day (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of American Memorial Day in Gregorian ;; year--the last Monday in May. (nth-kday last monday (gregorian-date may 31 g-year))) **) MemorialDay[year_] := NthKDay[OLast[], Monday[], Gregorian[May[], 31, year]] (** (defun election-day (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of American Election Day in Gregorian ;; year--the Tuesday after the first Monday in November. (nth-kday first tuesday (gregorian-date november 2 g-year))) **) ElectionDay[year_] := NthKDay[OFirst[], Tuesday[], Gregorian[November[], 2, year]] (** (defun daylight-savings-start (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of the start of American daylight savings ;; time in Gregorian year--the first Sunday in April. (nth-kday first sunday (gregorian-date april 1 g-year))) **) DaylightSavingsStart[year_] := NthKDay[OFirst[], Sunday[], Gregorian[April[], 1, year]] (** (defun daylight-savings-end (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of the end of American daylight savings time ;; in Gregorian year--the last Sunday in October. (nth-kday last sunday (gregorian-date october 31 g-year))) **) DaylightSavingsEnd[year_] := NthKDay[OLast[], Sunday[], Gregorian[October[], 31, year]] (** (defun christmas (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Christmas in Gregorian year. (fixed-from-gregorian (gregorian-date december 25 g-year))) **) Christmas[year_] := ToFixed[Gregorian[December[], 25, year]] (** (defun advent (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Advent in Gregorian year. (kday-nearest (fixed-from-gregorian (gregorian-date november 30 g-year)) sunday)) **) Advent[year_] := KDayNearest[ToFixed[Gregorian[November[], 30, year]], Sunday[]] (** (defun epiphany (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Epiphany in Gregorian year. (+ 12 (christmas (1- g-year)))) **) Epiphany[year_] := Christmas[year - 1] + 12 (** ;;;; Section: ISO Calendar **) (** (defun fixed-from-iso (i-date) ;; TYPE iso-date -> fixed-date ;; Fixed date equivalent to ISO (week day year). (let* ((week (iso-week i-date)) (day (iso-day i-date)) (year (iso-year i-date))) ;; Add fixed date of Sunday preceding date plus day ;; in week. (+ (nth-kday week sunday (gregorian-date december 28 (1- year))) day))) **) ToFixed[date_ISO] := NthKDay[CWeek[date], Sunday[], Gregorian[December[], 28, CYear[date] - 1]] + CDay[date] (** (defun iso-from-fixed (date) ;; TYPE fixed-date -> iso-date ;; ISO (week day year) corresponding to the fixed date. (let* ((approx ; Year may be one too small. (gregorian-year-from-fixed (- date 3))) (year (if (>= date (fixed-from-iso (iso-date 1 1 (1+ approx)))) (1+ approx) approx)) (week (1+ (quotient (- date (fixed-from-iso (iso-date 1 1 year))) 7))) (day (adjusted-mod date 7))) (iso-date week day year))) **) ISO[date_Integer] := Module[{approx, year, week, day}, approx = GregorianYearFromFixed[date - 3]; year = If[ date >= ToFixed[ISO[1, 1, approx + 1]], approx + 1, approx]; week = Quotient[date - ToFixed[ ISO[1, 1, year] ], 7] + 1; day = AdjustedMod[date, 7]; ISO[week, day, year] ] (** ;;;; Section: Julian Calendar **) (** (defconstant julian-epoch ;; TYPE fixed-date ;; Fixed date of start of the Julian calendar. (fixed-from-gregorian (gregorian-date december 30 0))) **) JulianEpoch[] = ToFixed[Gregorian[December[], 30, 0]] (** (defun bce (n) ;; TYPE standard-year -> julian-year ;; Negative value to indicate a BCE Julian year. (- n)) **) BCE[n_] := - n (** (defun ce (n) ;; TYPE standard-year -> julian-year ;; Positive value to indicate a CE Julian year. n) **) CE[n_] := n (** (defun julian-leap-year? (j-year) ;; TYPE julian-year -> boolean ;; True if year is a leap year on the Julian calendar. (= (mod j-year 4) (if (> j-year 0) 0 3))) **) LeapYearQ[Julian, year_] := Mod[year, 4] == If[year > 0, 0, 3] (** (defun fixed-from-julian (j-date) ;; TYPE julian-date -> fixed-date ;; Fixed date equivalent to the Julian date. (let* ((month (standard-month j-date)) (day (standard-day j-date)) (year (standard-year j-date)) (y (if (< year 0) (1+ year) ; No year zero year))) (+ (1- julian-epoch) ; Days before start of calendar ( * 365 (1- y)) ; Ordinary days since epoch. (quotient (1- y) 4); Leap days since epoch... (quotient ; Days in prior months this year... (- ( * 367 month) 362); ...assuming 30-day Feb 12) (if (<= month 2) ; Correct for 28- or 29-day Feb 0 (if (julian-leap-year? year) -1 -2)) day))) ; Days so far this month. **) ToFixed[date_Julian] := Module[{month, day, year, y}, month = CMonth[date]; day = CDay[date]; year = CYear[date]; y = If[year < 0, year + 1, year]; JulianEpoch[] - 1 + 365 (y - 1) + Quotient[y - 1, 4] + Quotient[367 month - 362, 12] + If[month <= 2, 0, If[LeapYearQ[Julian, year], -1, -2]] + day ] (** (defun julian-from-fixed (date) ;; TYPE fixed-date -> julian-date ;; Julian (month day year) corresponding to fixed date. (let* ((approx ; Nominal year. (quotient (+ ( * 4 (- date julian-epoch)) 1464) 1461)) (year (if (<= approx 0) (1- approx) ; No year 0. approx)) (prior-days; This year (- date (fixed-from-julian (julian-date january 1 year)))) (correction; To simulate a 30-day Feb (if (< date (fixed-from-julian (julian-date march 1 year))) 0 (if (julian-leap-year? year) 1 2))) (month ; Assuming a 30-day Feb (quotient (+ ( * 12 (+ prior-days correction)) 373) 367)) (day ; Calculate the day by subtraction. (1+ (- date (fixed-from-julian (julian-date month 1 year)))))) (julian-date month day year))) **) Julian[date_Integer] := Module[{approx, year, priorDays, correction, month}, approx = Quotient[4 (date - JulianEpoch[]) + 1464, 1461]; year = If[approx <= 0, approx - 1, approx]; priorDays = date - ToFixed[Julian[January[], 1, year]]; correction = If[date < ToFixed[Julian[March[], 1, year]], 0, If[LeapYearQ[Julian, year], 1, 2]]; month = Quotient[12 (priorDays + correction) + 373, 367]; day = date - ToFixed[Julian[month, 1, year]] + 1; Julian[month, day, year] ] (** (defun nicaean-rule-easter (j-year) ;; TYPE julian-year -> fixed-date ;; Fixed date of Easter in positive Julian year, according ;; to the rule of the Council of Nicaea. (let* ((shifted-epact ; Age of moon for April 5. (mod (+ 14 ( * 11 (mod j-year 19))) 30)) (paschal-moon ; Day after full moon on ; or after March 21. (- (fixed-from-julian (julian-date april 19 j-year)) shifted-epact))) ;; Return the Sunday following the Paschal moon (kday-after paschal-moon sunday))) **) NicaeanRuleEaster[jYear_] := Module[{shiftedEpact, paschalMoon}, shiftedEpact = Mod[14 + 11 Mod[jYear, 19], 30]; paschalMoon = ToFixed[Julian[April[], 19, jYear]] - shiftedEpact; KDayAfter[paschalMoon, Sunday[]] ] (** (defun easter (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Easter in Gregorian year. (let* ((century (1+ (quotient g-year 100))) (shifted-epact ; Age of moon for April 5... (mod (+ 14 ( * 11 (mod g-year 19)); ...by Nicaean rule (- ;...corrected for the Gregorian century rule (quotient ( * 3 century) 4)) (quotient; ...corrected for Metonic ; cycle inaccuracy. (+ 5 ( * 8 century)) 25)) 30)) (adjusted-epact ; Adjust for 29.5 day month. (if (or (= shifted-epact 0) (and (= shifted-epact 1) (< 10 (mod g-year 19)))) (1+ shifted-epact) shifted-epact)) (paschal-moon; Day after full moon on ; or after March 21. (- (fixed-from-gregorian (gregorian-date april 19 g-year)) adjusted-epact))) ;; Return the Sunday following the Paschal moon. (kday-after paschal-moon sunday))) **) Easter[gYear_] := Module[{century, shiftedEpact, adjustedEpact, paschalMoon}, century = 1 + Quotient[gYear, 100]; shiftedEpact = Mod[14 + 11 Mod[gYear, 19] - Quotient[3 century, 4] + Quotient[5 + 8 century, 25], 30]; adjustedEpact = If[shiftedEpact == 0 || (shiftedEpact == 1 && 10 < Mod[gYear, 19]), shiftedEpact + 1, shiftedEpact]; paschalMoon = ToFixed[Gregorian[April[], 19, gYear]] - adjustedEpact; KDayAfter[paschalMoon, Sunday[]] ] (** (defun pentecost (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Pentecost in Gregorian year. (+ (easter g-year) 49)) **) Pentecost[gYear_] := Easter[gYear] + 49 (** (defun julian-in-gregorian (j-month j-day g-year) ;; TYPE (julian-month julian-day gregorian-year) ;; TYPE -> list-of-fixed-dates ;; List of the fixed dates of Julian month, day ;; that occur in Gregorian year. (let* ((jan1 (fixed-from-gregorian (gregorian-date january 1 g-year))) (dec31 (fixed-from-gregorian (gregorian-date december 31 g-year))) (y (standard-year (julian-from-fixed jan1))) ;; The possible occurrences in one year are (date1 (fixed-from-julian (julian-date j-month j-day y))) (date2 (fixed-from-julian (julian-date j-month j-day (1+ y))))) (append (if ; date1 occurs in current year (<= jan1 date1 dec31) ;; Then that date; otherwise, none (list date1) nil) (if ; date2 occurs in current year (<= jan1 date2 dec31) ;; Then that date; otherwise, none (list date2) nil)))) **) JulianInGregorian[jMonth_, jDay_, gYear_] := Module[{jan1, dec31, y, date1, date2, result}, jan1 = ToFixed[Gregorian[January[], 1, gYear]]; dec31 = ToFixed[Gregorian[December[], 31, gYear]]; y = CYear[Julian[jan1]]; date1 = ToFixed[Julian[jMonth, jDay, y]]; date2 = ToFixed[Julian[jMonth, jDay, y + 1]]; result = {}; If[LessEqual[jan1, date1, dec31], result = Append[result, date1]]; If[LessEqual[jan1, date2, dec31], result = Append[result, date2]]; result ] (** (defun eastern-orthodox-christmas (g-year) ;; TYPE gregorian-year -> list-of-fixed-dates ;; List of zero or one fixed dates of Eastern Orthodox ;; Christmas in Gregorian year. (julian-in-gregorian december 25 g-year)) **) EasternOrthodoxChristmas[gYear_] := JulianInGregorian[December[], 25, gYear] (** ;;;; Section: Coptic and Ethiopic Calendars **) (** (defconstant coptic-epoch ;; TYPE fixed-date ;; Fixed date of start of the Coptic calendar. (fixed-from-julian (julian-date august 29 (ce 284)))) **) CopticEpoch[] = ToFixed[Julian[August[], 29, CE[284]]] (** (defun coptic-leap-year? (c-year) ;; TYPE coptic-year -> boolean ;; True if year is a leap year on the Coptic calendar. (= (mod c-year 4) 3)) **) LeapYearQ[Coptic, year_] := Mod[year, 4] == 3 (** (defun fixed-from-coptic (c-date) ;; TYPE coptic-date -> fixed-date ;; Fixed date of Coptic date. (let* ((month (standard-month c-date)) (day (standard-day c-date)) (year (standard-year c-date))) (+ coptic-epoch -1 ; Days before start of calendar ( * 365 (1- year)); Ordinary days in prior years (quotient year 4); Leap days in prior years ( * 30 (1- month)); Days in prior months this year day))) ; Days so far this month **) ToFixed[date_Coptic] := Module[{year}, year = CYear[date]; CopticEpoch[] - 1 + 365 (year - 1) + Quotient[year, 4] + 30 (CMonth[date] - 1) + CDay[date] ] (** (defun coptic-from-fixed (date) ;; TYPE fixed-date -> coptic-date ;; Coptic equivalent of fixed date. (let* ((year ; Calculate the year by cycle-of-years formula (quotient (+ ( * 4 (- date coptic-epoch)) 1463) 1461)) (month; Calculate the month by division. (1+ (quotient (- date (fixed-from-coptic (coptic-date 1 1 year))) 30))) (day ; Calculate the day by subtraction. (- date -1 (fixed-from-coptic (coptic-date month 1 year))))) (coptic-date month day year))) **) Coptic[date_Integer] := Module[{year, month, day}, year = Quotient[4 (date - CopticEpoch[]) + 1463, 1461]; month = Quotient[date - ToFixed[Coptic[1, 1, year]], 30] + 1; day = date + 1 - ToFixed[Coptic[month, 1, year]]; Coptic[month, day, year] ] (** (defconstant ethiopic-epoch ;; TYPE fixed-date ;; Fixed date of start of the Ethiopic calendar. (fixed-from-julian (julian-date august 29 (ce 7)))) **) EthiopicEpoch[] = ToFixed[Julian[August[], 29, CE[7]]] (** (defun fixed-from-ethiopic (e-date) ;; TYPE ethiopic-date -> fixed-date ;; Fixed date of Ethiopic date. (let* ((month (standard-month e-date)) (day (standard-day e-date)) (year (standard-year e-date))) (+ ethiopic-epoch (- (fixed-from-coptic (coptic-date month day year)) coptic-epoch)))) **) ToFixed[date_Ethiopic] := ToFixed[Coptic[CMonth[date], CDay[date], CYear[date]]] - CopticEpoch[] + EthiopicEpoch[] (** (defun ethiopic-from-fixed (date) ;; TYPE fixed-date -> ethiopic-date ;; Ethiopic equivalent of fixed date. (coptic-from-fixed (+ date (- coptic-epoch ethiopic-epoch)))) **) Ethiopic[date_Integer] := Apply[Ethiopic, Coptic[date + CopticEpoch[] - EthiopicEpoch[]]] (** (defun coptic-in-gregorian (c-month c-day g-year) ;; TYPE (coptic-month coptic-day gregorian-year) ;; TYPE -> list-of-fixed-dates ;; List of the fixed dates of Coptic month, day ;; that occur in Gregorian year. (let* ((jan1 (fixed-from-gregorian (gregorian-date january 1 g-year))) (dec31 (fixed-from-gregorian (gregorian-date december 31 g-year))) (y (standard-year (coptic-from-fixed jan1))) ;; The possible occurrences in one year are (date1 (fixed-from-coptic (coptic-date c-month c-day y))) (date2 (fixed-from-coptic (coptic-date c-month c-day (1+ y))))) (append (if ; date1 occurs in current year (<= jan1 date1 dec31) ;; Then that date; otherwise, none (list date1) nil) (if ; date2 occurs in current year (<= jan1 date2 dec31) ;; Then that date; otherwise, none (list date2) nil)))) **) CopticInGregorian[cMonth_, cDay_, gYear_] := Module[{jan1, dec31, y, date1, date2, result}, jan1 = ToFixed[Gregorian[January[], 1, gYear]]; dec31 = ToFixed[Gregorian[December[], 31, gYear]]; y = CYear[Coptic[jan1]]; date1 = ToFixed[Coptic[cMonth, cDay, y]]; date2 = ToFixed[Coptic[cMonth, cDay, y + 1]]; result = {}; If[LessEqual[jan1, date1, dec31], result = Append[result, date1]]; If[LessEqual[jan1, date2, dec31], result = Append[result, date2]]; result ] (** (defun coptic-christmas (g-year) ;; TYPE gregorian-year -> list-of-fixed-dates ;; List of zero or one fixed dates of Coptic Christmas ;; in Gregorian year. (coptic-in-gregorian 4 29 g-year)) **) CopticChristmas[gYear_] := CopticInGregorian[4, 29, gYear] (** ;;;; Section: Islamic Calendar **) (** (defconstant islamic-epoch ;; TYPE fixed-date ;; Fixed date of start of the Islamic calendar. (fixed-from-julian (julian-date july 16 (ce 622)))) **) IslamicEpoch[] = ToFixed[Julian[July[], 16, CE[622]]] (** (defun islamic-leap-year? (i-year) ;; TYPE islamic-year -> boolean ;; True if year is an Islamic leap year. (< (mod (+ 14 ( * 11 i-year)) 30) 11)) **) LeapYearQ[Islamic, year_] := Mod[11 year + 14, 30] < 11 (** (defun fixed-from-islamic (i-date) ;; TYPE islamic-date -> fixed-date ;; Fixed date equivalent to Islamic date. (let* ((month (standard-month i-date)) (day (standard-day i-date)) (year (standard-year i-date))) (+ day ; Days so far this month. (ceiling ; Days in prior months. ( * 29.5 (1- month))) ( * (1- year) 354) ; Nonleap days in prior years. (quotient ; Leap days in prior years. (+ 3 ( * 11 year)) 30) islamic-epoch -1))) ; Days before start of calendar. **) ToFixed[date_Islamic] := Module[{year}, year = CYear[date]; CDay[date] + Ceiling[29.5 (CMonth[date] - 1)] + (year - 1) 354 + Quotient[3 + 11 year, 30] + IslamicEpoch[] - 1 ] (** (defun islamic-from-fixed (date) ;; TYPE fixed-date -> islamic-date ;; Islamic date (month day year) corresponding to fixed ;; date. (let* ((year ; Divide elapsed days by average year length. (quotient (+ ( * 30 (- date islamic-epoch)) 10646) 10631)) (month ; Months alternate between 29 and 30 days (min 12 ; Last month can be longer (1+ (ceiling (/ (- date 29 (fixed-from-islamic (islamic-date 1 1 year))) 29.5))))) (day ; Calculate the day by subtraction. (1+ (- date (fixed-from-islamic (islamic-date month 1 year)))))) (islamic-date month day year))) **) Islamic[date_Integer] := Module[{year, month, day}, year = Quotient[30 (date - IslamicEpoch[]) + 10646, 10631]; month = Min[12, 1 + Ceiling[(date - 29 - ToFixed[Islamic[1, 1, year]]) / 29.5]]; day = 1 + date - ToFixed[Islamic[month, 1, year]]; Islamic[month, day, year] ] (** (defun islamic-in-gregorian (i-month i-day g-year) ;; TYPE (islamic-month islamic-day gregorian-year) ;; TYPE -> list-of-fixed-dates ;; List of the fixed dates of Islamic month, day ;; that occur in Gregorian year. (let* ((jan1 (fixed-from-gregorian (gregorian-date january 1 g-year))) (dec31 (fixed-from-gregorian (gregorian-date december 31 g-year))) (y (standard-year (islamic-from-fixed jan1))) ;; The possible occurrences in one year are (date1 (fixed-from-islamic (islamic-date i-month i-day y))) (date2 (fixed-from-islamic (islamic-date i-month i-day (1+ y)))) (date3 (fixed-from-islamic (islamic-date i-month i-day (+ 2 y))))) ;; Combine in one list those that occur in current year (append (if (<= jan1 date1 dec31) (list date1) nil) (if (<= jan1 date2 dec31) (list date2) nil) (if (<= jan1 date3 dec31) (list date3) nil)))) **) IslamicInGregorian[iMonth_, iDay_, gYear_] := Module[{jan1, dec31, y, date1, date2, date3, result}, jan1 = ToFixed[Gregorian[January[], 1, gYear]]; dec31 = ToFixed[Gregorian[December[], 31, gYear]]; y = CYear[Islamic[jan1]]; date1 = ToFixed[Islamic[iMonth, iDay, y]]; date2 = ToFixed[Islamic[iMonth, iDay, y + 1]]; date3 = ToFixed[Islamic[iMonth, iDay, y + 2]]; result = {}; If[LessEqual[jan1, date1, dec31], result = Append[result, date1]]; If[LessEqual[jan1, date2, dec31], result = Append[result, date2]]; If[LessEqual[jan1, date3, dec31], result = Append[result, date3]]; result ] (** (defun mawlid-an-nabi (g-year) ;; TYPE gregorian-year -> list-of-fixed-dates ;; List of fixed dates of Mawlid-an-Nabi occurring in ;; Gregorian year. (islamic-in-gregorian 3 12 g-year)) **) MawlidAnNabi[gYear_] := IslamicInGregorian[3, 12, gYear] (** ;;;; Section: Bahai Calendar **) (** (defconstant bahai-epoch ;; TYPE fixed-date ;; Fixed date of start of Bahai calendar. (fixed-from-gregorian (gregorian-date march 21 1844))) **) BahaiEpoch[] = ToFixed[Gregorian[March[], 21, 1844]] (** (defun fixed-from-bahai (b-date) ;; TYPE bahai-date -> fixed-date ;; Fixed date equivalent to the Bahai date b-date. (let* ((major (bahai-major b-date)) (cycle (bahai-cycle b-date)) (year (bahai-year b-date)) (month (bahai-month b-date)) (day (bahai-day b-date)) (g-year; Corresponding Gregorian year. (+ ( * 361 (1- major)) ( * 19 (1- cycle)) year -1 (gregorian-year-from-fixed bahai-epoch)))) (+ (fixed-from-gregorian ; Prior years. (gregorian-date march 20 g-year)) ( * 19 (1- month)) ; Elapsed months. ;; Subtract 14 or 15 if counted ayyam-i-ha. (if (/= month 20) 0 (if (gregorian-leap-year? (1+ g-year)) -14 -15)) day))) ; Days in current month. **) ToFixed[date_Bahai] := Module[{month, gYear}, month = CMonth[date]; gYear = 361 (CMajor[date] - 1) + 19 (CCycle[date] - 1) + CYear[date] - 1 + GregorianYearFromFixed[BahaiEpoch[]]; ToFixed[Gregorian[March[], 20, gYear]] + 19 (month - 1) + If[month =!= 20, 0, If[LeapYearQ[Gregorian, gYear + 1], -14, -15]] + CDay[date] ] (** (defun bahai-from-fixed (date) ;; TYPE fixed-date -> bahai-date ;; Bahai (month day cycle year) corresponding to fixed ;; date. (let* ((g-year (gregorian-year-from-fixed date)) (start ; 1844 (gregorian-year-from-fixed bahai-epoch)) (years ; Since start of Bahai calendar. (- g-year start (if (<= (fixed-from-gregorian (gregorian-date january 1 g-year)) date (fixed-from-gregorian (gregorian-date march 20 g-year))) 1 0))) (major (1+ (quotient years 361))) (cycle (1+ (quotient (mod years 361) 19))) (year (1+ (mod years 19))) (days; Since start of year (- date (fixed-from-bahai (bahai-date major cycle year 1 1)))) (month (if (>= date (fixed-from-bahai (bahai-date major cycle year 20 1))) 20 (1+ (quotient days 19)))) (day (- date -1 (fixed-from-bahai (bahai-date major cycle year month 1))))) (bahai-date major cycle year month day))) **) Bahai[date_Integer] := Module[{gYear, start, years, major, cycle, year, days, month, day}, gYear = GregorianYearFromFixed[date]; start = GregorianYearFromFixed[BahaiEpoch[]]; years = gYear - start - If[ToFixed[Gregorian[January[], 1, gYear]] <= date <= ToFixed[Gregorian[March[], 20, gYear]], 1, 0]; major = 1 + Quotient[years, 361]; cycle = 1 + Quotient[Mod[years, 361], 19]; year = 1 + Mod[years, 19]; days = date - ToFixed[Bahai[major, cycle, year, 1, 1]]; month = If[date >= ToFixed[Bahai[major, cycle, year, 20, 1]], 20, 1 + Quotient[days, 19]]; day = date + 1 - ToFixed[Bahai[major, cycle, year, month, 1]]; Bahai[major, cycle, year, month, day] ] (** (defun bahai-new-year (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Bahai New Year in Gregorian year. (fixed-from-gregorian (gregorian-date march 21 g-year))) **) BahaiNewYear[gYear_] := ToFixed[Gregorian[March[], 21, gYear]] (** (defun feast-of-ridvan (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Bahai New Year in Gregorian year. (let* ((years (- g-year (gregorian-year-from-fixed bahai-epoch))) (major (1+ (quotient years 361))) (cycle (1+ (quotient (mod years 361) 19))) (year (1+ (mod years 19)))) (fixed-from-bahai (bahai-date major cycle year 2 13)))) **) FeastOfRidvan[gYear_] := Module[{years, major, cycle, year}, years = gYear - GregorianYearFromFixed[BahaiEpoch[]]; major = 1 + Quotient[years, 361]; cycle = 1 + Quotient[Mod[years, 361], 19]; year = 1 + Mod[years, 19]; ToFixed[Bahai[major, cycle, year, 2, 13]] ] (** ;;;; Section: Persian Calendar **) (** (defun persian-leap-year? (p-year) ;; TYPE persian-year -> boolean ;; True if year is a leap year on the Persian calendar. (let* ((y ; Years since start of 2820-year cycles (if (< 0 p-year) (- p-year 474) (- p-year 473))); No year zero (year ; Equivalent year in the range 474...3263 (+ (mod y 2820) 474))) (< (mod ( * (+ year 38) 682) 2816) 682))) **) LeapYearQ[Persian, pYear_] := Module[{y, year}, y = If[0 < pYear, pYear - 474, pYear - 473]; year = Mod[y, 2820] + 474; Mod[(year + 38) 682, 2816] < 682 ] (** (defconstant persian-epoch ;; TYPE fixed-date ;; Fixed date of start of the Persian calendar. (fixed-from-julian (julian-date march 19 (ce 622)))) **) PersianEpoch[] = ToFixed[Julian[March[], 19, CE[622]]] (** (defun fixed-from-persian (p-date) ;; TYPE persian-date -> fixed-date ;; Fixed date equivalent to Persian date. (let* ((day (standard-day p-date)) (month (standard-month p-date)) (p-year (standard-year p-date)) (y ; Years since start of 2820-year cycle (if (< 0 p-year) (- p-year 474) (- p-year 473))); No year zero (year ; Equivalent year in the range 474...3263 (+ (mod y 2820) 474))) (+ (1- persian-epoch); Days before epoch ( * 1029983 ; Days in 2820-year cycles ; before Persian year 474 (quotient y 2820)) ( * 365 (1- year)) ; Nonleap days in prior years this ; 2820-year cycle (quotient ; Leap days in prior years this ; 2820-year cycle (- ( * 682 year) 110) 2816) (if (<= month 7) ; Days in prior months this year ( * 31 (1- month)) (+ ( * 30 (1- month)) 6)) day))) ; Days so far this month **) ToFixed[date_Persian] := Module[{day, month, pYear, y, year}, day = CDay[date]; month = CMonth[date]; pYear = CYear[date]; y = If[0 < pYear, pYear - 474, pYear - 473]; year = Mod[y, 2820] + 474; PersianEpoch[] - 1 + 1029983 Quotient[y, 2820] + 365 (year - 1) + Quotient[682 year - 110, 2816] + If[month <= 7, 31 (month - 1), 30 (month - 1) + 6] + day ] (** (defun persian-year-from-fixed (date) ;; TYPE fixed-date -> persian-year ;; Persian year corresponding to the fixed date. (let* ((d0 ; Prior days since start of 2820-year cycle ; beginning in Persian year 474 (- date (fixed-from-persian (persian-date 1 1 475)))) (n2820 ; Completed prior 2820-year cycles (quotient d0 1029983)) (d1 ; Prior days not in n2820--that is, days ; since start of last 2820-year cycle (mod d0 1029983)) (y2820 ; Years since start of last 2820-year cycle (if (= d1 1029982) ;; Last day of 2820-year cycle 2820 ;; Otherwise use cycle of years formula (quotient (+ ( * 2816 d1) 1031337) 1028522) ;; If ( * 2816 d1) causes integers that are ;; too-large, use instead: ;; (let ((a (floor d1 366)) ;; (b (mod d1 366))) ;; (+ 1 a (quotient ;; (+ ( * 2134 a) ( * 2816 b) 2815) ;; 1028522))) )) (year ; Years since Persian epoch (+ 474 ; Years before start of 2820-year cycles ( * 2820 n2820) ; Years in prior 2820-year cycles y2820))); Years since start of last 2820-year ; cycle (if (< 0 year) year (1- year)))); No year zero **) PersianYearFromFixed[date_] := Module[{d0, n2820, d1, y2820, year}, d0 = date - ToFixed[Persian[1, 1, 475]]; n2820 = Quotient[d0, 1029983]; d1 = Mod[d0, 1029983]; y2820 = If[d1 == 1029982, 2820, Quotient[2816 d1 + 1031337, 1028522]]; year = 474 + 2820 n2820 + y2820; If[0 < year, year, year - 1] ] (** (defun persian-from-fixed (date) ;; TYPE fixed-date -> persian-date ;; Persian (month day year) corresponding to fixed date. (let* ((year (persian-year-from-fixed date)) (day-of-year (1+ (- date (fixed-from-persian (persian-date 1 1 year))))) (month (if (<= day-of-year 186) (ceiling (/ day-of-year 31)) (ceiling (/ (- day-of-year 6) 30)))) (day ; Calculate the day by subtraction (- date (1- (fixed-from-persian (persian-date month 1 year)))))) (persian-date month day year))) **) Persian[date_Integer] := Module[{year, dayOfYear, month, day}, year = PersianYearFromFixed[date]; dayOfYear = 1 + date - ToFixed[Persian[1, 1, year]]; month = If[dayOfYear < 186, Ceiling[dayOfYear / 31], Ceiling[(dayOfYear - 6) / 30]]; day = date - (ToFixed[Persian[month, 1, year]] - 1); Persian[month, day, year] ] (** (defun naw-ruz (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Persian New Year (Naw-Ruz) in Gregorian ;; year. (let* ((persian-year (1+ (- g-year (gregorian-year-from-fixed persian-epoch))))) (fixed-from-persian (persian-date 1 1 (if (<= persian-year 0) ;; No Persian year 0 (1- persian-year) persian-year))))) **) NawRuz[gYear_] := Module[{persianYear}, persianYear = 1 + gYear - GregorianYearFromFixed[PersianEpoch[]]; ToFixed[Persian[1, 1, If[persianYear <= 0, persianYear - 1, persianYear]]] ] (** ;;;; Section: Hebrew Calendar **) (** (defconstant hebrew-epoch ;; TYPE fixed-date ;; Fixed date of start of the Hebrew calendar, that is, ;; Tishri 1, 1 AM. (fixed-from-julian (julian-date october 7 (bce 3761)))) **) HebrewEpoch[] = ToFixed[Julian[October[], 7, BCE[3761]]] (** (defun hebrew-leap-year? (h-year) ;; TYPE hebrew-year -> boolean ;; True if year is a leap year on Hebrew calendar. (< (mod (1+ ( * 7 h-year)) 19) 7)) **) LeapYearQ[Hebrew, hYear_] := Mod[1 + 7 hYear, 19] < 7 (** (defun last-month-of-hebrew-year (h-year) ;; TYPE hebrew-year -> hebrew-month ;; Last month of Hebrew year. (if (hebrew-leap-year? h-year) 13 12)) **) LastMonthOfHebrewYear[hYear_Integer] := If[LeapYearQ[Hebrew, hYear], 13, 12] (** (defun last-day-of-hebrew-month (h-month h-year) ;; TYPE (hebrew-month hebrew-year) -> hebrew-day ;; Last day of month in Hebrew year. (if (or (member h-month (list 2 4 6 10 13)) (and (= h-month 12) (not (hebrew-leap-year? h-year))) (and (= h-month 8) (not (long-heshvan? h-year))) (and (= h-month 9) (short-kislev? h-year))) 29 30)) **) LastDayOfHebrewMonth[hMonth_Integer, hYear_Integer] := If[MemberQ[{2, 4, 6, 10, 13}, hMonth] || ((hMonth == 12) && !LeapYearQ[Hebrew, hYear]) || ((hMonth == 8) && !LongHeshvanQ[hYear]) || ((hMonth == 9) && ShortKislevQ[hYear]), 29, 30] (** (defun hebrew-calendar-elapsed-days (h-year) ;; TYPE hebrew-year -> integer ;; Number of days elapsed from the (Sunday) noon prior ;; to the epoch of the Hebrew calendar to the mean ;; conjunction (molad) of Tishri of Hebrew year h-year, ;; or one day later. (let* ((months-elapsed ; Since start of Hebrew calendar. (quotient (- ( * 235 h-year) 234) 19)) (parts-elapsed; Fractions of days since prior noon. (+ 12084 ( * 13753 months-elapsed))) (day ; Whole days since prior noon. (+ ( * 29 months-elapsed) (quotient parts-elapsed 25920))) ;; If ( * 13753 months-elapsed) causes integers that ;; are to large, use instead: ;; (parts-elapsed ;; (+ 204 ( * 793 (mod months-elapsed 1080)))) ;; (hours-elapsed ;; (+ 11 ( * 12 months-elapsed) ;; ( * 793 (quotient months-elapsed 1080)) ;; (quotient parts-elapsed 1080))) ;; (day ;; (+ ( * 29 months-elapsed) ;; (quotient hours-elapsed 24))) ) (if (< (mod ( * 3 (1+ day)) 7) 3); Sun, Wed, or Fri (1+ day) ; Delay one day. day))) **) HebrewCalendarElapsedDays[hYear_Integer] := Module[{monthsElapsed, partsElapsed, day}, monthsElapsed = Quotient[235 hYear - 234, 19]; partsElapsed = 12084 + 13753 monthsElapsed; day = 29 monthsElapsed + Quotient[partsElapsed, 25920]; If[Mod[3 (day + 1), 7] < 3, day + 1, day] ] (** (defun hebrew-new-year-delay (h-year) ;; TYPE hebrew-year -> {0,1,2} ;; Delays to start of Hebrew year to keep ordinary year in ;; range 353-356 and leap year in range 383-386. (let* ((ny0 (hebrew-calendar-elapsed-days (1- h-year))) (ny1 (hebrew-calendar-elapsed-days h-year)) (ny2 (hebrew-calendar-elapsed-days (1+ h-year)))) (cond ((= (- ny2 ny1) 356) 2); Next year would be too long. ((= (- ny1 ny0) 382) 1); Previous year too short. (t 0)))) **) HebrewNewYearDelay[hYear_Integer] := Module[{ny0, ny1, ny2}, ny0 = HebrewCalendarElapsedDays[hYear - 1]; ny1 = HebrewCalendarElapsedDays[hYear]; ny2 = HebrewCalendarElapsedDays[hYear + 1]; Which[ny2 - ny1 == 356, 2, ny1 - ny0 == 382, 1, True, 0] ] (** (defun days-in-hebrew-year (h-year) ;; TYPE hebrew-year -> {353,354,355,383,384,385} ;; Number of days in Hebrew year. Calls fixed-from-hebrew ;; for value that does not in turn require ;; days-in-hebrew-year. (- (fixed-from-hebrew (hebrew-date 7 1 (1+ h-year))) (fixed-from-hebrew (hebrew-date 7 1 h-year)))) **) DaysInHebrewYear[hYear_Integer] := ToFixed[Hebrew[7, 1, hYear + 1]] - ToFixed[Hebrew[7, 1, hYear]]; (** (defun long-heshvan? (h-year) ;; TYPE hebrew-year -> boolean ;; True if Heshvan is long in Hebrew year. (= (mod (days-in-hebrew-year h-year) 10) 5)) **) LongHeshvanQ[hYear_Integer] := Mod[DaysInHebrewYear[hYear], 10] == 5 (** (defun short-kislev? (h-year) ;; TYPE hebrew-year -> boolean ;; True if Kislev is short in Hebrew year. (= (mod (days-in-hebrew-year h-year) 10) 3)) **) ShortKislevQ[hYear_Integer] := Mod[DaysInHebrewYear[hYear], 10] == 3 (** (defun fixed-from-hebrew (h-date) ;; TYPE hebrew-date -> fixed-date ;; Fixed date of Hebrew date. This function is designed ;; so that it works for Hebrew dates month, day, year even ;; if the month has fewer than day days--in that case the ;; function returns the (day-1)st day after month 1, year. ;; This property is required by the functions ;; hebrew-birthday and yahrzeit. (let* ((month (standard-month h-date)) (day (standard-day h-date)) (year (standard-year h-date))) (+ hebrew-epoch ; Days before fixed date 1. (hebrew-calendar-elapsed-days; Days in prior years. year) (hebrew-new-year-delay year) day -1 ; Days so far this month. (if ;; before Tishri (< month 7) ;; Then add days in prior months this year before ;; and after Nisan. (+ (sum (last-day-of-hebrew-month m year) m 7 (<= m (last-month-of-hebrew-year year))) (sum (last-day-of-hebrew-month m year) m 1 (< m month))) ;; Else add days in prior months this year (sum (last-day-of-hebrew-month m year) m 7 (< m month)))))) **) ToFixed[date_Hebrew] := Module[{month, day, year}, month = CMonth[date]; day = CDay[date]; year = CYear[date]; HebrewEpoch[] + HebrewCalendarElapsedDays[year] + HebrewNewYearDelay[year] + day - 1 + If[month < 7, MSum[(LastDayOfHebrewMonth[#, year])&, 7, (# <= LastMonthOfHebrewYear[year])&] + MSum[(LastDayOfHebrewMonth[#, year])&, 1, (# < month)&], MSum[(LastDayOfHebrewMonth[#, year])&, 7, (# < month)&]] ] (** (defun hebrew-from-fixed (date) ;; TYPE fixed-date -> hebrew-date ;; Hebrew (month day year) corresponding to fixed date. ;; The fraction can be approximated by 365.25. (let* ((approx ; Approximate year (may be off by 1) (quotient (- date hebrew-epoch) 35975351/98496)) ;; The value 35975351/98496, the average length of ;; a Hebrew year, can be approximated by 365.25 (year ; Search forward. (+ approx -1 ; Lower bound. (sum 1 y approx (>= date (fixed-from-hebrew (hebrew-date 7 1 y)))))) (start ; Starting month for search for month. (if (< date (fixed-from-hebrew (hebrew-date 1 1 year))) 7 ; Tishri 1)) ; Nisan (month ; Search forward from either Tishri or Nisan. (+ start (sum 1 m start (> date (fixed-from-hebrew (hebrew-date m (last-day-of-hebrew-month m year) year)))))) (day ; Calculate the day by subtraction. (1+ (- date (fixed-from-hebrew (hebrew-date month 1 year)))))) (hebrew-date month day year))) **) Hebrew[date_Integer] := Module[{approx, year, start, month, day}, approx = Quotient[date - HebrewEpoch[], 35975351/98496]; year = approx - 1 + MSum[(1)&, approx, (date >= ToFixed[Hebrew[7, 1, #]])&]; start = If[date < ToFixed[Hebrew[1, 1, year]], 7, 1]; month = start + MSum[(1)&, start, (date > ToFixed[Hebrew[#, LastDayOfHebrewMonth[#, year], year]])&]; day = 1 + date - ToFixed[Hebrew[month, 1, year]]; Hebrew[month, day, year] ] (** (defun yom-kippur (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Yom Kippur occurring in Gregorian year. (let* ((hebrew-year (1+ (- g-year (gregorian-year-from-fixed hebrew-epoch))))) (fixed-from-hebrew (hebrew-date 7 10 hebrew-year)))) **) YomKippur[gYear_Integer] := Module[{hebrewYear}, hebrewYear = 1 + gYear - GregorianYearFromFixed[HebrewEpoch[]]; ToFixed[Hebrew[7, 10, hebrewYear]] ] (** (defun passover (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Passover occurring in Gregorian year. (let* ((hebrew-year (- g-year (gregorian-year-from-fixed hebrew-epoch)))) (fixed-from-hebrew (hebrew-date 1 15 hebrew-year)))) **) Passover[gYear_Integer] := Module[{hebrewYear}, hebrewYear = gYear - GregorianYearFromFixed[HebrewEpoch[]]; ToFixed[Hebrew[1, 15, hebrewYear]] ] (** (defun omer (date) ;; TYPE fixed-date -> omer-count ;; Number of elapsed weeks and days in the omer at date. ;; Returns bogus if that date does not fall during the ;; omer. (let* ((c (- date (fixed-from-hebrew (hebrew-date 1 15 (standard-year (hebrew-from-fixed date))))))) (if (<= 1 c 49) (list (quotient c 7) (mod c 7)) bogus))) **) Omer[date_Integer] := Module[{c}, c = date - ToFixed[Hebrew[1, 15, CYear[Hebrew[date]]]]; If[1 <= c <= 49, {Quotient[c, 7], Mod[c, 7]}, Bogus] ] (** (defun purim (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Purim occurring in Gregorian year. (let* ((hebrew-year (- g-year (gregorian-year-from-fixed hebrew-epoch))) (last-month ; Adar or Adar II (last-month-of-hebrew-year hebrew-year))) (fixed-from-hebrew (hebrew-date last-month 14 hebrew-year)))) **) Purim[gYear_Integer] := Module[{hebrewYear, lastMonth}, hebrewYear = gYear - GregorianYearFromFixed[HebrewEpoch[]]; lastMonth = LastMonthOfHebrewYear[hebrewYear]; ToFixed[Hebrew[lastMonth, 14, hebrewYear]] ] (** (defun ta-anith-esther (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Ta'anith Esther occurring in ;; Gregorian year. (let* ((purim-date (purim g-year))) (if ; Purim is on Sunday (= (day-of-week-from-fixed purim-date) sunday) ;; Then prior Thursday (- purim-date 3) ;; Else previous day (1- purim-date)))) **) TaAnithEsther[gYear_Integer] := Module[{purimDate}, purimDate = Purim[gYear]; If[DayOfWeekFromFixed[purimDate] == Sunday[], purimDate - 3, purimDate - 1] ] (** (defun tisha-b-av (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Tisha B'Av occurring in Gregorian year. (let* ((hebrew-year (- g-year (gregorian-year-from-fixed hebrew-epoch))) (ninth-of-av (fixed-from-hebrew (hebrew-date 5 9 hebrew-year)))) (if ; Ninth of Av is Saturday (= (day-of-week-from-fixed ninth-of-av) saturday) ;; Then the next day (1+ ninth-of-av) ninth-of-av))) **) TishaBAv[gYear_Integer] := Module[{hebrewYear, ninthOfAv}, hebrewYear = gYear - GregorianYearFromFixed[HebrewEpoch[]]; ninthOfAv = ToFixed[Hebrew[5, 9, hebrewYear]]; If[DayOfWeekFromFixed[ninthOfAv] == Saturday[], ninthOfAv + 1, ninthOfAv] ] (** (defun birkath-ha-hama (g-year) ;; TYPE gregorian-year -> list-of-fixed-dates ;; List of fixed date of Birkath HaHama occurring in ;; Gregorian year, if it occurs. (let* ((mar26 (julian-in-gregorian march 26 g-year))) (if (and (not (equal mar26 nil)) (= (mod (standard-year (julian-from-fixed (first mar26))) 28) 21)) mar26 nil))) **) BirkathHaHama[gYear_Integer] := Module[{mar26}, mar26 = JulianInGregorian[March[], 26, gYear]; If[mar26 =!= {} && Mod[CYear[Julian[First[mar26]]], 28] == 21, mar26, {}] ] (** (defun sh-ela (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Sh'ela occurring in Gregorian year. (- (first (julian-in-gregorian march 26 (1+ g-year))) 124)) **) ShEla[gYear_Integer] := First[JulianInGregorian[March[], 26, gYear + 1]] - 124 (** (defun yom-ha-zikaron (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Yom HaZikaron occurring in Gregorian year. (let* ((hebrew-year (- g-year (gregorian-year-from-fixed hebrew-epoch))) (h; Ordinarily Iyar 4 (fixed-from-hebrew (hebrew-date 2 4 hebrew-year)))) (if (< wednesday (day-of-week-from-fixed h)) ;; But prior Wednesday if Iyar 5 is Friday or ;; Saturday (kday-before h wednesday) h))) **) YomHaZikaron[gYear_Integer] := Module[{hebrewYear, h}, hebrewYear = gYear - GregorianYearFromFixed[HebrewEpoch[]]; h = ToFixed[Hebrew[2, 4, hebrewYear]]; If[Wednesday[] < DayOfWeekFromFixed[h], KDayBefore[h, Wednesday[]], h] ] (** (defun hebrew-birthday (birthdate h-year) ;; TYPE (hebrew-date hebrew-year) -> fixed-date ;; Fixed date of the anniversary of Hebrew birthdate ;; occurring in Hebrew year. This function assumes ;; that the function fixed-from-hebrew works for Hebrew ;; dates month, day, year even if the month has fewer than ;; day days--in that case the function returns the ;; (day-1)st day after month 1, year. (let* ((birth-day (standard-day birthdate)) (birth-month (standard-month birthdate)) (birth-year (standard-year birthdate))) (if ; It's Adar in a normal Hebrew year or Adar II ; in a Hebrew leap year, (= birth-month (last-month-of-hebrew-year birth-year)) ;; Then use the same day in last month of Hebrew year. (fixed-from-hebrew (hebrew-date (last-month-of-hebrew-year h-year) birth-day h-year)) ;; Else use the normal anniversary of the birth date, ;; or the corresponding day in years without that date (fixed-from-hebrew (hebrew-date birth-month birth-day h-year))))) **) HebrewBirthday[birthDate_Hebrew, hYear_Integer] := Module[{birthDay, birthMonth, birthYear}, birthDay = CDay[birthDate]; birthMonth = CMonth[birthDate]; birthYear = CYear[birthDate]; If[birthMonth == LastMonthOfHebrewYear[birthYear], ToFixed[Hebrew[LastMonthOfHebrewYear[hYear], birthDay, hYear]], ToFixed[Hebrew[birthMonth, birthDay, hYear]]] ] (** (defun yahrzeit (death-date h-year) ;; TYPE (hebrew-date hebrew-year) -> fixed-date ;; Fixed date of the anniversary of Hebrew death-date ;; occurring in Hebrew year. This function assumes ;; that the function fixed-from-hebrew works for Hebrew ;; dates month, day, year even if the month has fewer than ;; day days--in that case the function returns the ;; (day-1)st day after month 1, year. (let* ((death-day (standard-day death-date)) (death-month (standard-month death-date)) (death-year (standard-year death-date))) (cond ;; If it's Heshvan 30 it depends on the first ;; anniversary; if that was not Heshvan 30, use ;; the day before Kislev 1. ((and (= death-month 8) (= death-day 30) (not (long-heshvan? (1+ death-year)))) (1- (fixed-from-hebrew (hebrew-date 9 1 h-year)))) ;; If it's Kislev 30 it depends on the first ;; anniversary; if that was not Kislev 30, use ;; the day before Teveth 1. ((and (= death-month 9) (= death-day 30) (short-kislev? (1+ death-year))) (1- (fixed-from-hebrew (hebrew-date 10 1 h-year)))) ;; If it's Adar II, use the same day in last ;; month of Hebrew year (Adar or Adar II). ((= death-month 13) (fixed-from-hebrew (hebrew-date (last-month-of-hebrew-year h-year) death-day h-year))) ;; If it's the 30th in Adar I and Hebrew year is not a ;; Hebrew leap year (so Adar has only 29 days), use the ;; last day in Shevat. ((and (= death-day 30) (= death-month 12) (not (hebrew-leap-year? h-year))) (fixed-from-hebrew (hebrew-date 11 30 h-year))) ;; In all other cases, use the normal anniversary of ;; the date of death. (t (fixed-from-hebrew (hebrew-date death-month death-day h-year)))))) **) Yahrzeit[deathDate_Hebrew, hYear_Integer] := Module[{deathDay, deathMonth, deathYear}, deathDay = CDay[deathDate]; deathMonth = CMonth[deathDate]; deathYear = CMonth[deathDate]; Which[ deathMonth == 8 && deathDay == 30 && !LongHeshvanQ[deathYear + 1], 1 - ToFixed[Hebrew[9, 1, hYear]], deathMonth == 9 && deathDay == 30 && ShortKislevQ[deathYear + 1], 1 - ToFixed[Hebrew[10, 1, hYear]], deathMonth == 13, ToFixed[Hebrew[LastMonthOfHebrewYear[hYear], deathDay, hYear]], deathDay == 30 && deathMonth == 12 && !LeapYearQ[Hebrew, hYear], ToFixed[Hebrew[11, 30, hYear]], True, ToFixed[Hebrew[deathMonth, deathDay, hYear]] ] ] (** ;;;; Section: Mayan Calendars **) (** (defconstant mayan-epoch ;; TYPE fixed-date ;; Fixed date of start of the Mayan calendar, according ;; to the Goodman-Martinez-Thompson correlation. ;; That is, August 11, -3113. (fixed-from-jd 584282.5)) **) MayanEpoch[] = FixedFromJD[584282.5] (** (defun fixed-from-mayan-long-count (count) ;; TYPE mayan-long-count-date -> fixed-date ;; Fixed date corresponding to the Mayan long count ;; count, which is a list (baktun katun tun uinal kin). (let* ((baktun (mayan-baktun count)) (katun (mayan-katun count)) (tun (mayan-tun count)) (uinal (mayan-uinal count)) (kin (mayan-kin count))) (+ mayan-epoch ; Fixed date at Mayan 0.0.0.0.0 ( * baktun 144000); Baktun. ( * katun 7200) ; Katun. ( * tun 360) ; Tun. ( * uinal 20) ; Uinal. kin))) ; Kin (days). **) ToFixed[date_MayanLongCount] := MayanEpoch[] + CBaktun[date] 144000 + CKatun[date] 7200 + CTun[date] 360 + CUinal[date] 20 + CKin[date] (** (defun mayan-long-count-from-fixed (date) ;; TYPE fixed-date -> mayan-long-count-date ;; Mayan long count date of fixed date. (let* ((long-count (- date mayan-epoch)) (baktun (quotient long-count 144000)) (day-of-baktun (mod long-count 144000)) (katun (quotient day-of-baktun 7200)) (day-of-katun (mod day-of-baktun 7200)) (tun (quotient day-of-katun 360)) (day-of-tun (mod day-of-katun 360)) (uinal (quotient day-of-tun 20)) (kin (mod day-of-tun 20))) (mayan-long-count-date baktun katun tun uinal kin))) **) MayanLongCount[date_Integer] := Module[{longCount, baktun, dayOfBaktun, katun, dayOfKatun, tun, dayOfTun, uinal, kin}, longCount = date - MayanEpoch[]; baktun = Quotient[longCount, 144000]; dayOfBaktun = Mod[longCount, 144000]; katun = Quotient[dayOfBaktun, 7200]; dayOfKatun = Mod[dayOfBaktun, 7200]; tun = Quotient[dayOfKatun, 360]; dayOfTun = Mod[dayOfKatun, 360]; uinal = Quotient[dayOfTun, 20]; kin = Mod[dayOfTun, 20]; MayanLongCount[baktun, katun, tun, uinal, kin] ] (** (defconstant mayan-haab-at-epoch ;; TYPE mayan-haab-date ;; Haab date at long count 0.0.0.0.0. (mayan-haab-date 8 18)) **) MayanHaabAtEpoch[] = MayanHaab[8, 18] (** (defun mayan-haab-from-fixed (date) ;; TYPE fixed-date -> mayan-haab-date ;; Mayan haab date of fixed date. (let* ((long-count (- date mayan-epoch)) (day-of-haab (mod (+ long-count (mayan-haab-day mayan-haab-at-epoch) ( * 20 (1- (mayan-haab-month mayan-haab-at-epoch)))) 365)) (day (mod day-of-haab 20)) (month (1+ (quotient day-of-haab 20)))) (mayan-haab-date day month))) **) MayanHaab[date_Integer] := Module[{longCount, dayOfHaab, day, month}, longCount = date - MayanEpoch[]; dayOfHaab = Mod[longCount + CDay[MayanHaabAtEpoch[]] + 20 (CMonth[MayanHaabAtEpoch[]] - 1), 365]; day = Mod[dayOfHaab, 20]; month = 1 + Quotient[dayOfHaab, 20]; MayanHaab[day, month] ] (** (defun mayan-haab-difference (h-date1 h-date2) ;; TYPE (mayan-haab-date mayan-haab-date) -> integer ;; Number of days from Mayan haab date h-date1 to the next ;; occurrence of Mayan haab date h-date2. (let* ((day1 (mayan-haab-day h-date1)) (day2 (mayan-haab-day h-date2)) (month1 (mayan-haab-month h-date1)) (month2 (mayan-haab-month h-date2))) (mod (+ ( * 20 (- month2 month1)) (- day2 day1)) 365))) **) MayanHaabDifference[hDate1_MayanHaab, hDate2_MayanHaab] := Mod[20 (CMonth[hDate2] - CMonth[hDate1]) + (CDay[hDate2] - CDay[hDate1]), 365] (** (defun mayan-haab-on-or-before (haab date) ;; TYPE (mayan-haab-date fixed-date) -> fixed-date ;; Fixed date of latest date on or before fixed date ;; that is Mayan haab date haab. (- date (mod (- date (mayan-haab-difference (mayan-haab-from-fixed 0) haab)) 365))) **) MayanHaabOnOrBefore[haab_MayanHaab, date_Integer] := date - Mod[date - MayanHaabDifference[MayanHaab[0], haab], 365] (** (defconstant mayan-tzolkin-at-epoch ;; TYPE mayan-tzolkin-date ;; Tzolkin date at long count 0.0.0.0.0. (mayan-tzolkin-date 4 20)) **) MayanTzolkinAtEpoch[] = MayanTzolkin[4, 20] (** (defun mayan-tzolkin-from-fixed (date) ;; TYPE fixed-date -> mayan-tzolkin-date ;; Mayan tzolkin date of fixed date. (let* ((long-count (- date mayan-epoch)) (number (adjusted-mod (+ long-count (mayan-tzolkin-number mayan-tzolkin-at-epoch)) 13)) (name (adjusted-mod (+ long-count (mayan-tzolkin-name mayan-tzolkin-at-epoch)) 20))) (mayan-tzolkin-date number name))) **) MayanTzolkin[date_Integer] := Module[{longCount, number, name}, longCount = date - MayanEpoch[]; number = AdjustedMod[longCount + CNumber[MayanTzolkinAtEpoch[]], 13]; name = AdjustedMod[longCount + CName[MayanTzolkinAtEpoch[]], 20]; MayanTzolkin[number, name] ] (** (defun mayan-tzolkin-difference (t-date1 t-date2) ;; TYPE (mayan-tzolkin-date mayan-tzolkin-date) -> integer ;; Number of days from Mayan tzolkin date t-date1 to the ;; next occurrence of Mayan tzolkin date t-date2. (let* ((number1 (mayan-tzolkin-number t-date1)) (number2 (mayan-tzolkin-number t-date2)) (name1 (mayan-tzolkin-name t-date1)) (name2 (mayan-tzolkin-name t-date2)) (number-difference (- number2 number1)) (name-difference (- name2 name1))) (mod (+ number-difference ( * 13 (mod ( * 3 (- number-difference name-difference)) 20))) 260))) **) MayanTzolkinDifference[date1_MayanTzolkin, date2_MayanTzolkin] := Module[{numberDifference, nameDifference}, numberDifference = CNumber[date2] - CNumber[date1]; nameDifference = CName[date2] - CName[date1]; Mod[numberDifference + 13 Mod[3 (numberDifference - nameDifference), 20], 260] ] (** (defun mayan-tzolkin-on-or-before (tzolkin date) ;; TYPE (mayan-tzolkin-date fixed-date) -> fixed-date ;; Fixed date of latest date on or before fixed date ;; that is Mayan tzolkin date tzolkin. (- date (mod (- date (mayan-tzolkin-difference (mayan-tzolkin-from-fixed 0) tzolkin)) 260))) **) MayanTzolkinOnOrBefore[tzolkin_MayanTzolkin, date_Integer] := date - Mod[date - MayanTzolkinDifference[MayanTzolkin[0], tzolkin], 260] (** (defun mayan-haab-tzolkin-on-or-before (haab tzolkin date) ;; TYPE (mayan-haab-date mayan-tzolkin-date fixed-date) ;; TYPE -> fixed-date ;; Fixed date of latest date on or before date ;; that is Mayan haab date haab and tzolkin date tzolkin; ;; returns bogus if such a haab-tzolkin combination is ;; impossible. (let* ((haab-difference (mayan-haab-difference (mayan-haab-from-fixed 0) haab)) (tzolkin-difference (mayan-tzolkin-difference (mayan-tzolkin-from-fixed 0) tzolkin)) (diff (- tzolkin-difference haab-difference))) (if (= (mod diff 5) 0) (- date (mod (- date (+ haab-difference ( * 365 diff))) 18980)) bogus))); haab-tzolkin combination is impossible. **) MayanHaabTzolkinOnOrBefore[haab_MayanHaab, tzolkin_MayanTzolkin, date_Integer] := Module[{haabDifference, tzolkinDifference, diff}, haabDifference = MayanHaabDifference[MayanHaab[0], haab]; tzolkinDifference = MayanTzolkinDifference[MayanTzolkin[0], tzolkin]; diff = tzolkinDifference - haabDifference; If[Mod[diff, 5] == 0, date - Mod[date - (haabDifference + 365 diff), 18980], Bogus] ] (** ;;;; Section: Old Hindu Calendars **) (** (defconstant hindu-epoch ;; TYPE fixed-date ;; Fixed date of start of the Hindu calendar (Kali Yuga). (fixed-from-julian (julian-date february 18 (bce 3102)))) **) HinduEpoch[] = ToFixed[Julian[February[], 18, BCE[3102]]] (** (defun hindu-day-count (date) ;; TYPE fixed-date -> integer ;; Elapsed days (Ahargana) to date since Hindu epoch (K.Y.). (- date hindu-epoch)) **) HinduDayCount[date_] := date - HinduEpoch[] (** (defconstant arya-sidereal-year ;; TYPE rational ;; Length of Old Hindu solar year. 1577917500/4320000) **) AryaSiderealYear[] = 1577917500/4320000 (** (defconstant arya-solar-month ;; TYPE rational ;; Length of Old Hindu solar month. (/ arya-sidereal-year 12)) **) AryaSolarMonth[] = AryaSiderealYear[] / 12 (** (defun old-hindu-solar-from-fixed (date) ;; TYPE fixed-date -> hindu-solar-date ;; Old Hindu solar date equivalent to fixed date. (let* ((rise ; Sunrise on Hindu date. (+ (hindu-day-count date) 1/4)) (year ; Elapsed years. (quotient rise arya-sidereal-year)) (month (1+ (mod (quotient rise arya-solar-month) 12))) (day (1+ (floor (mod rise arya-solar-month))))) (hindu-solar-date month day year))) **) OldHinduSolar[date_Integer] := Module[{rise, year, month, day}, rise = HinduDayCount[date] + 1/4; year = Quotient[rise, AryaSiderealYear[]]; month = 1 + Mod[Quotient[rise, AryaSolarMonth[]], 12]; day = 1 + Floor[Mod[rise, AryaSolarMonth[]]]; OldHinduSolar[month, day, year] ] (** (defun fixed-from-old-hindu-solar (s-date) ;; TYPE hindu-solar-date -> fixed-date ;; Fixed date corresponding to Old Hindu solar date. (let* ((month (standard-month s-date)) (day (standard-day s-date)) (year (standard-year s-date))) (floor (+ hindu-epoch ; Since start of era. ( * year arya-sidereal-year) ; Days in elapsed years ( * (1- month) arya-solar-month) ; ...in months. day -1/4)))) ; Midnight of day. **) ToFixed[date_OldHinduSolar] := Floor[ HinduEpoch[] + CYear[date] AryaSiderealYear[] + (CMonth[date] - 1) AryaSolarMonth[] + CDay[date] - 1/4 ] (** (defconstant arya-lunar-month ;; TYPE rational ;; Length of Old Hindu lunar month. 1577917500/53433336) **) AryaLunarMonth[] = 1577917500/53433336 (** (defconstant arya-lunar-day ;; TYPE rational ;; Length of Old Hindu lunar day. (/ arya-lunar-month 30)) **) AryaLunarDay[] = AryaLunarMonth[] / 30 (** (defun old-hindu-lunar-leap-year? (l-year) ;; TYPE old-hindu-lunar-year -> boolean ;; True if year is a leap year on the ;; old Hindu calendar. (>= (mod (- ( * l-year arya-sidereal-year) arya-solar-month) arya-lunar-month) (- arya-lunar-month (mod arya-sidereal-year arya-lunar-month)))) **) LeapYearQ[OldHinduLunar, lYear_] := Mod[lYear AryaSiderealYear[] - AryaSolarMonth[], AryaLunarMonth[]] >= AryaLunarMonth[] - Mod[AryaSiderealYear[], AryaLunarMonth[]] (** (defun old-hindu-lunar-from-fixed (date) ;; TYPE fixed-date -> old-hindu-lunar-date ;; Old Hindu lunar date equivalent to fixed date. (let* ((rise ; Sunrise on Hindu date. (+ (hindu-day-count date) 1/4)) (new-moon ; Beginning of lunar month. (- rise (mod rise arya-lunar-month))) (leap ; If lunar contained in solar. (and (>= (- arya-solar-month arya-lunar-month) (mod new-moon arya-solar-month)) (> (mod new-moon arya-solar-month) 0))) (month ; Next solar month's name. (1+ (mod (ceiling (/ new-moon arya-solar-month)) 12))) (day ; Lunar days since beginning of lunar month. (1+ (mod (quotient rise arya-lunar-day) 30))) (year ; Solar year at end of lunar month(s). (1- (ceiling (/ (+ new-moon arya-solar-month) arya-sidereal-year))))) (old-hindu-lunar-date month leap day year))) **) OldHinduLunar[date_Integer] := Module[{rise, newMoon, leap, month, day, year}, rise = HinduDayCount[date] + 1/4; newMoon = rise - Mod[rise, AryaLunarMonth[]]; leap = AryaSolarMonth[] - AryaLunarMonth[] >= Mod[newMoon, AryaSolarMonth[]] && Mod[newMoon, AryaSolarMonth[]] > 0; month = 1 + Mod[Ceiling[newMoon / AryaSolarMonth[]], 12]; day = 1 + Mod[Quotient[rise, AryaLunarDay[]], 30]; year = Ceiling[(newMoon + AryaSolarMonth[]) / AryaSiderealYear[]] - 1; OldHinduLunar[month, leap, day, year] ] (** (defun fixed-from-old-hindu-lunar (l-date) ;; TYPE old-hindu-lunar-date -> fixed-date ;; Fixed date corresponding to Old Hindu lunar date. (let* ((year (old-hindu-lunar-year l-date)) (month (old-hindu-lunar-month l-date)) (leap (old-hindu-lunar-leap l-date)) (day (old-hindu-lunar-day l-date)) (mina ; One solar month before solar new year. ( * (1- ( * 12 year)) arya-solar-month)) (lunar-new-year ; New moon after mina. ( * arya-lunar-month (1+ (quotient mina arya-lunar-month))))) (floor (+ hindu-epoch lunar-new-year ( * arya-lunar-month (if ; If there was a leap month this year. (and (not leap) (<= (ceiling (/ (- lunar-new-year mina) (- arya-solar-month arya-lunar-month))) month)) month (1- month))) ( * (1- day) arya-lunar-day) ; Lunar days. 3/4)))) ; Add one if lunar day begins after sunrise. **) ToFixed[date_OldHinduLunar] := Module[{month, mina, lunarNewYear}, month = CMonth[date]; mina = (12 CYear[date] - 1) AryaSolarMonth[]; lunarNewYear = AryaLunarMonth[] (Quotient[mina, AryaLunarMonth[]] + 1); Floor[ HinduEpoch[] + lunarNewYear + AryaLunarMonth[] If[!CLeap[date] && Ceiling[(lunarNewYear - mina) / (AryaSolarMonth[] - AryaLunarMonth[])] <= month, month, month - 1] + (CDay[date] - 1) AryaLunarDay[] + 3/4 ] ] (** (defconstant arya-jovian-period ;; TYPE rational ;; Number of days in one revolution of Jupiter around the ;; Sun. 1577917500/364224) **) AryaJovianPeriod[] = 1577917500/364224 (** (defun jovian-year (date) ;; TYPE fixed-date -> {1-60} ;; Year of Jupiter cycle at fixed date. (1+ (mod (quotient (hindu-day-count date) (/ arya-jovian-period 12)) 60))) **) JovianYear[date_Integer] := Mod[Quotient[HinduDayCount[date], AryaJovianPeriod[] / 12], 60] + 1 (** ;;;; Section: Time and Astronomy **) (** (defun time-from-moment (moment) ;; TYPE moment -> (hour minute second) ;; Time of day (hour minute second) from moment moment. (let* ((hour (floor (mod ( * moment 24) 24))) (minute (floor (mod ( * moment 24 60) 60))) (second (double-float (mod ( * moment 24 60 60) 60)))) (time-of-day hour minute second))) **) TimeOfDay[moment_] := Module[{hour, minute, second}, hour = Floor[Mod[moment 24, 24]]; minute = Floor[Mod[moment 24 60, 60]]; second = Mod[moment 24 60 60, 60]; TimeOfDay[hour, minute, second] ] (* ADDED *) ToMoment[time_TimeOfDay] := CHour[time] / 24 + CMinute[time] / (24 60) + CSecond[time] / (24 60 60) (** (defun local-from-universal (u-time zone) ;; TYPE (moment minute) -> moment ;; Local time from u-time in universal time at time-zone ;; zone. (+ u-time (/ zone 24d0 60d0))) **) LocalFromUniversal[uTime_, zone_] := uTime + zone / (24 60) (** (defun universal-from-local (l-time zone) ;; TYPE (moment minute) -> moment ;; Universal time from l-time in local time at time-zone ;; zone. (- l-time (/ zone 24d0 60d0))) **) UniversalFromLocal[lTime_, zone_] := lTime - zone / (24 60) (** (defun location-offset (longitude zone) ;; TYPE (angle minute) -> minute ;; Offset of location at longitude ;; from standard time at zone. (- ( * 4 longitude) zone)) **) LocationOffset[longitude_, zone_] := 4 longitude - zone (** (defun local-from-standard (s-time offset) ;; TYPE (moment minute) -> moment ;; Local time from standard s-time at distance ;; offset (in minutes) from time zone. (+ s-time (/ offset 24d0 60d0))) **) LocalFromStandard[sTime_, offset_] := sTime + offset / (24 60) (** (defun standard-from-local (l-time offset) ;; TYPE (moment minute) -> moment ;; Standard time from local l-time at distance ;; offset (in minutes) from time zone. (- l-time (/ offset 24d0 60d0))) **) StandardFromLocal[lTime_, offset_] := lTime - offset / (24 60) (** (defun degrees (theta) ;; TYPE real -> angle ;; Normalize angle theta to range 0-360 degrees. (mod theta 360)) **) Degrees[theta_] := Mod[theta, 360] (** (defun radians-to-degrees (theta) ;; TYPE radian -> angle ;; Convert angle theta from radians to degrees. (degrees (/ theta pi 1/180))) **) RadiansToDegrees[theta_] := Degrees[theta / Pi 180] (** (defun degrees-to-radians (theta) ;; TYPE real -> radian ;; Convert angle theta from degrees to radians. ( * (degrees theta) pi 1/180)) **) DegreesToRadians[theta_] := Degrees[theta] Pi / 180 (** (defun sin-degrees (theta) ;; TYPE angle -> amplitude ;; Sine of theta (given in degrees). (sin (degrees-to-radians theta))) **) SinDegrees[theta_] := Sin[DegreesToRadians[theta]] (** (defun cosine-degrees (theta) ;; TYPE angle -> amplitude ;; Cosine of theta (given in degrees). (cos (degrees-to-radians theta))) **) CosineDegrees[theta_] := Cos[DegreesToRadians[theta]] (** (defun tangent-degrees (theta) ;; TYPE angle -> real ;; Tangent of theta (given in degrees). (tan (degrees-to-radians theta))) **) TangentDegrees[theta_] := Tan[DegreesToRadians[theta]] (** (defun arctan-degrees (x quad) ;; TYPE (real quadrant) -> angle ;; Arctangent of x in degrees in quadrant quad. (let* ((deg (radians-to-degrees (atan x)))) (if (or (= quad 1) (= quad 4)) deg (+ deg 180)))) **) ArcTanDegrees[x_, quad_] := Module[{deg}, deg = RadiansToDegrees[ArcTan[x]]; If[quad == 1 || quad == 4, deg, deg + 180] ] (** (defun arcsin-degrees (x) ;; TYPE amplitude -> angle ;; Arcsine of x in degrees. (radians-to-degrees (asin x))) **) ArcSinDegrees[x_] := RadiansToDegrees[ArcSin[x]] (** (defun arccos-degrees (x) ;; TYPE amplitude -> angle ;; Arccosine of x in degrees. (radians-to-degrees (acos x))) **) ArcCosDegrees[x_] := RadiansToDegrees[ArcCos[x]] (** (defun local-from-apparent (moment) ;; TYPE moment -> moment ;; Local time from sundial time. (- moment (equation-of-time moment))) **) LocalFromApparent[moment_] := moment - EquationOfTime[moment] (** (defun apparent-from-local (moment) ;; TYPE moment -> moment ;; Sundial time at local time. (+ moment (equation-of-time moment))) **) ApparentFromLocal[moment_] := moment + EquationOfTime[moment] (** (defun solar-moment (date latitude longitude rise-or-set) ;; TYPE (fixed-date angle angle real) -> moment ;; Local time (fraction of day) of sunrise/sunset at ;; latitude, longitude (in nonpolar regions) for fixed ;; date. rise-or-set is -0.25 for sunrise and +0.25 for ;; sunset. (let* ((approx ; Approximate time of event. (+ (day-number (gregorian-from-fixed date)) 0.5 rise-or-set (/ longitude -360.0d0))) (anomaly ; Anomaly of sun. (- ( * 0.9856d0 approx) 3.289d0)) (sun ; Longitude of sun. (degrees (+ anomaly ( * 1.916d0 (sin-degrees anomaly)) 282.634d0 ( * 0.020d0 (sin-degrees ( * 2 anomaly)))))) (right-ascension ; Right ascension of sun. (arctan-degrees ( * (cosine-degrees 23.441884d0) (tangent-degrees sun)) (1+ (quotient sun 90)))) ; Quadrant. (declination ; Declination of sun. (arcsin-degrees ( * (sin-degrees 23.441884d0) (sin-degrees sun)))) (local ( * (signum rise-or-set) (arccos-degrees (/ (- (cosine-degrees 90.833333d0) ( * (sin-degrees declination) (sin-degrees latitude))) (cosine-degrees declination) (cosine-degrees latitude)))))) (mod (- (/ (+ local right-ascension) 360) 0.27592d0 ( * 0.00273792d0 approx)) 1))) **) SolarMoment[date_Integer, latitude_, longitude_, riseOrSet_] := Module[{approx, anomaly, sun, rightAscension, declination, local}, approx = DayNumber[Gregorian[date]] + 0.5 + riseOrSet + (longitude / -360.0); anomaly = 0.9856 approx - 3.289; sun = Degrees[anomaly + 1.916 SinDegrees[anomaly] + 282.634 + 0.020 SinDegrees[2 anomaly]]; rightAscension = ArcTanDegrees[CosineDegrees[23.441884] TangentDegrees[sun], 1 + Quotient[sun, 90]]; declination = ArcSinDegrees[SinDegrees[23.441884] SinDegrees[sun]]; local = Sign[riseOrSet] ArcCosDegrees[ (CosineDegrees[90.833333] - SinDegrees[declination] SinDegrees[latitude]) / CosineDegrees[declination] / CosineDegrees[latitude]]; N[Mod[(local + rightAscension) / 360 - 0.27592 - 0.00273792 approx, 1]] ] (** (defun sunrise (date latitude longitude) ;; TYPE (fixed-date angle angle) -> moment ;; Local time (fraction of day) of sunrise at latitude, ;; longitude (in nonpolar regions) for fixed date. (solar-moment date latitude longitude -0.25d0)) **) Sunrise[date_Integer, latitude_, longitude_] := SolarMoment[date, latitude, longitude, -0.25] (** (defun sunset (date latitude longitude) ;; TYPE (fixed-date angle angle) -> moment ;; Local time (fraction of day) of sunset at latitude, ;; longitude (in nonpolar regions) for fixed date. (solar-moment date latitude longitude 0.25d0)) **) Sunset[date_Integer, latitude_, longitude_] := SolarMoment[date, latitude, longitude, 0.25] (** (defun universal-from-ephemeris (jd) ;; TYPE julian-day-number -> julian-day-number ;; Universal time from Ephemeris time. (- jd (ephemeris-correction (moment-from-jd jd)))) **) UniversalFromEphemeris[jd_] := jd - EphemerisCorrection[MomentFromJD[jd]] (** (defun ephemeris-from-universal (jd) ;; TYPE julian-day-number -> julian-day-number ;; Ephemeris time at Universal time. (+ jd (ephemeris-correction (moment-from-jd jd)))) **) EphemerisFromUniversal[jd_] := jd + EphemerisCorrection[MomentFromJD[jd]] (** (defconstant j2000 ;; TYPE julian-day-number ;; Julian day number (2451545) of Gregorian year 2000. (jd-from-moment (+ 0.5d0 (fixed-from-gregorian (gregorian-date january 1 2000))))) **) J2000[] = JDFromMoment[0.5 + ToFixed[Gregorian[January[], 1, 2000]]] (** (defun julian-centuries (moment) ;; TYPE moment -> moment ;; Julian centuries since j2000 at Universal time. (/ (- (ephemeris-from-universal moment) j2000) 36525d0)) **) JulianCenturies[moment_] := (EphemerisFromUniversal[moment] - J2000[]) / 36525 (** (defconstant mean-tropical-year ;; TYPE real 365.242199d0) **) MeanTropicalYear[] = 365.242199 (** (defconstant mean-synodic-month ;; TYPE real 29.530588853d0) **) MeanSynodicMonth[] = 29.530588853 (** (defun ephemeris-correction (moment) ;; TYPE moment -> fraction-of-day ;; Ephemeris Time minus Universal Time (in days) for ;; fixed time. Adapted from "Astronomical Algorithms" ;; by Jean Meeus, Willmann-Bell, Inc., 1991. (let* ((year (gregorian-year-from-fixed moment)) (theta (/ (gregorian-date-difference (gregorian-date january 1 1900) (gregorian-date july 1 year)) 36525d0)) (coeff-19th (list -0.00002d0 0.000297d0 0.025184d0 -0.181133d0 0.553040d0 -0.861938d0 0.677066d0 -0.212591d0)) (coeff-18th (list -0.000009d0 0.003844d0 0.083563d0 0.865736d0 4.867575d0 15.845535d0 31.332267d0 38.291999d0 28.316289d0 11.636204d0 2.043794d0))) (cond ((<= 1988 year 2019) (/ (- year 1933) 24d0 60d0 60d0)) ((<= 1900 year 1987) (poly theta coeff-19th)) ((<= 1800 year 1899) (poly theta coeff-18th)) ((<= 1620 year 1799) (/ (poly (- year 1600) (list 196.58333d0 -4.0675d0 0.0219167d0)) 24d0 60d0 60d0)) (t (let* ((x (+ 0.5d0 (gregorian-date-difference (gregorian-date january 1 1810) (gregorian-date january 1 year))))) (/ (- (/ ( * x x) 41048480d0) 15) 24d0 60d0 60d0)))))) **) EphemerisCorrection[moment_] := Module[{year, theta, coeff19th, coeff18th}, year = GregorianYearFromFixed[moment]; theta = GregorianDateDifference[Gregorian[January[], 1, 1900], Gregorian[July[], 1, year]] / 36525; coeff19th = {-0.00002, 0.000297, 0.025184, -0.181133, 0.553040, -0.861938, 0.677066, -0.212591}; coeff18th = {-0.000009, 0.003844, 0.083563, 0.865736, 4.867575, 15.845535, 31.332267, 38.291999, 28.316289, 11.636204, 2.043794}; Which[ 1988 <= year <= 2019, (year - 1933) / (24 60 60), 1900 <= year <= 1987, Poly[theta, coeff19th], 1800 <= year <= 1899, Poly[theta, coeff18th], 1620 <= year <= 1799, Poly[year - 1600, {196.58333, -4.0675, 0.0219167}] / (24 60 60), True, x = 0.5 + GregorianDateDifference[Gregorian[January[], 1, 1810], Gregorian[January[], 1, year]]; (x^2 / 41048480 - 15) / (24 60 60)] ] (** (defun equation-of-time (jd) ;; TYPE moment -> fraction-of-day ;; Equation of time (in days) for julian day number jd. ;; Adapted from "Astronomical Algorithms" by Jean Meeus, ;; Willmann-Bell, Inc., 1991. (let* ((c (/ (- jd j2000) 36525d0)) (longitude (poly c (list 280.46645d0 36000.76983d0 0.0003032d0))) (anomaly (poly c (list 357.52910d0 35999.05030d0 -0.0001559d0 -0.00000048d0))) (inclination (poly c (list 23.43929111d0 -0.013004167d0 -0.00000016389d0 0.0000005036d0))) (eccentricity (poly c (list 0.016708617d0 -0.000042037d0 -0.0000001236d0))) (y (expt (tangent-degrees (/ inclination 2)) 2))) (/ (+ ( * y (sin-degrees ( * 2 longitude))) ( * -2 eccentricity (sin-degrees anomaly)) ( * 4 eccentricity y (sin-degrees anomaly) (cosine-degrees ( * 2 longitude))) ( * -0.5 y y (sin-degrees ( * 4 longitude))) ( * -1.25 eccentricity eccentricity (sin-degrees ( * 2 anomaly)))) 2 pi))) **) EquationOfTime[jd_] := Module[{c, longitude, anomaly, inclination, eccentricity, y}, c = (jd - J2000[]) / 36525; longitude = Poly[c, {280.46645, 36000.76983, 0.0003032}]; anomaly = Poly[c, {357.52910, 35999.05030, -0.0001559, -0.00000048}]; inclination = Poly[c, {23.43929111, -0.013004167, -0.00000016389, 0.0000005036}]; eccentricity = Poly[c, {0.016708617, -0.000042037, -0.0000001236}]; y = TangentDegrees[inclination / 2]^2; N[(y SinDegrees[2 longitude] + -2 eccentricity SinDegrees[anomaly] + 4 eccentricity y SinDegrees[anomaly] CosineDegrees[2 longitude] + -0.5 y^2 SinDegrees[4 longitude] + -1.25 eccentricity^2 SinDegrees[2 anomaly]) / (2 Pi)] ] (** (defun solar-longitude (jd) ;; TYPE julian-day-number -> angle ;; Longitude of sun on astronomical (julian) day number jd. ;; Adapted from "Planetary Programs and Tables from -4000 ;; to +2800" by Pierre Bretagnon and Jean-Louis Simon, ;; Willmann-Bell, Inc., 1986. (let* ((c ; Ephemeris time in Julian centuries (julian-centuries jd)) (coefficients (list 403406 195207 119433 112392 3891 2819 1721 0 660 350 334 314 268 242 234 158 132 129 114 99 93 86 78 72 68 64 46 38 37 32 29 28 27 27 25 24 21 21 20 18 17 14 13 13 13 12 10 10 10 10)) (multipliers (list 0.01621043d0 628.30348067d0 628.30821524d0 628.29634302d0 1256.605691d0 1256.60984d0 628.324766d0 0.00813d0 1256.5931d0 575.3385d0 -0.33931d0 7771.37715d0 786.04191d0 0.05412d0 393.02098d0 -0.34861d0 1150.67698d0 157.74337d0 52.9667d0 588.4927d0 52.9611d0 -39.807d0 522.3769d0 550.7647d0 2.6108d0 157.7385d0 1884.9103d0 -77.5655d0 2.6489d0 1179.0627d0 550.7575d0 -79.6139d0 1884.8981d0 21.3219d0 1097.7103d0 548.6856d0 254.4393d0 -557.3143d0 606.9774d0 21.3279d0 1097.7163d0 -77.5282d0 1884.9191d0 2.0781d0 294.2463d0 -0.0799d0 469.4114d0 -0.6829d0 214.6325d0 1572.084d0)) (addends (list 4.721964d0 5.937458d0 1.115589d0 5.781616d0 5.5474d0 1.512d0 4.1897d0 1.163d0 5.415d0 4.315d0 4.553d0 5.198d0 5.989d0 2.911d0 1.423d0 0.061d0 2.317d0 3.193d0 2.828d0 0.52d0 4.65d0 4.35d0 2.75d0 4.5d0 3.23d0 1.22d0 0.14d0 3.44d0 4.37d0 1.14d0 2.84d0 5.96d0 5.09d0 1.72d0 2.56d0 1.92d0 0.09d0 5.98d0 4.03d0 4.47d0 0.79d0 4.24d0 2.01d0 2.65d0 4.98d0 0.93d0 2.21d0 3.59d0 1.5d0 2.55d0)) (longitude (+ 4.9353929d0 ( * 628.33196168d0 c) ( * 0.0000001d0 (sigma ((x coefficients) (y addends) (z multipliers)) ( * x (sin (+ y ( * z c))))))))) (radians-to-degrees (+ longitude (aberration c) (nutation c))))) **) SolarLongitude[jd_] := Module[{c, coefficients, multipliers, addends, longitude, x, y, z}, c = JulianCenturies[jd]; coefficients = {403406, 195207, 119433, 112392, 3891, 2819, 1721, 0, 660, 350, 334, 314, 268, 242, 234, 158, 132, 129, 114, 99, 93, 86, 78, 72, 68, 64, 46, 38, 37, 32, 29, 28, 27, 27, 25, 24, 21, 21, 20, 18, 17, 14, 13, 13, 13, 12, 10, 10, 10, 10}; multipliers = {0.01621043, 628.30348067, 628.30821524, 628.29634302, 1256.605691, 1256.60984, 628.324766, 0.00813, 1256.5931, 575.3385, -0.33931, 7771.37715, 786.04191, 0.05412, 393.02098, -0.34861, 1150.67698, 157.74337, 52.9667, 588.4927, 52.9611, -39.807, 522.3769, 550.7647, 2.6108, 157.7385, 1884.9103, -77.5655, 2.6489, 1179.0627, 550.7575, -79.6139, 1884.8981, 21.3219, 1097.7103, 548.6856, 254.4393, -557.3143, 606.9774, 21.3279, 1097.7163, -77.5282, 1884.9191, 2.0781, 294.2463, -0.0799, 469.4114, -0.6829, 214.6325, 1572.084}; addends = {4.721964, 5.937458, 1.115589, 5.781616, 5.5474, 1.512, 4.1897, 1.163, 5.415, 4.315, 4.553, 5.198, 5.989, 2.911, 1.423, 0.061, 2.317, 3.193, 2.828, 0.52, 4.65, 4.35, 2.75, 4.5, 3.23, 1.22, 0.14, 3.44, 4.37, 1.14, 2.84, 5.96, 5.09, 1.72, 2.56, 1.92, 0.09, 5.98, 4.03, 4.47, 0.79, 4.24, 2.01, 2.65, 4.98, 0.93, 2.21, 3.59, 1.5, 2.55}; longitude = 4.9353929 + 628.33196168 c + 0.0000001 Sigma[{{x, coefficients}, {y, addends}, {z, multipliers}}, x Sin[z c + y]]; N[RadiansToDegrees[longitude + Aberration[c] + Nutation[c]]] ] (** (defun nutation (c) ;; TYPE julian-centuries -> radian ;; Longitudinal nutation in radians at c Julian centuries. (let* ((A (poly c (list 124.90d0 -1934.134d0 0.002063d0))) (B (poly c (list 201.11d0 72001.5377d0 0.00057d0)))) (+ ( * -.0000834d0 (sin-degrees A)) ( * -.0000064d0 (sin-degrees B))))) **) Nutation[c_] := Module[{a, b}, a = Poly[c, {124.90, -1934.134, 0.002063}]; b = Poly[c, {201.11, 72001.5377, 0.00057}]; -.0000834 SinDegrees[a] + -.0000064 SinDegrees[b] ] (** (defun aberration (c) ;; TYPE julian-centuries -> radian ;; Aberration in radians at c Julian centuries. (- ( * 0.0000017d0 (cosine-degrees (+ 177.63d0 ( * 35999.01848d0 c)))) 0.0000973d0)) **) Aberration[c_] := 0.0000017 CosineDegrees[177.63 + 35999.01848 c] - 0.0000973; (** (defun date-next-solar-longitude (jd l) ;; TYPE (julian-day-number angle) -> julian-day-number ;; Julian day number of the first date at or after julian ;; day number jd (in Greenwich) when the solar longitude ;; will be a multiple of l degrees; l must be a proper ;; divisor of 360. (let* ((next (double-float (degrees ( * l (ceiling (/ (solar-longitude jd) l))))))) (binary-search start jd end (+ jd ( * (/ l 360) 400d0)) x (if (= next 0); Discontinuity at next=0 ;; Then test for drop in longitude (>= l (solar-longitude x)) ;; Else test if we are past the desired ;; longitude (>= (solar-longitude x) next)) (>= 0.00001d0 (- end start))))) **) DateNextSolarLongitude[jd_, l_] := Module[{next}, next = Degrees[l Ceiling[SolarLongitude[jd] / l]]; BinarySearch[jd, jd + l / 360 * 400, Function[{x}, If[next == 0, l >= SolarLongitude[x], SolarLongitude[x] >= next]], Function[{lo, hi}, 0.00001 >= hi - lo]] ] (** (defun lunar-longitude (u-time) ;; TYPE moment -> angle ;; Longitude of moon (in degrees) at u-time (Universal time). ;; Adapted from "Astronomical Algorithms" by Jean Meeus, ;; Willmann-Bell, Inc., 1991. (let* ((c (julian-centuries u-time)) (mean-moon (degrees (poly c (list 218.3164591d0 481267.88134236d0 -.0013268d0 1/538841 -1/65194000)))) (elongation (degrees (poly c (list 297.8502042d0 445267.1115168d0 -.00163d0 1/545868 -1/113065000)))) (solar-anomaly (degrees (poly c (list 357.5291092d0 35999.0502909d0 -.0001536d0 1/24490000)))) (lunar-anomaly (degrees (poly c (list 134.9634114d0 477198.8676313d0 0.008997d0 1/69699 -1/14712000)))) (moon-from-node (degrees (poly c (list 93.2720993d0 483202.0175273d0 -.0034029d0 -1/3526000 1/863310000)))) (e (poly c (list 1 -0.002516d0 -0.0000074d0))) (args-lunar-elongation (list 0 2 2 0 0 0 2 2 2 2 0 1 0 2 0 0 4 0 4 2 2 1 1 2 2 4 2 0 2 2 1 2 0 0 2 2 2 4 0 3 2 4 0 2 2 2 4 0 4 1 2 0 1 3 4 2 0 1 2 2)) (args-solar-anomaly (list 0 0 0 0 1 0 0 -1 0 -1 1 0 1 0 0 0 0 0 0 1 1 0 1 -1 0 0 0 1 0 -1 0 -2 1 2 -2 0 0 -1 0 0 1 -1 2 2 1 -1 0 0 -1 0 1 0 1 0 0 -1 2 1 0 0)) (args-lunar-anomaly (list 1 -1 0 2 0 0 -2 -1 1 0 -1 0 1 0 1 1 -1 3 -2 -1 0 -1 0 1 2 0 -3 -2 -1 -2 1 0 2 0 -1 1 0 -1 2 -1 1 -2 -1 -1 -2 0 1 4 0 -2 0 2 1 -2 -3 2 1 -1 3 -1)) (args-moon-from-node (list 0 0 0 0 0 2 0 0 0 0 0 0 0 -2 2 -2 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 -2 2 0 2 0 0 0 0 0 0 -2 0 0 0 0 -2 -2 0 0 0 0 0 0 0 -2)) (sine-coefficients (list 6288774 1274027 658314 213618 -185116 -114332 58793 57066 53322 45758 -40923 -34720 -30383 15327 -12528 10980 10675 10034 8548 -7888 -6766 -5163 4987 4036 3994 3861 3665 -2689 -2602 2390 -2348 2236 -2120 -2069 2048 -1773 -1595 1215 -1110 -892 -810 759 -713 -700 691 596 549 537 520 -487 -399 -381 351 -340 330 327 -323 299 294 0)) (longitude ( * 1/1000000 (sigma ((v sine-coefficients) (w args-lunar-elongation) (x args-solar-anomaly) (y args-lunar-anomaly) (z args-moon-from-node)) ( * v (expt e (abs x)) (sin-degrees (+ ( * w elongation) ( * x solar-anomaly) ( * y lunar-anomaly) ( * z moon-from-node))))))) (venus ( * 3958/1000000 (sin-degrees (+ 119.75d0 ( * c 131.849d0))))) (jupiter ( * 318/1000000 (sin-degrees (+ 53.09d0 ( * c 479264.29d0))))) (flat-earth ( * 1962/1000000 (sin-degrees (- mean-moon moon-from-node))))) (degrees (+ mean-moon longitude venus jupiter flat-earth (radians-to-degrees (nutation c)))))) **) LunarLongitude[uTime_] := Module[{c, meanMoon, elongation, solarAnomaly, lunarAnomaly, moonFromNode, e, argsLunarElongation, argsSolarAnomaly, argsLunarAnomaly, argsMoonFromNode, sineCoefficients, longitude, venus, jupiter, flatEarth, v, w, x, y, z}, c = JulianCenturies[uTime]; meanMoon = Degrees[Poly[c, {218.3164591, 481267.88134236, -.0013268, 1/538841, -1/65194000}]]; elongation = Degrees[Poly[c, {297.8502042, 445267.1115168, -.00163, 1/545868, -1/113065000}]]; solarAnomaly = Degrees[Poly[c, {357.5291092, 35999.0502909, -.0001536, 1/24490000}]]; lunarAnomaly = Degrees[Poly[c, {134.9634114, 477198.8676313, 0.008997, 1/69699, -1/14712000}]]; moonFromNode = Degrees[Poly[c, {93.2720993, 483202.0175273, -.0034029, -1/3526000, 1/863310000}]]; e = Poly[c, {1, -0.002516, -0.0000074}]; argsLunarElongation = { 0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 1, 0, 2, 0, 0, 4, 0, 4, 2, 2, 1, 1, 2, 2, 4, 2, 0, 2, 2, 1, 2, 0, 0, 2, 2, 2, 4, 0, 3, 2, 4, 0, 2, 2, 2, 4, 0, 4, 1, 2, 0, 1, 3, 4, 2, 0, 1, 2, 2}; argsSolarAnomaly = { 0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, -2, 1, 2, -2, 0, 0, -1, 0, 0, 1, -1, 2, 2, 1, -1, 0, 0, -1, 0, 1, 0, 1, 0, 0, -1, 2, 1, 0, 0}; argsLunarAnomaly = { 1, -1, 0, 2, 0, 0, -2, -1, 1, 0, -1, 0, 1, 0, 1, 1, -1, 3, -2, -1, 0, -1, 0, 1, 2, 0, -3, -2, -1, -2, 1, 0, 2, 0, -1, 1, 0, -1, 2, -1, 1, -2, -1, -1, -2, 0, 1, 4, 0, -2, 0, 2, 1, -2, -3, 2, 1, -1, 3, -1}; argsMoonFromNode = { 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, -2, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, -2, 2, 0, 2, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0, -2}; sineCoefficients = { 6288774, 1274027, 658314, 213618, -185116, -114332, 58793, 57066, 53322, 45758, -40923, -34720, -30383, 15327, -12528, 10980, 10675, 10034, 8548, -7888, -6766, -5163, 4987, 4036, 3994, 3861, 3665, -2689, -2602, 2390, -2348, 2236, -2120, -2069, 2048, -1773, -1595, 1215, -1110, -892, -810, 759, -713, -700, 691, 596, 549, 537, 520, -487, -399, -381, 351, -340, 330, 327, -323, 299, 294, 0}; longitude = 1/1000000 Sigma[{ {v, sineCoefficients}, {w, argsLunarElongation}, {x, argsSolarAnomaly}, {y, argsLunarAnomaly}, {z, argsMoonFromNode}}, v e^Abs[x] SinDegrees[ w elongation + x solarAnomaly + y lunarAnomaly + z moonFromNode] ]; venus = 3958/1000000 SinDegrees[119.75 + c 131.849]; jupiter = 318/1000000 SinDegrees[53.09 + c 479264.29]; flatEarth = 1962/1000000 SinDegrees[meanMoon - moonFromNode]; N[Degrees[meanMoon + longitude + venus + jupiter + flatEarth + RadiansToDegrees[Nutation[c]]]] ] (** (defun new-moon-time (k) ;; TYPE integer -> julian-day-number ;; Astronomical (julian) day number (at Greenwich) of k-th ;; new moon after (or before) the new moon of January 6, ;; 2000. Adapted from "Astronomical Algorithms" by Jean ;; Meeus, Willmann-Bell, Inc., 1991. (let* ((c (/ k 1236.85d0)) (JDE (poly c (list 2451550.09765d0 ( * mean-synodic-month 1236.85d0) 0.0001337d0 -0.000000150d0 0.00000000073d0))) (e (poly c (list 1 -0.002516d0 -0.0000074d0))) (solar-anomaly (poly c (list 2.5534d0 ( * 29.10535669d0 1236.85d0) -0.0000218d0 -0.00000011d0))) (lunar-anomaly (poly c (list 201.5643d0 ( * 385.81693528d0 1236.85d0) 0.0107438d0 0.00001239d0 -0.000000058d0))) (moon-argument (poly c (list 160.7108d0 ( * 390.67050274d0 1236.85d0) -0.0016341d0 -0.00000227d0 0.000000011d0))) (omega (poly c (list 124.7746d0 ( * -1.56375580d0 1236.85d0) 0.0020691d0 0.00000215d0))) (e-factor (list 0 1 0 0 1 1 2 0 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0)) (solar-coeff (list 0 1 0 0 -1 1 2 0 0 1 0 1 1 -1 2 0 3 1 0 1 -1 -1 1 0)) (lunar-coeff (list 1 0 2 0 1 1 0 1 1 2 3 0 0 2 1 2 0 1 2 1 1 1 3 4)) (moon-coeff (list 0 0 0 2 0 0 0 -2 2 0 0 2 -2 0 0 -2 0 -2 2 2 2 -2 0 0)) (sine-coeff (list -0.40720d0 0.17241d0 0.01608d0 0.01039d0 0.00739d0 -0.00514d0 0.00208d0 -0.00111d0 -0.00057d0 0.00056d0 -0.00042d0 0.00042d0 0.00038d0 -0.00024d0 -0.00007d0 0.00004d0 0.00004d0 0.00003d0 0.00003d0 -0.00003d0 0.00003d0 -0.00002d0 -0.00002d0 0.00002d0)) (correction (+ ( * -.00017d0 (sin-degrees omega)) (sigma ((v sine-coeff) (w e-factor) (x solar-coeff) (y lunar-coeff) (z moon-coeff)) ( * v (expt e w) (sin-degrees (+ ( * x solar-anomaly) ( * y lunar-anomaly) ( * z moon-argument))))))) (add-const (list 299.77d0 251.88d0 251.83d0 349.42d0 84.66d0 141.74d0 207.14d0 154.84d0 34.52d0 207.19d0 291.34d0 161.72d0 239.56d0 331.55d0)) (add-coeff (list 0.107408d0 0.016321d0 26.641886d0 36.412478d0 18.206239d0 53.303771d0 2.453732d0 7.306860d0 27.261239d0 0.121824d0 1.844379d0 24.198154d0 25.513099d0 3.592518d0)) (add-extra (list -0.009173d0 0 0 0 0 0 0 0 0 0 0 0 0 0)) (add-factor (list 0.000325d0 0.000165d0 0.000164d0 0.000126d0 0.000110d0 0.000062d0 0.000060d0 0.000056d0 0.000047d0 0.000042d0 0.000040d0 0.000037d0 0.000035d0 0.000023d0)) (additional (sigma ((i add-const) (j add-coeff) (n add-extra) (l add-factor)) ( * l (sin-degrees (+ i ( * j k) ( * n c c))))))) (universal-from-ephemeris (+ JDE correction additional)))) **) NewMoonTime[k_] := Module[{c, JDE, e, solarAnomaly, lunarAnomaly, moonArgument, omega, eFactor, solarCoeff, lunarCoeff, moonCoeff, sineCoeff, correction, addConst, addCoeff, addExtra, addFactor, additional, v, w, x, y, z, i, j, n, l}, c = k / 1236.85; JDE = Poly[c, {2451550.09765, MeanSynodicMonth[] * 1236.85, 0.0001337, -0.000000150, 0.00000000073}]; e = Poly[c, {1, -0.002516, -0.0000074}]; solarAnomaly = Poly[c, {2.5534, 29.10535669 * 1236.85, -0.0000218, -0.00000011}]; lunarAnomaly = Poly[c, {201.5643, 385.81693528 * 1236.85, 0.0107438, 0.00001239, -0.000000058}]; moonArgument = Poly[c, {160.7108, 390.67050274 * 1236.85, -0.0016341, -0.00000227, 0.000000011}]; omega = Poly[c, {124.7746, -1.56375580 * 1236.85, 0.0020691, 0.00000215}]; eFactor = {0, 1, 0, 0, 1, 1, 2, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; solarCoeff = {0, 1, 0, 0, -1, 1, 2, 0, 0, 1, 0, 1, 1, -1, 2, 0, 3, 1, 0, 1, -1, -1, 1, 0}; lunarCoeff = {1, 0, 2, 0, 1, 1, 0, 1, 1, 2, 3, 0, 0, 2, 1, 2, 0, 1, 2, 1, 1, 1, 3, 4}; moonCoeff = {0, 0, 0, 2, 0, 0, 0, -2, 2, 0, 0, 2, -2, 0, 0, -2, 0, -2, 2, 2, 2, -2, 0, 0}; sineCoeff = {-0.40720, 0.17241, 0.01608, 0.01039, 0.00739, -0.00514, 0.00208, -0.00111, -0.00057, 0.00056, -0.00042, 0.00042, 0.00038, -0.00024, -0.00007, 0.00004, 0.00004, 0.00003, 0.00003, -0.00003, 0.00003, -0.00002, -0.00002, 0.00002}; correction = -.00017 SinDegrees[omega] + Sigma[{ {v, sineCoeff}, {w, eFactor}, {x, solarCoeff}, {y, lunarCoeff}, {z, moonCoeff}}, v e^w SinDegrees[x solarAnomaly + y lunarAnomaly + z moonArgument] ]; addConst = {299.77, 251.88, 251.83, 349.42, 84.66, 141.74, 207.14, 154.84, 34.52, 207.19, 291.34, 161.72, 239.56, 331.55}; addCoeff = {0.107408, 0.016321, 26.641886, 36.412478, 18.206239, 53.303771, 2.453732, 7.306860, 27.261239, 0.121824, 1.844379, 24.198154, 25.513099, 3.592518}; addExtra = {-0.009173, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; addFactor = { 0.000325, 0.000165, 0.000164, 0.000126, 0.000110, 0.000062, 0.000060, 0.000056, 0.000047, 0.000042, 0.000040, 0.000037, 0.000035, 0.000023}; additional = Sigma[{ {i, addConst}, {j, addCoeff}, {n, addExtra}, {l, addFactor}}, l SinDegrees[i + j k + n c^2] ]; UniversalFromEphemeris[N[JDE + correction + additional]] ] (** (defun new-moon-at-or-after (jd) ;; TYPE julian-day-number -> julian-day-number ;; Astronomical (julian) day number of first new moon at ;; or after astronomical (julian) day number jd (in ;; Greenwich). The fractional part is the time of day. (let* ((date; Gregorian date, Universal Time (gregorian-from-fixed (floor (moment-from-jd jd)))) (approx; Approximate number of new moons since ; January 6, 2000 (1- (floor ( * (- (+ (standard-year date) (/ (day-number date) 365.25d0)) 2000.0d0) 12.3685d0)))) ; New moons per year. (error (sum 1 k approx (< (new-moon-time k) jd)))) (new-moon-time (+ approx error)))) **) NewMoonAtOrAfter[jd_] := Module[{date, approx, error}, date = Gregorian[Floor[MomentFromJD[jd]]]; approx = Floor[((CYear[date] + DayNumber[date] / 365.25) - 2000.0) * 12.3685] - 1; error = MSum[(1)&, approx, (NewMoonTime[#] < jd)&]; NewMoonTime[approx + error] ] (** (defun new-moon-before (jd) ;; TYPE julian-day-number -> julian-day-number ;; Astronomical (julian) day number of last new moon ;; before astronomical (julian) day number jd (in ;; Greenwich). The fractional part is the time of day. (new-moon-at-or-after (- (new-moon-at-or-after jd) 45))) **) NewMoonBefore[jd_] := NewMoonAtOrAfter[NewMoonAtOrAfter[jd] - 45]; (** (defconstant new ;; TYPE phase ;; Excess of lunar longitude over solar longitude at new ;; moon. 0) **) New[] = 0 (** (defconstant first-quarter ;; TYPE phase ;; Excess of lunar longitude over solar longitude at first ;; quarter moon. 90) **) FirstQuarter[] = 90 (** (defconstant full ;; TYPE phase ;; Excess of lunar longitude over solar longitude at full ;; moon. 180) **) FullMoon[] = 180 (** (defconstant last-quarter ;; TYPE phase ;; Excess of lunar longitude over solar longitude at last ;; quarter moon. 270) **) LastQuarter[] = 270 (** (defun lunar-solar-angle (jd) ;; TYPE julian-day-number -> phase ;; Lunar phase, as an angle in degrees, at astronomical ;; (julian) day number jd. An angle of 0 means a new ;; moon, 90 degrees means the first quarter, 180 means a ;; full moon, and 270 degrees means the last quarter. (degrees (- (lunar-longitude jd) (solar-longitude jd)))) **) LunarSolarAngle[jd_] := Degrees[LunarLongitude[jd] - SolarLongitude[jd]]; (** (defun lunar-phase-at-or-before (phase jd) ;; TYPE (phase julian-day-number) -> julian-day-number ;; Astronomical (julian) day number of last time moon was ;; in phase phase (in degrees) at or before astronomical ;; (julian) day number jd (in Greenwich). The fractional ;; part is the time of day. (let* ((close (< (degrees (- (lunar-solar-angle jd) phase)) 40)) (yesterday (1- jd)) (orig (+ 2451550.26d0 ( * mean-synodic-month (/ phase 360)))) (epsilon 0.000001d0) (tau (- yesterday (mod (- yesterday orig) mean-synodic-month)))) (binary-search l (if close (- jd 4) (- tau 2)) u (if close jd (+ tau 2)) x (<= phase (lunar-solar-angle x) (+ phase 90)) (< (- u l) epsilon)))) **) LunarPhaseAtOrBefore[phase_, jd_] := Module[{close, yesterday, orig, epsilon, tau}, close = Degrees[LunarSolarAngle[jd] - phase] < 40; yesterday = jd - 1; orig = 2451550.26 + MeanSynodicMonth[] * (phase / 360); epsilon = 0.000001; tau = yesterday - Mod[yesterday - orig, MeanSynodicMonth[]]; BinarySearch[If[close, jd - 4, tau - 2], If[close, jd, tau + 2], Function[{x}, phase <= LunarSolarAngle[x] <= phase + 90], Function[{l, u}, u - l < epsilon] ] ] (** (defun new-moon-at-or-before (jd) ;; TYPE julian-day-number -> julian-day-number ;; Astronomical (julian) day number of last new moon ;; at or before astronomical (julian) day number jd (in ;; Greenwich). The fractional part is the time of day. (lunar-phase-at-or-before new jd)) **) NewMoonAtOrBefore[jd_] := LunarPhaseAtOrBefore[New[], jd] (** (defun full-moon-at-or-before (jd) ;; TYPE julian-day-number -> julian-day-number ;; Astronomical (julian) day number of last full moon ;; at or before astronomical (julian) day number jd (in ;; Greenwich). The fractional part is the time of day. (lunar-phase-at-or-before full jd)) **) FullMoonAtOrBefore[jd_] := LunarPhaseAtOrBefore[FullMoon[], jd] (** (defun first-quarter-moon-at-or-before (jd) ;; TYPE julian-day-number -> julian-day-number ;; Astronomical (julian) day number of last first-quarter ;; moon at or before astronomical (julian) day number jd ;; (in Greenwich). The fractional part is the time of ;; day. (lunar-phase-at-or-before first-quarter jd)) **) FirstQuarterMoonAtOrBefore[jd_] := LunarPhaseAtOrBefore[FirstQuarter[], jd] (** (defun last-quarter-moon-at-or-before (jd) ;; TYPE julian-day-number -> julian-day-number ;; Astronomical (julian) day number of last last-quarter ;; moon at or before astronomical (julian) day number jd ;; (in Greenwich). The fractional part is the time of ;; day. (lunar-phase-at-or-before last-quarter jd)) **) LastQuarterMoonAtOrBefore[jd_] := LunarPhaseAtOrBefore[LastQuarter[], jd] (** (defun sidereal-from-jd (jd) ;; TYPE julian-day-number -> julian-day-number ;; Sidereal Time from julian-day-number. ;; Adapted from "Astronomical Algorithms" ;; by Jean Meeus, Willmann-Bell, Inc., 1991. (let* ((c (/ (- jd j2000) 36525)) (sidereal-coeff (list 280.46061837d0 ( * 36525 360.98564736629d0) 0.000387933d0 1/38710000))) (/ (poly c sidereal-coeff) 360))) **) SiderealFromJD[jd_] := Module[{c, siderealCoeff}, c = (jd - J2000[]) / 36525; siderealCoeff = {280.46061837, 36525 * 360.98564736629, 0.000387933, 1/38710000}; Poly[c, siderealCoeff] / 360 ] (** ;;;; Section: French Revolutionary Calendar **) (** (defconstant french-epoch ;; TYPE fixed-date ;; Fixed date of start of the French Revolutionary ;; calendar. (fixed-from-gregorian (gregorian-date september 22 1792))) **) FrenchEpoch[] = ToFixed[Gregorian[September[], 22, 1792]] (** (defconstant french-time-zone ;; TYPE minute ;; The difference (in minutes) of the Paris time zone ;; from Universal Time. (+ 9 21/60)) **) FrenchTimeZone[] = 9 + 21/60 (** (defun french-autumnal-equinox-on-or-before (date) ;; TYPE fixed-date -> fixed-date ;; Fixed date, in Paris local time, of autumnal equinox ;; on or before fixed date. (let* ((theta; solar longitude at end of date in Paris (solar-longitude (universal-from-local (jd-from-moment (1+ date)) french-time-zone))) (d-prime; date shortly before previous autumnal ; equinox (if (< 150d0 theta 180d0); if date is in August or ; in September before the ; autumnal equinox... (- date 370); ...then pick a date before the ; prior autumnal equinox ; ...otherwise, if date is beyond the autumnal ; equinox, pick a date before the last autumnal ; equinox (- date (mod (- date 260) mean-tropical-year))))) (fixed-from-jd (apparent-from-local (local-from-universal (date-next-solar-longitude (universal-from-local (jd-from-moment d-prime) french-time-zone) 90) french-time-zone))))) **) FrenchAutumnalEquinoxOnOrBefore[date_] := Module[{theta, dPrime}, theta = SolarLongitude[UniversalFromLocal[JDFromMoment[date + 1], FrenchTimeZone[]]]; dPrime = If[150 < theta < 180, date - 370, date - Mod[date - 260, MeanTropicalYear[]]]; FixedFromJD[ ApparentFromLocal[ LocalFromUniversal[ DateNextSolarLongitude[ UniversalFromLocal[ JDFromMoment[dPrime], FrenchTimeZone[]], 90], FrenchTimeZone[]] ] ] ] (** (defun fixed-from-french (f-date) ;; TYPE french-date -> fixed-date ;; Fixed date of French Revolutionary date. (let* ((month (standard-month f-date)) (day (standard-day f-date)) (year (standard-year f-date)) (new-year (french-autumnal-equinox-on-or-before (ceiling (+ french-epoch ( * mean-tropical-year (1- year))))))) (+ (1- new-year) ; Days in prior years ( * 30 (1- month)); Days in prior months day))) ; Days this month **) ToFixed[date_French] := Module[{month, day, year, newYear}, month = CMonth[date]; day = CDay[date]; year = CYear[date]; newYear = FrenchAutumnalEquinoxOnOrBefore[ Ceiling[FrenchEpoch[] + MeanTropicalYear[] * (year - 1)]]; (newYear - 1) + 30 (month - 1) + day ] (** (defun french-from-fixed (date) ;; TYPE fixed-date -> french-date ;; French Revolutionary date of fixed date. (let* ((new-year (french-autumnal-equinox-on-or-before date)) (year (1+ (round (/ (- new-year french-epoch) mean-tropical-year)))) (month (1+ (quotient (- date new-year) 30))) (day (1+ (mod (- date new-year) 30)))) (french-date month day year))) **) French[date_Integer] := Module[{newYear, year, month, day}, newYear = FrenchAutumnalEquinoxOnOrBefore[date]; year = Round[(newYear - FrenchEpoch[]) / MeanTropicalYear[]] + 1; month = Quotient[date - newYear, 30] + 1; day = Mod[date - newYear, 30] + 1; French[month, day, year] ] (** (defun modified-french-leap-year? (f-year) ;; TYPE french-year -> boolean ;; True if f-year is a leap year on the French Revolutionary ;; calendar. (and (= (mod f-year 4) 0) (not (member (mod f-year 400) (list 100 200 300))) (not (= (mod f-year 4000) 0)))) **) LeapYearQ[ModifiedFrench, fYear_] := Mod[fYear, 4] == 0 && !MemberQ[{100, 200, 300}, Mod[fYear, 400]] && !Mod[fYear, 4000] == 0 (** (defun fixed-from-modified-french (f-date) ;; TYPE french-date -> fixed-date ;; Fixed date of French Revolutionary date. (let* ((month (first f-date)) (day (second f-date)) (year (third f-date))) (+ french-epoch -1; Days before start of calendar. ( * 365 (1- year)); Ordinary days in prior years. ; Leap days in prior years. (quotient (1- year) 4) (- (quotient (1- year) 100)) (quotient (1- year) 400) (- (quotient (1- year) 4000)) ( * 30 (1- month)); Days in prior months this year. day))); Days this month. **) ToFixed[date_ModifiedFrench] := Module[{year}, year = CYear[date]; FrenchEpoch[] - 1 + 365 (year - 1) + Quotient[year - 1, 4] - Quotient[year - 1, 100] + Quotient[year - 1, 400] - Quotient[year - 1, 4000] + 30 (CMonth[date] - 1) + CDay[date] ] (** (defun modified-french-from-fixed (date) ;; TYPE fixed-date -> french-year ;; French Revolutionary date (month day year) of fixed ;; date. (let* ((approx ; Approximate year (may be off by 1). (quotient (- date french-epoch) 1460969/4000)) (year ; Search forward. (+ approx -1 ; Lower bound (sum 1 y approx (>= date (fixed-from-modified-french (french-date 1 1 y)))))) (month ; Calculate the month by division. (1+ (quotient (- date (fixed-from-modified-french (french-date 1 1 year))) 30))) (day ; Calculate the day by subtraction. (1+ (- date (fixed-from-modified-french (french-date month 1 year)))))) (french-date month day year))) **) ModifiedFrench[date_Integer] := Module[{approx, year, month, day}, approx = Quotient[date - FrenchEpoch[], 1460969/4000]; year = approx - 1 + MSum[(1)&, approx, (date >= ToFixed[ModifiedFrench[1, 1, #]])&]; month = Quotient[date - ToFixed[ModifiedFrench[1, 1, year]], 30] + 1; day = date - ToFixed[ModifiedFrench[month, 1, year]] + 1; ModifiedFrench[month, day, year] ] (** ;;;; Section: Chinese Calendar **) (** (defun chinese-time-zone (date) ;; TYPE fixed-date -> real ;; The difference (in minutes) of the Beijing time zone ;; from Universal Time on fixed date. (let* ((year (gregorian-year-from-fixed date))) (if (< year 1929) (+ 465 40/60) 480))) **) ChineseTimeZone[date_] := Module[{year}, year = GregorianYearFromFixed[date]; If[year < 1929, 465 + 40/60, 480] ] (** (defun chinese-date-next-solar-longitude (d l) ;; (fixed-date angle) -> fixed-date ;; Fixed date (Beijing time) of the first date on or ;; after fixed date d (Beijing time) when the solar ;; longitude will be a multiple of l degrees. (fixed-from-jd (local-from-universal (date-next-solar-longitude (universal-from-local (jd-from-moment d) (chinese-time-zone d)) l) (chinese-time-zone d)))) **) ChineseDateNextSolarLongitude[d_, l_] := FixedFromJD[ LocalFromUniversal[ DateNextSolarLongitude[ UniversalFromLocal[ JDFromMoment[d], ChineseTimeZone[d]], l ], ChineseTimeZone[d]] ] (** (defun current-major-solar-term (date) ;; TYPE fixed-date -> integer ;; Last Chinese major solar term (zhongqi) before fixed ;; date. (let ((s (solar-longitude (universal-from-local (jd-from-moment date) (chinese-time-zone date))))) (adjusted-mod (+ 2 (quotient s 30)) 12))) **) CurrentMajorSolarTerm[date_] := Module[{s}, s = SolarLongitude[ UniversalFromLocal[ JDFromMoment[date], ChineseTimeZone[date]]]; AdjustedMod[2 + Quotient[s, 30], 12] ] (** (defun major-solar-term-on-or-after (date) ;; TYPE fixed-date -> fixed-date ;; Fixed date (in Beijing) of the first Chinese major ;; solar term (zhongqi) on or after fixed date. The ;; major terms begin when the sun's longitude is a ;; multiple of 30 degrees. (chinese-date-next-solar-longitude date 30)) **) MajorSolarTermOnOrAfter[date_] := ChineseDateNextSolarLongitude[date, 30] (** (defun current-minor-solar-term (date) ;; TYPE fixed-date -> integer ;; Last Chinese minor solar term (jieqi) before date. (let* ((s (solar-longitude (universal-from-local (jd-from-moment date) (chinese-time-zone date))))) (adjusted-mod (+ 3 (quotient (- s 15) 30)) 12))) **) CurrentMinorSolarTerm[date_] := Module[{s}, s = SolarLongitude[ UniversalFromLocal[ JDFromMoment[date], ChineseTimeZone[date]]]; AdjustedMod[3 + Quotient[s - 15, 30], 12] ] (** (defun minor-solar-term-on-or-after (date) ;; TYPE fixed-date -> fixed-date ;; Fixed date (in Beijing) of the first Chinese minor ;; solar term (jieqi) on or after fixed date. The ;; minor terms begin when the sun's longitude is a ;; multiple of 30 degrees. (let* ((d (chinese-date-next-solar-longitude date 15)) (s (solar-longitude (universal-from-local (jd-from-moment d) (chinese-time-zone d))))) (if (= (mod (round s) 30) 0) (chinese-date-next-solar-longitude d 15) d))) **) MinorSolarTermOnOrAfter[date_] := Module[{d, s}, d = ChineseDateNextSolarLongitude[date, 15]; s = SolarLongitude[ UniversalFromLocal[ JDFromMoment[date], ChineseTimeZone[date]]]; If[Mod[Round[s], 30] == 0, ChineseDateNextSolarLongitude[d, 15], d] ] (** (defun chinese-new-moon-before (date) ;; TYPE fixed-date -> fixed-date ;; Fixed date (Beijing) of first new moon before ;; fixed date. (fixed-from-jd (local-from-universal (new-moon-before (universal-from-local (jd-from-moment date) (chinese-time-zone date))) (chinese-time-zone date)))) **) ChineseNewMoonBefore[date_] := FixedFromJD[ LocalFromUniversal[ NewMoonBefore[ UniversalFromLocal[ JDFromMoment[date], ChineseTimeZone[date]]], ChineseTimeZone[date]] ] (** (defun chinese-new-moon-on-or-after (date) ;; TYPE fixed-date -> fixed-date ;; Fixed date (Beijing) of first new moon on or after ;; fixed date. (fixed-from-jd (local-from-universal (new-moon-at-or-after (universal-from-local (jd-from-moment date) (chinese-time-zone date))) (chinese-time-zone date)))) **) ChineseNewMoonOnOrAfter[date_] := FixedFromJD[ LocalFromUniversal[ NewMoonAtOrAfter[ UniversalFromLocal[ JDFromMoment[date], ChineseTimeZone[date]]], ChineseTimeZone[date]] ] (** (defconstant chinese-epoch ;; TYPE fixed-date ;; Fixed date of start of the Chinese calendar. (fixed-from-gregorian (gregorian-date february 15 -2636))) **) ChineseEpoch[] = ToFixed[Gregorian[February[], 15, -2636]] (** (defun no-major-solar-term? (date) ;; TYPE fixed-date -> boolean ;; True if Chinese lunar month starting on date ;; has no major solar term. (= (current-major-solar-term date) (current-major-solar-term (chinese-new-moon-on-or-after (1+ date))))) **) NoMajorSolarTermQ[date_] := CurrentMajorSolarTerm[date] == CurrentMajorSolarTerm[ChineseNewMoonOnOrAfter[date + 1]] (** (defun fixed-from-chinese (c-date) ;; TYPE chinese-date -> fixed-date ;; Fixed date of Chinese date (cycle year month leap ;; day). (let* ((cycle (chinese-cycle c-date)) (year (chinese-year c-date)) (month (chinese-month c-date)) (leap (chinese-leap c-date)) (day (chinese-day c-date)) (g-year; Gregorian year at start of Chinese year (+ ( * (1- cycle) 60); years in prior cycles (1- year) ; prior years this cycle ;; Gregorian year at start of calendar (gregorian-year-from-fixed chinese-epoch))) (new-year (chinese-new-year g-year)) (p; new moon before date--a month too early if ; there was prior leap month that year (chinese-new-moon-on-or-after (+ new-year ( * (1- month) 29)))) (d (chinese-from-fixed p)) (prior-new-moon (if ; If the months match... (and (= month (chinese-month d)) (equal leap (chinese-leap d))) p; ...that's the right month ;; otherwise, there was a prior leap month that ;; year, so we want the next month (chinese-new-moon-on-or-after (1+ p))))) (+ prior-new-moon day -1))) **) ToFixed[date_Chinese] := Module[{month, gYear, newYear, p, d, priorNewMoon}, month = CMonth[date]; gYear = 60 (CCycle[date] - 1) + (CYear[date] - 1) + GregorianYearFromFixed[ChineseEpoch[]]; newYear = ChineseNewYear[gYear]; p = ChineseNewMoonOnOrAfter[newYear + 29 (month - 1)]; d = Chinese[p]; priorNewMoon = If[CMonth[date] == CMonth[d] && CLeap[date] == CLeap[d], p, ChineseNewMoonOnOrAfter[p + 1]]; priorNewMoon + CDay[date] - 1 ] (** (defun chinese-new-year (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Chinese New Year in Gregorian year. (let* ((s1; prior solstice (major-solar-term-on-or-after (fixed-from-gregorian (gregorian-date december 15 (1- g-year))))) (s2; following solstice (major-solar-term-on-or-after (fixed-from-gregorian (gregorian-date december 15 g-year)))) (m1 ; month after last 11th month--either 12 or ; leap 11 (chinese-new-moon-on-or-after (1+ s1))) (m2 ; month after m2--either 1 or leap 12 (chinese-new-moon-on-or-after (1+ m1))) (m11 ; next 11th month (chinese-new-moon-before (1+ s2)))) (if ; Either m1 or m2 is a leap month if there are 13 ; new moons (12 full lunar months) and either m1 or ; m2 has no major solar term (and (= (round (/ (- m11 m1) mean-synodic-month)) 12) (or (no-major-solar-term? m1) (no-major-solar-term? m2))) (chinese-new-moon-on-or-after (1+ m2)) m2))) **) ChineseNewYear[gYear_] := Module[{s1, s2, m1, m2, m11}, s1 = MajorSolarTermOnOrAfter[ ToFixed[Gregorian[December[], 15, gYear - 1]]]; s2 = MajorSolarTermOnOrAfter[ ToFixed[Gregorian[December[], 15, gYear]]]; m1 = ChineseNewMoonOnOrAfter[s1 + 1]; m2 = ChineseNewMoonOnOrAfter[m1 + 1]; m11 = ChineseNewMoonBefore[s2 + 1]; If[Round[(m11 - m1) / MeanSynodicMonth[]] == 12 && (NoMajorSolarTermQ[m1] || NoMajorSolarTermQ[m2]), ChineseNewMoonOnOrAfter[m2 + 1], m2] ] (** (defun chinese-from-fixed (date) ;; TYPE fixed-date -> chinese-date ;; Chinese date (cycle year month leap day) of fixed ;; date. (let* ((g-year (gregorian-year-from-fixed date)) (s1; Prior solstice for most dates (major-solar-term-on-or-after (fixed-from-gregorian (gregorian-date december 15 (1- g-year))))) (s2; Following solstice for most dates--can be the ; prior solstice if the date is at the end of ; the Gregorian year, just after the solstice (major-solar-term-on-or-after (fixed-from-gregorian (gregorian-date december 15 g-year)))) (m1 ; month after last 11th month (if (and (<= s1 date) (< date s2)) (chinese-new-moon-on-or-after (1+ s1)) ;; Date is at end of the Gregorian year, just ;; after the solstice, so we need the solstice ;; of that year (chinese-new-moon-on-or-after (1+ s2)))) (m2 ; next 11th month (if (and (<= s1 date) (< date s2)) (chinese-new-moon-before (1+ s2)) ;; Date is at end of the Gregorian year, just ;; after the solstice, so we need the solstice ;; of the following year (chinese-new-moon-before (1+ (major-solar-term-on-or-after (fixed-from-gregorian (gregorian-date december 15 (1+ g-year)))))))) (m ; start of month containing date (chinese-new-moon-before (1+ date))) (leap-year; if there are 13 new moons (12 full ; lunar months) (= (round (/ (- m2 m1) mean-synodic-month)) 12)) (month ; month number (adjusted-mod (- ;; ordinal position of month in year (round (/ (- m m1) mean-synodic-month)) ;; minus 1 during or after a leap month (if (and leap-year (prior-leap-month? m1 m)) 1 0)) 12)) (leap-month ; it's a leap month if... (and leap-year; ...there are 13 months (no-major-solar-term? m); no major solar term (not (prior-leap-month?; and no prior leap month m1 (chinese-new-moon-before m))))) (elapsed-years ; since the epoch (+ (- (gregorian-year-from-fixed date) (gregorian-year-from-fixed chinese-epoch)) (if (or (< month 11) (> date (fixed-from-gregorian (gregorian-date july 1 g-year)))) 1 0))) (cycle (1+ (quotient (1- elapsed-years) 60))) (year (adjusted-mod elapsed-years 60)) (day (1+ (- date m)))) (chinese-date cycle year month leap-month day))) **) Chinese[date_Integer] := Module[{gYear, s1, s2, m1, m2, m, leapYear, month, leapMonth, elapsedYears, cycle, year, day}, gYear = GregorianYearFromFixed[date]; s1 = MajorSolarTermOnOrAfter[ToFixed[Gregorian[December[], 15, gYear - 1]]]; s2 = MajorSolarTermOnOrAfter[ToFixed[Gregorian[December[], 15, gYear]]]; m1 = If[s1 <= date < s2, ChineseNewMoonOnOrAfter[s1 + 1], ChineseNewMoonOnOrAfter[s2 + 1]]; m2 = If[s1 <= date < s2, ChineseNewMoonBefore[s2 + 1], ChineseNewMoonBefore[ MajorSolarTermOnOrAfter[ ToFixed[Gregorian[December[], 15, gYear]]]] + 1 ]; m = ChineseNewMoonBefore[date + 1]; leapYear = Round[(m2 - m1) / MeanSynodicMonth[]] == 12; month = AdjustedMod[ Round[(m - m1) / MeanSynodicMonth[]] - If[leapYear && PriorLeapMonthQ[m1, m], 1, 0], 12]; leapMonth = leapYear && NoMajorSolarTermQ[m] && !PriorLeapMonthQ[m1, ChineseNewMoonBefore[m]]; elapsedYears = GregorianYearFromFixed[date] - GregorianYearFromFixed[ChineseEpoch[]] + If[month < 11 || date > ToFixed[Gregorian[July[], 1, gYear]], 1, 0]; cycle = Quotient[elapsedYears - 1, 60] + 1; year = AdjustedMod[elapsedYears, 60]; day = date - m + 1; Chinese[cycle, year, month, leapMonth, day] ] (** (defun prior-leap-month? (m-prime m) ;; TYPE (fixed-date fixed-date) -> boolean ;; True if there is a Chinese leap month at or after lunar ;; month m-prime and at or before lunar month m. (and (>= m m-prime) (or (prior-leap-month? m-prime (chinese-new-moon-before m)) (no-major-solar-term? m)))) **) PriorLeapMonthQ[mPrime_, m_] := m >= mPrime && (PriorLeapMonthQ[mPrime, ChineseNewMoonBefore[m]] || NoMajorSolarTermQ[m]) (** (defun chinese-sexagesimal-name (n) ;; TYPE integer -> ({1--10} {1-12}) ;; The n-th name of the Chinese sexagesimal cycle. (list (adjusted-mod n 10) (adjusted-mod n 12))) **) ChineseSexagesimalName[n_] := {AdjustedMod[n, 10], AdjustedMod[n, 12]} (** (defun chinese-name-of-year (y) ;; TYPE chinese-year -> ({1--10} {1-12}) ;; Sexagesimal name for Chinese year y of any cycle. (chinese-sexagesimal-name y)) **) ChineseNameOfYear[y_] := ChineseSexagesimalName[y] (** (defun chinese-name-of-month (y m) ;; TYPE chinese-month -> ({1--10} {1-12}) ;; Sexagesimal name for Chinese month m in year y of any ;; cycle. (chinese-sexagesimal-name (+ ( * 12 y) m 44))) **) ChineseNameOfMonth[y_, m_] := ChineseSexagesimalName[12 y + m + 44] (** (defun chinese-name-of-day (date) ;; TYPE chinese-date -> ({1--10} {1-12}) ;; Sexagesimal name for Chinese date. (chinese-sexagesimal-name (+ (fixed-from-chinese date) 15))) **) ChineseNameOfDay[date_Chinese] := ChineseSexagesimalName[ToFixed[date] + 15] (** (defun dragon-festival (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of the Dragon Festival occurring in ;; Gregorian year. (let* ((elapsed-years (1+ (- g-year (gregorian-year-from-fixed chinese-epoch)))) (cycle (1+ (quotient (1- elapsed-years) 60))) (year (adjusted-mod elapsed-years 60))) (fixed-from-chinese (chinese-date cycle year 5 false 5)))) **) DragonFestival[gYear_] := Module[{elapsedYears, cycle, year}, elapsedYears = gYear - GregorianYearFromFixed[ChineseEpoch[]] + 1; cycle = Quotient[elapsedYears - 1, 60] + 1; year = AdjustedMod[elapsedYears, 60]; ToFixed[Chinese[cycle, year, 5, False, 5]] ] (** ;;;; Section: Modern Hindu Calendars **) (** (defun hindu-sine-table (entry) ;; TYPE integer -> arcminute ;; This simulates the Hindu sine table. ;; entry is an angle given as a multiplier of 225'. (let* ((exact ( * 3438 (sin-degrees ( * entry 225/60)))) (error ( * 0.215 (signum exact) (signum (- (abs exact) 1716))))) (round (+ exact error)))) **) HinduSineTable[entry_] := Module[{exact, error}, exact = 3438 SinDegrees[entry 225/60]; error = 0.215 Sign[exact] Sign[ Abs[exact] - 1716 ]; N[Round[exact + error]] ] (** (defun hindu-sine (theta) ;; TYPE arcminute -> [-3438:3438] ;; theta is in minutes of arc. ;; Linear interpolation in Hindu table is used. (let* ((entry (/ theta 225)) ; Interpolate in table. (fraction (mod entry 1))) (+ ( * fraction (hindu-sine-table (ceiling entry))) ( * (- 1 fraction) (hindu-sine-table (floor entry)))))) **) HinduSine[theta_] := Module[{entry, fraction}, entry = theta / 225; fraction = Mod[entry, 1]; fraction HinduSineTable[Ceiling[entry]] + (1 - fraction) HinduSineTable[Floor[entry]] ] (** (defun hindu-arcsin (units) ;; TYPE [-3438:3438] -> arcminute ;; Inverse of Hindu sine function. (if (< units 0) (- (hindu-arcsin (- units))) (let* ((pos (sum 1 k 0 (> units (hindu-sine-table k)))) (val ; Lower value in table. (hindu-sine-table (1- pos)))) ( * 225 (+ pos -1 ; Interpolate. (/ (- units val) (- (hindu-sine-table pos) val))))))) **) HinduArcsin[units_] := If[units < 0, -HinduArcsin[-units], Module[{pos, val}, pos = MSum[(1)&, 0, (units > HinduSineTable[#])&]; val = HinduSineTable[pos - 1]; 225 (pos - 1 + (units - val) / (HinduSineTable[pos] - val)) ] ] (** (defun mean-position (ky-time period) ;; TYPE (hindu-moment rational) -> arcminute ;; Position in minutes of arc at ky-time ;; in uniform circular orbit of period days. ( * 21600 (mod (/ ky-time period) 1))) **) MeanPosition[kyTime_, period_] := 21600 Mod[kyTime / period, 1] (** (defconstant hindu-sidereal-year ;; TYPE rational ;; Mean length of Hindu sidereal year. (+ 365 279457/1080000)) **) HinduSiderealYear[] = 365 + 279457/1080000 (** (defconstant hindu-creation ;; TYPE integer ;; Days from creation to onset of Hindu epoch. ( * 1955880000 hindu-sidereal-year)) **) HinduCreation[] = 1955880000 HinduSiderealYear[] (** (defconstant hindu-sidereal-month ;; TYPE rational ;; Mean length of Hindu sidereal month. (+ 27 4644439/14438334)) **) HinduSiderealMonth[] = 27 + 4644439/14438334 (** (defconstant hindu-synodic-month ;; TYPE rational ;; Mean time from new moon to new moon. (+ 29 7087771/13358334)) **) HinduSynodicMonth[] = 29 + 7087771/13358334 (** (defconstant hindu-anomalistic-year ;; TYPE rational ;; Time from aphelion to aphelion. (/ 1577917828000 (- 4320000000 387))) **) HinduAnomalisticYear[] = 1577917828000 / (4320000000 - 387) (** (defconstant hindu-anomalistic-month ;; TYPE rational ;; Time from apogee to apogee, with bija correction. (/ 1577917828 (- 57753336 488199))) **) HinduAnomalisticMonth[] = 1577917828 / (57753336 - 488199) (** (defun true-position (ky-time period size anomalistic change) ;; TYPE (hindu-moment rational positive-integer rational ;; TYPE non-negative-integer) -> arcminute ;; Longitudinal position (in arcminutes) at ky-time since ;; Hindu epoch. period is period of mean motion. size is ;; ratio of radii of epicycle and deferent. anomalistic ;; is the period of retrograde revolution about epicycle. ;; change is maximum decrease in epicycle size. (let* ((long ; Position of epicycle center (mean-position ky-time period)) (days (+ ky-time hindu-creation)) ; since creation (offset ; Sine of anomaly (hindu-sine (mean-position days anomalistic))) (contraction ( * (abs offset) change size 1/3438)) (equation ; Equation of center (hindu-arcsin ( * offset (- size contraction))))) (mod (- long equation) 21600))) **) TruePosition[kyTime_, period_, size_, anomalistic_, change_] := Module[{long, days, offset, contraction, equation}, long = MeanPosition[kyTime, period]; days = kyTime + HinduCreation[]; offset = HinduSine[MeanPosition[days, anomalistic]]; contraction = Abs[offset] change size 1/3438; equation = HinduArcsin[offset (size - contraction)]; Mod[long - equation, 21600] ] (** (defun hindu-solar-longitude (ky-time) ;; TYPE hindu-moment -> arcminute ;; Solar longitude in arcminutes at ky-time since epoch. (true-position ky-time hindu-sidereal-year 14/360 hindu-anomalistic-year 1/42)) **) HinduSolarLongitude[kyTime_] := TruePosition[kyTime, HinduSiderealYear[], 14/360, HinduAnomalisticYear[], 1/42] (** (defun hindu-zodiac (ky-time) ;; TYPE hindu-moment -> hindu-solar-month ;; Zodiacal sign of the sun, as integer in range 1..12, ;; at ky-time. (1+ (quotient (hindu-solar-longitude ky-time) 1800))) **) HinduZodiac[kyTime_] := Quotient[HinduSolarLongitude[kyTime], 1800] + 1 (** (defun hindu-lunar-longitude (ky-time) ;; TYPE hindu-moment -> arcminute ;; Lunar longitude in arcminutes at ky-time. (true-position ky-time hindu-sidereal-month 32/360 hindu-anomalistic-month 1/42)) **) HinduLunarLongitude[kyTime_] := TruePosition[kyTime, HinduSiderealMonth[], 32/360, HinduAnomalisticMonth[], 1/42] (** (defun hindu-lunar-phase (ky-time) ;; TYPE hindu-moment -> arcminute ;; Longitudinal distance between the sun and moon ;; in arcminutes at ky-time since Hindu epoch. (mod (- (hindu-lunar-longitude ky-time) (hindu-solar-longitude ky-time)) 21600)) **) HinduLunarPhase[kyTime_] := Mod[HinduLunarLongitude[kyTime] - HinduSolarLongitude[kyTime], 21600] (** (defun lunar-day (ky-time) ;; TYPE hindu-moment -> hindu-lunar-day ;; Phase of moon (tithi) at ky-time, ;; as an integer in the range 1..30. (1+ (quotient (hindu-lunar-phase ky-time) 720))) **) LunarDay[kyTime_] := Quotient[HinduLunarPhase[kyTime], 720] + 1 (** (defun hindu-new-moon (ky-time) ;; TYPE hindu-moment -> hindu-moment ;; Last new moon preceding ky-time since epoch. (let* ((tomorrow (1+ ky-time)) (estimate ; Can be off by 2/3 day. (- tomorrow (mod tomorrow hindu-synodic-month))) (try (binary-search ; Search for phase start. l (- estimate 2/3) u (+ estimate 2/3) x (< (hindu-lunar-phase x) 10800) (or (< ky-time l) ; Wrong new moon (and (<= u ky-time) (= (hindu-zodiac l) (hindu-zodiac u))))))) (if (> try ky-time); Check whether before/after ky-time. (hindu-new-moon ; It's previous month. (- (floor ky-time) 20)) try))) **) HinduNewMoon[kyTime_] := Module[{tomorrow, estimate, try}, tomorrow = kyTime + 1; estimate = tomorrow - Mod[tomorrow, HinduSynodicMonth[]]; try = BinarySearch[estimate - 2/3, estimate + 2/3, Function[{x}, HinduLunarPhase[x] < 10800], Function[{l, u}, kyTime < l || (u <= kyTime && HinduZodiac[l] == HinduZodiac[u])] ]; If[try > kyTime, HinduNewMoon[Floor[kyTime] - 20], try ] ] (** (defun lunar-day-start (ky-time k accuracy) ;; TYPE (hindu-moment rational real) -> hindu-moment ;; Time lunar-day (tithi) number k begins before or at ;; ky-time in days since onset of Hindu epoch. k can be ;; fractional (for karanas). accuracy is in fraction of ;; (civil) day. (let* ((tomorrow (1+ ky-time)) (part ; Mean time from start of month. ( * (1- k) 1/30 hindu-synodic-month)) (estimate ; Mean occurrence of lunar-day. (- tomorrow (mod (- tomorrow part) hindu-synodic-month))) (try (binary-search ; Search for phase. l (- estimate 2/3) u (+ estimate 2/3) x (let* ((diff-x (- (hindu-lunar-phase x) ( * (1- k) 720)))) (or (< 0 diff-x 10800) ; Past the lunar-day (< diff-x -10800))) ; Over the edge (or (< ky-time l) (< (- u l) accuracy))))) (if (> try ky-time) ; Too late. (lunar-day-start ; Previous month. (- ky-time 20) k accuracy) try))) **) LunarDayStart[kyTime_, k_, accuracy_] := Module[{tomorrow, part, estimate, try}, tomorrow = kyTime + 1; part = (k - 1) 1/30 HinduSynodicMonth[]; estimate = tomorrow - Mod[tomorrow - part, HinduSynodicMonth[]]; try = BinarySearch[estimate - 2/3, estimate + 2/3, Function[{x}, Module[{diffX}, diffX = HinduLunarPhase[x] - (k - 1) 720; 0 < diffX < 10800 || diffX < -10800 ] ], Function[{l, u}, kyTime < l || (u - l) < accuracy ] ]; If[try > kyTime, LunarDayStart[kyTime - 20, k, accuracy], try ] ] (** (defun hindu-solar-from-fixed (date) ;; TYPE fixed-date -> hindu-solar-date ;; Hindu solar date equivalent to fixed date. (let* ((ky-time (hindu-day-count date)) ; Hindu date KY. (rise ; Sunrise on Hindu date. (hindu-sunrise ky-time)) (month (hindu-zodiac rise)) (year (- (hindu-calendar-year rise) hindu-solar-era)) (approx ; 3 days before start of mean month. (- ky-time 3 (quotient (mod (hindu-solar-longitude rise) 1800) 60))) (begin (+ approx ; Search forward for beginning... (sum 1 i approx ; ... of month. (/= (hindu-zodiac (hindu-sunrise i)) month)))) (day (- ky-time begin -1))) (gregorian-date month day year))) **) HinduSolar[date_Integer] := Module[{kyTime, rise, month, year, approx, begin, day}, kyTime = HinduDayCount[date]; rise = HinduSunrise[kyTime]; month = HinduZodiac[rise]; year = HinduCalendarYear[rise] - HinduSolarEra[]; approx = kyTime - 3 - Quotient[Mod[HinduSolarLongitude[rise], 1800], 60]; begin = approx + MSum[(1)&, approx, (HinduZodiac[HinduSunrise[#]] =!= month)&]; day = kyTime - begin + 1; HinduSolar[month, day, year] ] (** (defun hindu-solar-precedes? (s-date1 s-date2) ;; TYPE (hindu-solar-date hindu-solar-date) -> boolean ;; True if Hindu solar s-date1 precedes s-date2. (let* ((month1 (standard-month s-date1)) (month2 (standard-month s-date2)) (day1 (standard-day s-date1)) (day2 (standard-day s-date2)) (year1 (standard-year s-date1)) (year2 (standard-year s-date2))) (or (< year1 year2) (and (= year1 year2) (or (< month1 month2) (and (= month1 month2) (< day1 day2))))))) **) HinduSolarPrecedesQ[date1_HinduSolar, date2_HinduSolar] := Module[{month1, month2, day1, day2, year1, year2}, month1 = CMonth[date1]; month2 = CMonth[date2]; day1 = CDay[date1]; day2 = CDay[date2]; year1 = CYear[date1]; year2 = CYear[date2]; year1 < year2 || year1 == year2 && (month1 < month2 || month1 == month2 && day1 < day2) ] (** (defun hindu-calendar-year (ky-time) ;; TYPE hindu-moment -> hindu-solar-year ;; Determine solar year at given time ky-time since Hindu ;; epoch. (let* ((mean; Mean solar longitude (arcminutes) at ky-time. ( * 21600 (mod (/ ky-time hindu-sidereal-year) 1))) (real; True longitude. (hindu-solar-longitude ky-time)) (year; Mean year. (quotient ky-time hindu-sidereal-year))) (cond ((> real 20000 1000 mean); Really previous year. (1- year)) ((> mean 20000 1000 real); Really next year. (1+ year)) (t year)))) **) HinduCalendarYear[kyTime_] := Module[{mean, real, year}, mean = 21600 Mod[kyTime / HinduSiderealYear[], 1]; real = HinduSolarLongitude[kyTime]; year = Quotient[kyTime, HinduSiderealYear[]]; Which[real > 20000 && 1000 > mean, year - 1, mean > 20000 && 1000 > real, year + 1, True, year] ] (** (defconstant hindu-solar-era ;; TYPE standard-year ;; Years from Kali Yuga until Saka era. 3179) **) HinduSolarEra[] = 3179 (** (defun fixed-from-hindu-solar (s-date) ;; TYPE hindu-solar-date -> fixed-date ;; Fixed date corresponding to Hindu solar date s-date ;; (Saka era). Returns bogus for nonexistent date. (let* ((month (standard-month s-date)) (day (standard-day s-date)) (year (standard-year s-date)) (approx; Approximate date from below ; by adding days... (floor (+ ( * (+ year hindu-solar-era ; in months... (/ (1- month) 12)) hindu-sidereal-year) ; ... and years hindu-epoch ; and days before fixed date 0. day -9))); Potential discrepancy of mean date. (try (+ approx ; Search forward to correct date, (sum 1 i approx; or just past it. (hindu-solar-precedes? (hindu-solar-from-fixed i) s-date))))) (if (equal (hindu-solar-from-fixed try) s-date) try bogus))); Date nonexistent on Hindu solar calendar. **) ToFixed[date_HinduSolar] := Module[{approx, try}, approx = Floor[(CYear[date] + HinduSolarEra[] + (CMonth[date] - 1) / 12) HinduSiderealYear[] + HinduEpoch[] + CDay[date] - 9]; try = approx + MSum[(1)&, approx, (HinduSolarPrecedesQ[HinduSolar[#], date])&]; If[HinduSolar[try] === date, try, Bogus] ] (** (defconstant hindu-lunar-era ;; TYPE standard-year ;; Years from Kali Yuga until Vikrama era 3044) **) HinduLunarEra[] = 3044 (** (defun hindu-lunar-from-fixed (date) ;; TYPE fixed-date -> hindu-lunar-date ;; Hindu lunar date equivalent to fixed date. (let* ((ky-time (hindu-day-count date)) ; Hindu date. (rise (hindu-sunrise ky-time)) ; Sunrise that day. (day (lunar-day rise)); Day of month. (leapday ; If previous day the same. (= day (lunar-day (hindu-sunrise (1- ky-time))))) (last-new-moon (hindu-new-moon rise)) (next-new-moon (hindu-new-moon (+ (floor last-new-moon) 35))) (solar-month ; Solar month name. (hindu-zodiac last-new-moon)) (leapmonth ; If begins and ends in same sign. (= solar-month (hindu-zodiac next-new-moon))) (month ; Month of lunar year. (adjusted-mod (1+ solar-month) 12)) (year ; Solar year at next new moon. (- (hindu-calendar-year next-new-moon) hindu-lunar-era ; Era ;; If month is leap, it belongs to next month's ;; year. (if (and leapmonth (= month 1)) -1 0)))) (hindu-lunar-date month leapmonth day leapday year))) **) HinduLunar[date_Integer] := Module[{kyTime, rise, day, leapDay, lastNewMoon, nextNewMoon, solarMonth, leapMonth, month, year}, kyTime = HinduDayCount[date]; rise = HinduSunrise[kyTime]; day = LunarDay[rise]; leapDay = day == LunarDay[HinduSunrise[kyTime - 1]]; lastNewMoon = HinduNewMoon[rise]; nextNewMoon = HinduNewMoon[Floor[lastNewMoon] + 35]; solarMonth = HinduZodiac[lastNewMoon]; leapMonth = solarMonth == HinduZodiac[nextNewMoon]; month = AdjustedMod[solarMonth + 1, 12]; year = HinduCalendarYear[nextNewMoon] - HinduLunarEra[] - If[leapMonth && month == 1, -1, 0]; HinduLunar[month, leapMonth, day, leapDay, year] ] (** (defun hindu-lunar-precedes? (l-date1 l-date2) ;; TYPE (hindu-lunar-date hindu-lunar-date) -> boolean ;; True if Hindu lunar l-date1 precedes l-date2. (let* ((month1 (hindu-lunar-month l-date1)) (month2 (hindu-lunar-month l-date2)) (leap1 (hindu-lunar-leap-month l-date1)) (leap2 (hindu-lunar-leap-month l-date2)) (day1 (hindu-lunar-day l-date1)) (day2 (hindu-lunar-day l-date2)) (leapday1 (hindu-lunar-leap-day l-date1)) (leapday2 (hindu-lunar-leap-day l-date2)) (year1 (hindu-lunar-year l-date1)) (year2 (hindu-lunar-year l-date2))) (or (< year1 year2) (and (= year1 year2) (or (< month1 month2) (and (= month1 month2) (or (and leap1 (not leap2)) (and (equal leap1 leap2) (or (< day1 day2) (and (= day1 day2) (not leapday1) leapday2)))))))))) **) HinduLunarPrecedesQ[date1_HinduLunar, date2_HinduLunar] := Module[{month1, month2, leap1, leap2, day1, day2, leapday1, leapday2, year1, year2}, month1 = CMonth[date1]; month2 = CMonth[date2]; leap1 = CLeapMonth[date1]; leap2 = CLeapMonth[date2]; day1 = CDay[date1]; day2 = CDay[date2]; leapday1 = CLeapDay[date1]; leapday2 = CLeapDay[date2]; year1 = CYear[date1]; year2 = CYear[date2]; year1 < year2 || year1 == year2 && (month1 < month2 || month1 == month2 && (leap1 && !leap2 || leap1 == leap2 && (day1 < day2 || day1 == day2 && !leapday1 && leapday2))) ] (** (defun fixed-from-hindu-lunar (l-date) ;; TYPE hindu-lunar-date -> fixed-date ;; Fixed date corresponding to Hindu lunar date l-date; ;; returns bogus if no such date exists. (let* ((year (hindu-lunar-year l-date)) (month (hindu-lunar-month l-date)) (leap (hindu-lunar-leap-month l-date)) (leap-day (hindu-lunar-leap-day l-date)) (day (hindu-lunar-day l-date)) (ky-year (+ year hindu-lunar-era)) (mean; Approximate date. (fixed-from-old-hindu-lunar (old-hindu-lunar-date month leap day ky-year))) (approx; Check if one month off either way. (cond ((hindu-lunar-precedes? (hindu-lunar-from-fixed (+ mean 15)) l-date) (+ mean hindu-synodic-month)) ((hindu-lunar-precedes? l-date (hindu-lunar-from-fixed (- mean 15))) (- mean hindu-synodic-month)) (t mean))) (try ; Search for correct date ; or just before or after it. (floor (binary-search l (- approx 4) u (+ approx 4) d (not (hindu-lunar-precedes? (hindu-lunar-from-fixed (floor d)) l-date)) (<= (- u l) 2))))) (cond ((equal (hindu-lunar-from-fixed try) l-date) try) ((equal (hindu-lunar-from-fixed (1+ try)) l-date) (1+ try)) ((equal (hindu-lunar-from-fixed (1- try)) l-date) (1- try)) (t bogus)))); Nonexistent Hindu lunar calendar date. **) ToFixed[date_HinduLunar] := Module[{kyYear, mean, approx, try}, kyYear = CYear[date] + HinduLunarEra[]; mean = ToFixed[OldHinduLunar[CMonth[date], CLeapMonth[date], CDay[date], kyYear]]; approx = Which[ HinduLunarPrecedesQ[HinduLunar[mean + 15], date], mean + HinduSynodicMonth[], HinduLunarPrecedesQ[date, HinduLunar[mean - 15]], mean - HinduSynodicMonth[], True, mean ]; try = Floor[BinarySearch[approx - 4, approx + 4, Function[{d}, !HinduLunarPrecedesQ[HinduLunar[Floor[d]], date] ], Function[{l, u}, u - l <= 2 ]] ]; Which[ HinduLunar[try] === date, try, HinduLunar[try + 1] === date, try + 1, HinduLunar[try - 1] === date, try - 1, True, Bogus ] ] (** (defun hindu-equation-of-time (ky-time) ;; TYPE hindu-moment -> hindu-moment ;; Time from mean to true midnight. ;; (This is a gross approximation to the correct value.) (let* ((offset (hindu-sine (mean-position (+ hindu-creation ky-time) hindu-anomalistic-year))) (equation-sun ; Sun's equation of center ; Arcsin is not needed since small ( * offset (- (/ (abs offset) 3713040) 14/360)))) ( * (daily-motion ky-time) equation-sun hindu-sidereal-year 1/21600 1/21600))) **) HinduEquationOfTime[kyTime_] := Module[{offset, equationSun}, offset = HinduSine[MeanPosition[HinduCreation[] + kyTime, HinduAnomalisticYear[]]]; equationSun = offset (Abs[offset] / 3713040 - 14/360); DailyMotion[kyTime] equationSun HinduSiderealYear[] 1/21600 1/21600 ] (** (defun ascensional-difference (ky-time latitude) ;; TYPE (hindu-moment arcminute) -> arcminute ;; Difference between right and oblique ascension ;; of sun (in arcminutes) at ky-time at latitude ;; (in arcminutes). (let* ((sin-decl ( * 1397/3438 (hindu-sine (tropical-longitude ky-time)))) (diurnal-radius (hindu-sine (- 5400 (hindu-arcsin sin-decl)))) (tan ; Tangent of latitude as rational number. (/ (hindu-sine latitude) (hindu-sine (+ 5400 latitude)))) (earth-sine ( * sin-decl tan))) (hindu-arcsin ( * -3438 (/ earth-sine diurnal-radius))))) **) AscensionalDifference[kyTime_, latitude_] := Module[{sinDecl, diurnalRadius, tan}, sinDecl = 1397/3438 HinduSine[TropicalLongitude[kyTime]]; diurnalRadius = HinduSine[5400 - HinduArcsin[sinDecl]]; tan = HinduSine[latitude] / HinduSine[5400 + latitude]; earthSine = sinDecl tan; HinduArcsin[-3438 earthSine / diurnalRadius] ] (** (defun tropical-longitude (ky-time) ;; TYPE hindu-moment -> arcminute ;; Hindu tropical longitude at ky-time days since epoch. ;; Assumes precession with maximum of 27 degrees (1620 ;; minutes) and period of 7200 sidereal years (= ;; 1577917828/600 days). (let* ((midnight (floor ky-time)) ; Whole days. (precession (- 1620 (abs (- 1620 (mod ( * 6480 600/1577917828 midnight) 6480)))))) (mod (- (hindu-solar-longitude ky-time) precession) 21600))) **) TropicalLongitude[kyTime_] := Module[{midnight, precession}, midnight = Floor[kyTime]; precession = 1620 - Abs[1620 - Mod[6480 600/1577917828 midnight, 6480]]; Mod[HinduSolarLongitude[kyTime] - precession, 21600] ] (** (defun rising-sign (ky-time) ;; TYPE hindu-moment -> integer ;; Tabulated speed of rising of current zodiacal sign. (nth (mod (quotient (tropical-longitude ky-time) 1800) 6) (list 1670 1795 1935 1935 1795 1670))) **) RisingSign[kyTime_] := Part[{1670, 1795, 1935, 1935, 1795, 1670}, Mod[Quotient[TropicalLongitude[kyTime], 1800], 6] + 1] (** (defun daily-motion (ky-time) ;; TYPE hindu-moment -> arcminute ;; Sidereal daily motion (in arcminutes) of sun at ky-time. (let* ((mean-motion ; mean daily motion in arcminutes. (/ 21600 hindu-sidereal-year)) (anomaly (mean-position (+ hindu-creation ky-time) hindu-anomalistic-year)) (epicycle ; Current size of epicycle (- 14/360 (/ (abs (hindu-sine anomaly)) 3713040))) (entry (quotient anomaly 225)) (sine-table-step ; Marginal change in anomaly (- (hindu-sine-table (1+ entry)) (hindu-sine-table entry))) (equation-of-motion-factor ( * sine-table-step -1/225 epicycle))) ( * mean-motion (1+ equation-of-motion-factor)))) **) DailyMotion[kyTime_] := Module[{meanMotion, anomaly, epicycle, entry, sineTableStep, equationOfMotionFactor}, meanMotion = 21600 / HinduSiderealYear[]; anomaly = MeanPosition[HinduCreation[] + kyTime, HinduAnomalisticYear[]]; epicycle = 14/360 - Abs[HinduSine[anomaly]] / 3713040; entry = Quotient[anomaly, 225]; sineTableStep = HinduSineTable[entry + 1] - HinduSineTable[entry]; equationOfMotionFactor = sineTableStep (-1/225) epicycle; meanMotion (equationOfMotionFactor + 1) ] (** (defun solar-sidereal-difference (ky-time) ;; TYPE hindu-moment -> arcminute ;; Difference between solar and sidereal day at ky-time. ( * (daily-motion ky-time) (rising-sign ky-time) 1/1800)) **) SolarSiderealDifference[kyTime_] := DailyMotion[kyTime] RisingSign[kyTime] 1/1800 (** (defun hindu-sunrise (ky-time) ;; TYPE hindu-moment -> hindu-moment ;; Sunrise at Ujjain (latitude 1389 minutes) ;; ky-time (whole) days since Hindu epoch. (+ ky-time 1/4 ; Mean sunrise. (hindu-equation-of-time ky-time) ; Apparent midnight. ( * ; Convert sidereal arcminutes to fraction of civil 1577917828/1582237828 1/21600 ; ... day. (+ (ascensional-difference ky-time 1389) (/ (solar-sidereal-difference ky-time) 4))))) **) HinduSunrise[kyTime_] := kyTime + 1/4 + HinduEquationOfTime[kyTime] + 1577917828/1582237828 1/21600 * (AscensionalDifference[kyTime, 1389] + SolarSiderealDifference[kyTime] / 4) (** (defun sunrise-at-ujjain (ky-time) ;; TYPE hindu-moment -> hindu-moment ;; Astronomical sunrise at Ujjain (latitude 1389 minutes) ;; on ky-time (whole) days since Hindu epoch. (let* ((d (+ ky-time hindu-epoch)) (latitude 1389/60) (longitude 4546/60)) (+ ky-time (sunrise d latitude longitude)))) **) SunriseAtUjjain[kyTime_] := Module[{d, latitude, longitude}, d = kyTime + HinduEpoch[]; latitude = 1389/60; longitude = 4546/60; kyTime + Sunrise[d, latitude, longitude] ] (** (defun lunar-mansion (date) ;; TYPE fixed-date -> positive-integer ;; Hindu lunar mansion (nakshatra) at sunrise on fixed date. (let* ((rise (hindu-sunrise (hindu-day-count date)))) (1+ (quotient (hindu-lunar-longitude rise) 800)))) **) LunarMansion[date_] := Module[{rise}, rise = HinduSunrise[HinduDayCount[date]]; Quotient[HinduLunarLongitude[rise], 800] + 1 ] (** (defun samkranti (g-year m) ;; TYPE (gregorian-year hindu-solar-month) -> moment ;; Fixed moment of start of m-th solar month ;; of Hindu year beginning in Gregorian g-year. (let* ((diff (+ hindu-solar-era (standard-year (gregorian-from-fixed hindu-epoch)))) (h-year (- g-year diff)) ; Since K.Y. (ny (- (fixed-from-hindu-solar ; First of month. (hindu-solar-date m 1 h-year)) hindu-epoch)) (begin (binary-search ; Search for exact time. start (- ny 7/8) end (+ ny 3/8) x (if (= m 1) ; Look for 360-0 break. (< (hindu-solar-longitude x) 1800) (>= (hindu-solar-longitude x) ( * (1- m) 1800))) (>= 1/1000000 (- end start))))) (+ begin hindu-epoch))) **) Samkranti[gYear_, m_] := Module[{diff, hYear, ny, begin}, diff = HinduSolarEra[] + CYear[Gregorian[HinduEpoch[]]]; hYear = gYear - diff; ny = ToFixed[HinduSolar[m, 1, hYear]] - HinduEpoch[]; begin = BinarySearch[ny - 7/8, ny + 3/8, Function[{x}, If[m == 1, HinduSolarLongitude[x] < 1800, HinduSolarLongitude[x] >= (m - 1) 1800 ] ], Function[{start, end}, end - start < 1/1000000 ] ]; begin + HinduEpoch[] ] (** (defun mesha-samkranti (g-year) ;; TYPE gregorian-year -> moment ;; Fixed moment of Mesha samkranti (Vernal equinox) ;; in Gregorian g-year. (samkranti g-year 1)) **) MeshaSamkranti[gYear_] := Samkranti[gYear, 1] (** (defun hindu-lunar-new-year (g-year) ;; TYPE gregorian-year -> fixed-date ;; Fixed date of Hindu luni-solar new year ;; in Gregorian g-year. (let* ((mesha ; Fixed date of vernal equinox. (mesha-samkranti g-year)) (m1 ; Prior new moon. (lunar-day-start (hindu-day-count mesha) 1 1/100000)) (m0 ; New moon before that. (lunar-day-start (- m1 27) 1 1/100000)) (new-moon ; First new moon of year. (if ; leap month... (= (hindu-zodiac m0) (hindu-zodiac m1)) m0 m1)) (h-day (floor new-moon)) (rise ; Sunrise that day. (hindu-sunrise h-day))) (+ hindu-epoch h-day ;; Next day if new moon after sunrise, ;; unless lunar day ends before next sunrise. (if (or (< new-moon rise) (= (lunar-day (hindu-sunrise (1+ h-day))) 2)) 0 1)))) **) HinduLunarNewYear[gYear_] := Module[{mesha, m1, m0, newMoon, hDay, rise}, mesha = MeshaSamkranti[gYear]; m1 = LunarDayStart[HinduDayCount[mesha], 1, 1/100000]; m0 = LunarDayStart[m1 - 27, 1, 1/100000]; newMoon = If[HinduZodiac[m0] == HinduZodiac[m1], m0, m1]; hDay = Floor[newMoon]; rise = HinduSunrise[hDay]; HinduEpoch[] + hDay + If[newMoon < rise || LunarDay[HinduSunrise[hDay + 1]] == 2, 0, 1] ] (** (defun karana (n) ;; TYPE {1-60} -> {0-10} ;; Number (0-10) of the name of the n-th (1-60) Hindu karana. (cond ((= n 1) 0) ((> n 57) (- n 50)) (t (adjusted-mod (1- n) 7)))) **) Karana[n_] := Which[n == 1, 0, n > 57, n - 50, True, AdjustedMod[n - 1, 7] ] (** (defun yoga (ky-time) ;; TYPE hindu-moment -> {1-27} ;; Hindu yoga at ky-time days since Hindu epoch. (1+ (floor (mod (/ (+ (hindu-solar-longitude ky-time) (hindu-lunar-longitude ky-time)) 800) 27)))) **) Yoga[kyTime_] := Floor[ Mod[(HinduSolarLongitude[kyTime] + HinduLunarLongitude[kyTime]) / 800, 27] ] + 1 (** (defun sacred-wednesdays-in-gregorian (g-year) ;; TYPE gregorian-year -> list-of-fixed-dates ;; List of Wednesdays in Gregorian year g-year ;; that are day 8 of Hindu lunar months. (sacred-wednesdays (fixed-from-gregorian ; From beginning of year. (gregorian-date january 1 g-year)) (fixed-from-gregorian ; To end. (gregorian-date december 31 g-year)))) **) SacredWednesdaysInGregorian[gYear_] := SacredWednesdays[ ToFixed[Gregorian[January[], 1, gYear]], ToFixed[Gregorian[December[], 31, gYear]] ] (** (defun sacred-wednesdays (start end) ;; TYPE gregorian-year -> list-of-fixed-dates ;; List of Wednesdays between fixed dates start and end ;; (inclusive) that are day 8 of Hindu lunar months. (if (> start end) nil (let* ((wed (kday-on-or-after start wednesday)) (h-date (hindu-lunar-from-fixed wed))) (append (if (= (hindu-lunar-day h-date) 8) (list wed) nil) (sacred-wednesdays (1+ wed) end))))) **) SacredWednesdays[start_, end_] := Module[{result = {}, wed, hDate}, wed = KDayOnOrAfter[start, Wednesday[]]; While[wed <= end, hDate = HinduLunar[wed]; Print[hDate]; If[CDay[hDate] == 8, result = Append[result, wed]]; wed += 7; ]; result ] (* ADDED *) ConvertDateTo[date_, calendar_] := calendar[ToFixed[date]] (* These test tables (except where specified otherwise) were transcribed from Appendix C of _Calendrical Calculations_ *) TestFixedDates[] = { -214193, -61387, 25469, 49217, 171307, 210155, 253427, 369740, 400085, 434355, 452605, 470160, 473837, 507850, 524156, 544676, 567118, 569477, 601716, 613424, 626596, 645554, 664224, 671401, 694799, 704424, 708842, 709409, 709580, 727274, 728714, 744313, 764652} TestDaysOfWeek[] = { Sunday, Wednesday, Wednesday, Sunday, Wednesday, Monday, Saturday, Sunday, Sunday, Friday, Saturday, Friday, Sunday, Sunday, Wednesday, Saturday, Saturday, Saturday, Wednesday, Sunday, Friday, Sunday, Monday, Wednesday, Sunday, Sunday, Monday, Monday, Thursday, Tuesday, Sunday, Wednesday, Sunday} TestJulianDays[] = { 1507231.5, 1660037.5, 1746893.5, 1770641.5, 1892731.5, 1931579.5, 1974851.5, 2091164.5, 2121509.5, 2155779.5, 2174029.5, 2191584.5, 2195261.5, 2229274.5, 2245580.5, 2266100.5, 2288542.5, 2290901.5, 2323140.5, 2334848.5, 2348020.5, 2366978.5, 2385648.5, 2392825.5, 2416223.5, 2425848.5, 2430266.5, 2430833.5, 2431004.5, 2448698.5, 2450138.5, 2465737.5, 2486076.5} TestDates[Gregorian] = { {7, 24, -586}, {12, 5, -168}, {9, 24, 70}, {10, 2, 135}, {1, 8, 470}, {5, 20, 576}, {11, 10, 694}, {4, 25, 1013}, {5, 24, 1096}, {3, 23, 1190}, {3, 10, 1240}, {4, 2, 1288}, {4, 27, 1298}, {6, 12, 1391}, {2, 3, 1436}, {4, 9, 1492}, {9, 19, 1553}, {3, 5, 1560}, {6, 10, 1648}, {6, 30, 1680}, {7, 24, 1716}, {6, 19, 1768}, {8, 2, 1819}, {3, 27, 1839}, {4, 19, 1903}, {8, 25, 1929}, {9, 29, 1941}, {4, 19, 1943}, {10, 7, 1943}, {3, 17, 1992}, {2, 25, 1996}, {11, 10, 2038}, {7, 18, 2094}} TestDates[ISO] = { {29, 7, -586}, {49, 3, -168}, {39, 3, 70}, {39, 7, 135}, {2, 3, 470}, {21, 1, 576}, {45, 6, 694}, {16, 7, 1013}, {21, 7, 1096}, {12, 5, 1190}, {10, 6, 1240}, {14, 5, 1288}, {17, 7, 1298}, {23, 7, 1391}, {5, 3, 1436}, {14, 6, 1492}, {38, 6, 1553}, {9, 6, 1560}, {24, 3, 1648}, {26, 7, 1680}, {30, 5, 1716}, {24, 7, 1768}, {31, 1, 1819}, {13, 3, 1839}, {16, 7, 1903}, {34, 7, 1929}, {40, 1, 1941}, {16, 1, 1943}, {40, 4, 1943}, {12, 2, 1992}, {8, 7, 1996}, {45, 3, 2038}, {28, 7, 2094}} TestDates[Julian] = { {7, 30, -587}, {12, 8, -169}, {9, 26, 70}, {10, 3, 135}, {1, 7, 470}, {5, 18, 576}, {11, 7, 694}, {4, 19, 1013}, {5, 18, 1096}, {3, 16, 1190}, {3, 3, 1240}, {3, 26, 1288}, {4, 20, 1298}, {6, 4, 1391}, {1, 25, 1436}, {3, 31, 1492}, {9, 9, 1553}, {2, 24, 1560}, {5, 31, 1648}, {6, 20, 1680}, {7, 13, 1716}, {6, 8, 1768}, {7, 21, 1819}, {3, 15, 1839}, {4, 6, 1903}, {8, 12, 1929}, {9, 16, 1941}, {4, 6, 1943}, {9, 24, 1943}, {3, 4, 1992}, {2, 12, 1996}, {10, 28, 2038}, {7, 5, 2094}} (* Easter information from as of 9/9/97 *) TestDates[Easter] = { {4, 12, 1998}, {4, 4, 1999}, {4, 23, 2000}, {4, 15, 2001}, {3, 31, 2002}, {4, 20, 2003}, {4, 11, 2004}, {3, 27, 2005}, {4, 16, 2006}, {4, 8, 2007}, {3, 23, 2008}, {4, 12, 2009}, {4, 4, 2010}} TestDates[Coptic] = { {12, 6, -870}, {4, 12, -451}, {1, 29, -213}, {2, 5, -148}, {5, 12, 186}, {9, 23, 292}, {3, 11, 411}, {8, 24, 729}, {9, 23, 812}, {7, 20, 906}, {7, 7, 956}, {7, 30, 1004}, {8, 25, 1014}, {10, 10, 1107}, {5, 29, 1152}, {8, 5, 1208}, {1, 12, 1270}, {6, 29, 1276}, {10, 6, 1364}, {10, 26, 1396}, {11, 19, 1432}, {10, 14, 1484}, {11, 27, 1535}, {7, 19, 1555}, {8, 11, 1619}, {12, 19, 1645}, {1, 19, 1658}, {8, 11, 1659}, {1, 26, 1660}, {7, 8, 1708}, {6, 17, 1712}, {3, 1, 1755}, {11, 11, 1810}} TestDates[Ethiopic] = { {12, 7, -593}, {4, 13, -174}, {1, 29, 64}, {2, 6, 129}, {5, 13, 463}, {9, 24, 569}, {3, 11, 688}, {8, 25, 1006}, {9, 24, 1089}, {7, 21, 1183}, {7, 8, 1233}, {8, 1, 1281}, {8, 26, 1291}, {10, 10, 1384}, {5, 30, 1429}, {8, 6, 1485}, {1, 13, 1547}, {6, 30, 1553}, {10, 7, 1641}, {10, 27, 1673}, {11, 20, 1709}, {10, 15, 1761}, {11, 27, 1812}, {7, 19, 1832}, {8, 11, 1896}, {12, 20, 1922}, {1, 20, 1935}, {8, 11, 1936}, {1, 27, 1937}, {7, 9, 1985}, {6, 18, 1989}, {3, 1, 2032}, {11, 12, 2087}} TestDates[Islamic] = { {12, 9, -1245}, {2, 23, -813}, {4, 1, -568}, {4, 6, -501}, {10, 17, -157}, {6, 3, -47}, {7, 13, 75}, {10, 5, 403}, {5, 22, 489}, {2, 7, 586}, {8, 7, 637}, {2, 20, 687}, {7, 7, 697}, {7, 1, 793}, {7, 6, 839}, {6, 1, 897}, {9, 30, 960}, {5, 27, 967}, {5, 18, 1058}, {6, 2, 1091}, {8, 4, 1128}, {2, 3, 1182}, {10, 10, 1234}, {1, 11, 1255}, {1, 21, 1321}, {3, 19, 1348}, {9, 8, 1360}, {4, 13, 1362}, {10, 7, 1362}, {9, 13, 1412}, {10, 5, 1416}, {10, 12, 1460}, {3, 5, 1518}} TestDates[Bahai] = { {-6, 6, 3, 7, 12}, {-5, 9, 3, 14, 13}, {-4, 2, 13, 10, 17}, {-4, 6, 2, 11, 6}, {-3, 4, 13, 16, 9}, {-3, 10, 6, 4, 4}, {-3, 16, 10, 13, 7}, {-2, 14, 6, 2, 17}, {-2, 18, 13, 4, 8}, {-1, 4, 12, 1, 3}, {-1, 7, 4, 20, 9}, {-1, 9, 15, 1, 13}, {-1, 10, 6, 2, 19}, {-1, 15, 4, 5, 8}, {-1, 17, 10, 17, 16}, {0, 1, 10, 2, 1}, {0, 4, 14, 10, 12}, {0, 5, 1, 20, 4}, {0, 9, 14, 5, 6}, {0, 11, 8, 6, 7}, {0, 13, 6, 7, 12}, {0, 16, 1, 5, 15}, {0, 18, 14, 8, 2}, {0, 19, 15, 1, 7}, {1, 4, 3, 2, 11}, {1, 5, 10, 9, 6}, {1, 6, 3, 11, 3}, {1, 6, 5, 2, 11}, {1, 6, 5, 11, 11}, {1, 8, 15, 20, 16}, {1, 8, 19, 18, 19}, {1, 11, 5, 13, 7}, {1, 14, 4, 7, 6}} TestDates[Persian] = { {5, 1, -1208}, {9, 14, -790}, {7, 2, -552}, {7, 9, -487}, {10, 18, -153}, {2, 30, -46}, {8, 19, 73}, {2, 5, 392}, {3, 3, 475}, {1, 3, 569}, {12, 20, 618}, {1, 14, 667}, {2, 8, 677}, {3, 22, 770}, {11, 13, 814}, {1, 21, 871}, {6, 28, 932}, {12, 14, 938}, {3, 21, 1027}, {4, 10, 1059}, {5, 2, 1095}, {3, 30, 1147}, {5, 10, 1198}, {1, 7, 1218}, {1, 29, 1282}, {6, 3, 1308}, {7, 7, 1320}, {1, 29, 1322}, {7, 14, 1322}, {12, 27, 1370}, {12, 6, 1374}, {8, 19, 1417}, {4, 28, 1473}} TestDates[Hebrew] = { {5, 10, 3174}, {9, 25, 3593}, {7, 3, 3831}, {7, 9, 3896}, {10, 18, 4230}, {3, 4, 4336}, {8, 13, 4455}, {2, 6, 4773}, {2, 23, 4856}, {1, 7, 4950}, {13, 8, 5000}, {1, 21, 5048}, {2, 7, 5058}, {4, 1, 5151}, {11, 7, 5196}, {1, 3, 5252}, {7, 1, 5314}, {12, 27, 5320}, {3, 20, 5408}, {4, 3, 5440}, {5, 5, 5476}, {4, 4, 5528}, {5, 11, 5579}, {1, 12, 5599}, {1, 22, 5663}, {5, 19, 5689}, {7, 8, 5702}, {1, 14, 5703}, {7, 8, 5704}, {13, 12, 5752}, {12, 5, 5756}, {8, 12, 5799}, {5, 5, 5854}} (* Jewish Holiday information from as of 9/9/97 *) TestDates[YomKippur] = { {10, 11, 1997}, {9, 30, 1998}, {9, 20, 1999}, {10, 9, 2000}, {9, 27, 2001}, {9, 16, 2002}, {10, 6, 2003}, {9, 25, 2004}, {10, 13, 2005}} TestDates[Passover] = { {4, 11, 1998}, {4, 1, 1999}, {4, 20, 2000}, {4, 8, 2001}, {3, 28, 2002}, {4, 17, 2003}, {4, 6, 2004}, {4, 24, 2005}, {4, 13, 2006}} TestDates[Purim] = { {3, 12, 1998}, {3, 2, 1999}, {3, 21, 2000}, {3, 9, 2001}, {2, 26, 2002}, {3, 18, 2003}, {3, 7, 2004}, {3, 25, 2005}, {3, 14, 2006}} TestDates[MayanLongCount] = { {6, 8, 3, 13, 9}, {7, 9, 8, 3, 15}, {8, 1, 9, 8, 11}, {8, 4, 15, 7, 19}, {9, 1, 14, 10, 9}, {9, 7, 2, 8, 17}, {9, 13, 2, 12, 9}, {10, 9, 5, 14, 2}, {10, 13, 10, 1, 7}, {10, 18, 5, 4, 17}, {11, 0, 15, 17, 7}, {11, 3, 4, 13, 2}, {11, 3, 14, 16, 19}, {11, 8, 9, 7, 12}, {11, 10, 14, 12, 18}, {11, 13, 11, 12, 18}, {11, 16, 14, 1, 0}, {11, 17, 0, 10, 19}, {12, 1, 10, 2, 18}, {12, 3, 2, 12, 6}, {12, 4, 19, 4, 18}, {12, 7, 11, 16, 16}, {12, 10, 3, 14, 6}, {12, 11, 3, 13, 3}, {12, 14, 8, 13, 1}, {12, 15, 15, 8, 6}, {12, 16, 7, 13, 4}, {12, 16, 9, 5, 11}, {12, 16, 9, 14, 2}, {12, 18, 18, 16, 16}, {12, 19, 2, 16, 16}, {13, 1, 6, 4, 15}, {13, 4, 2, 13, 14}} TestDates[MayanHaab] = { {12, 11}, {3, 5}, {9, 4}, {12, 5}, {12, 14}, {5, 4}, {7, 14}, {5, 8}, {15, 10}, {15, 8}, {15, 8}, {10, 10}, {17, 11}, {5, 15}, {6, 9}, {6, 13}, {18, 3}, {7, 12}, {6, 18}, {9, 1}, {1, 3}, {19, 1}, {14, 4}, {16, 16}, {14, 18}, {4, 7}, {2, 9}, {4, 19}, {10, 9}, {4, 18}, {4, 17}, {8, 12}, {7, 7}} TestDates[MayanTzolkin] = { {5, 9}, {9, 15}, {12, 11}, {9, 19}, {3, 9}, {7, 17}, {2, 9}, {4, 2}, {7, 7}, {9, 17}, {7, 7}, {12, 2}, {10, 19}, {2, 12}, {6, 18}, {12, 18}, {3, 20}, {9, 19}, {8, 18}, {3, 6}, {6, 18}, {10, 16}, {12, 6}, {13, 3}, {11, 1}, {3, 6}, {1, 4}, {9, 11}, {11, 2}, {12, 16}, {9, 16}, {8, 15}, {2, 14}} TestDates[OldHinduSolar] = { {5, 19, 2515}, {9, 26, 2933}, {7, 11, 3171}, {7, 17, 3236}, {10, 19, 3570}, {2, 28, 3677}, {8, 17, 3795}, {1, 26, 4114}, {2, 24, 4197}, {12, 20, 4290}, {12, 7, 4340}, {12, 30, 4388}, {1, 24, 4399}, {3, 7, 4492}, {10, 28, 4536}, {1, 3, 4593}, {6, 12, 4654}, {11, 27, 4660}, {3, 1, 4749}, {3, 21, 4781}, {4, 13, 4817}, {3, 8, 4869}, {4, 20, 4920}, {12, 13, 4939}, {1, 4, 5004}, {5, 11, 5030}, {6, 15, 5042}, {1, 4, 5044}, {6, 23, 5044}, {12, 2, 5092}, {11, 11, 5096}, {7, 26, 5139}, {4, 2, 5195}} TestDates[OldHinduLunar] = { {6, False, 11, 2515}, {9, False, 26, 2933}, {8, False, 3, 3171}, {8, False, 9, 3236}, {11, True, 19, 3570}, {3, False, 5, 3677}, {9, False, 15, 3795}, {2, False, 7, 4114}, {2, False, 24, 4197}, {1, False, 9, 4291}, {12, False, 9, 4340}, {1, False, 23, 4389}, {2, False, 8, 4399}, {4, False, 2, 4492}, {11, False, 7, 4536}, {1, False, 3, 4593}, {7, False, 2, 4654}, {11, False, 29, 4660}, {3, False, 20, 4749}, {4, False, 4, 4781}, {5, False, 6, 4817}, {4, False, 5, 4869}, {5, False, 12, 4920}, {1, True, 13, 4940}, {1, False, 23, 5004}, {5, False, 21, 5030}, {7, False, 9, 5042}, {1, False, 15, 5044}, {7, False, 9, 5044}, {12, False, 14, 5092}, {12, False, 7, 5096}, {8, False, 14, 5139}, {4, False, 6, 5195}} TestDates[Sunrise] = { 0.21536906, 0.28179130, 0.24240510, 0.24590144, 0.29243450, 0.20816046, 0.26632230, 0.22286705, 0.20679373, 0.25051906, 0.26105749, 0.24078204, 0.22143083, 0.20403811, 0.28638499, 0.23475474, 0.24028035, 0.26528995, 0.20404552, 0.20680618, 0.21580853, 0.20456163, 0.21939959, 0.24696211, 0.22740816, 0.22970884, 0.24457271, 0.22740816, 0.24817497, 0.25495173, 0.27248681, 0.26632230, 0.21280751} TestDates[Sunset] = { 0.79339037, 0.70544950, 0.74634083, 0.73912182, 0.71698393, 0.78724026, 0.71105605, 0.77472226, 0.78903246, 0.75920974, 0.75343803, 0.76433481, 0.77569054, 0.79564710, 0.73318076, 0.76759416, 0.75091598, 0.75092375, 0.79536448, 0.79831158, 0.79296897, 0.79743449, 0.78902734, 0.76108161, 0.77184187, 0.77298475, 0.74180613, 0.77184187, 0.73473767, 0.75684059, 0.74617249, 0.71105605, 0.79556694} TestDates[French] = { {11, 5, -2378}, {3, 14, -1959}, {1, 2, -1721}, {1, 10, -1656}, {4, 19, -1322}, {9, 1, -1216}, {2, 19, -1097}, {8, 5, -779}, {9, 5, -696}, {7, 2, -602}, {6, 20, -552}, {7, 13, -504}, {8, 8, -494}, {9, 23, -401}, {5, 14, -356}, {7, 20, -300}, {13, 2, -239}, {6, 15, -232}, {9, 22, -144}, {10, 12, -112}, {11, 6, -76}, {10, 1, -24}, {11, 14, 27}, {7, 6, 47}, {7, 28, 111}, {12, 7, 137}, {1, 7, 150}, {7, 29, 151}, {1, 15, 152}, {6, 27, 200}, {6, 6, 204}, {2, 20, 247}, {10, 30, 302}} TestDates[ModifiedFrench] = { {11, 4, -2378}, {3, 13, -1959}, {1, 2, -1721}, {1, 10, -1656}, {4, 18, -1322}, {9, 1, -1216}, {2, 19, -1097}, {8, 4, -779}, {9, 5, -696}, {7, 1, -602}, {6, 20, -552}, {7, 13, -504}, {8, 8, -494}, {9, 23, -401}, {5, 13, -356}, {7, 19, -300}, {13, 1, -239}, {6, 14, -232}, {9, 22, -144}, {10, 12, -112}, {11, 6, -76}, {10, 1, -24}, {11, 14, 27}, {7, 6, 47}, {7, 29, 111}, {12, 7, 137}, {1, 7, 150}, {7, 29, 151}, {1, 15, 152}, {6, 27, 200}, {6, 7, 204}, {2, 20, 247}, {11, 1, 302}} TestDates[Chinese] = { {35, 11, 6, False, 12}, {42, 9, 10, False, 27}, {46, 7, 8, False, 4}, {47, 12, 8, False, 9}, {52, 46, 11, False, 20}, {54, 33, 4, False, 5}, {56, 31, 10, False, 15}, {61, 50, 3, False, 7}, {63, 13, 4, False, 24}, {64, 47, 2, False, 9}, {65, 37, 2, False, 9}, {66, 25, 2, False, 23}, {66, 35, 3, False, 9}, {68, 8, 5, False, 2}, {68, 53, 1, False, 8}, {69, 49, 3, False, 4}, {70, 50, 8, False, 2}, {70, 57, 1, False, 29}, {72, 25, 4, True, 20}, {72, 57, 6, False, 5}, {73, 33, 6, False, 6}, {74, 25, 5, False, 5}, {75, 16, 6, False, 12}, {75, 36, 2, False, 13}, {76, 40, 3, False, 22}, {77, 6, 7, False, 21}, {77, 18, 8, False, 9}, {77, 20, 3, False, 15}, {77, 20, 9, False, 9}, {78, 9, 2, False, 14}, {78, 13, 1, False, 7}, {78, 55, 10, False, 14}, {79, 51, 6, False, 7}} TestDates[HinduSolar] = { {5, 18, -664}, {9, 25, -246}, {7, 8, -8}, {7, 15, 57}, {10, 20, 391}, {2, 30, 498}, {8, 15, 616}, {1, 27, 935}, {2, 25, 1018}, {12, 22, 1111}, {12, 9, 1161}, {1, 1, 1210}, {1, 26, 1220}, {3, 7, 1313}, {10, 29, 1357}, {1, 4, 1414}, {6, 9, 1475}, {11, 28, 1481}, {3, 2, 1570}, {3, 21, 1602}, {4, 12, 1638}, {3, 9, 1690}, {4, 19, 1741}, {12, 15, 1760}, {1, 6, 1825}, {5, 9, 1851}, {6, 13, 1863}, {1, 6, 1865}, {6, 20, 1865}, {12, 3, 1913}, {11, 12, 1917}, {7, 23, 1960}, {4, 1, 2016}} TestDates[HinduLunar] = { {6, False, 11, False, -529}, {9, False, 27, False, -111}, {8, False, 3, False, 127}, {8, False, 9, False, 192}, {11, False, 19, False, 526}, {3, False, 5, False, 633}, {9, False, 15, False, 751}, {2, False, 6, False, 1070}, {3, True, 23, False, 1153}, {1, False, 8, False, 1247}, {1, False, 8, False, 1297}, {1, False, 22, False, 1345}, {2, False, 8, False, 1355}, {4, False, 1, False, 1448}, {11, False, 7, False, 1492}, {2, True, 3, False, 1549}, {7, False, 2, False, 1610}, {11, False, 28, True, 1616}, {3, False, 20, False, 1705}, {4, False, 4, False, 1737}, {5, False, 6, False, 1773}, {4, False, 5, False, 1825}, {5, False, 11, False, 1876}, {1, False, 13, False, 1896}, {1, False, 22, False, 1960}, {5, False, 20, False, 1986}, {7, False, 9, False, 1998}, {1, False, 14, False, 2000}, {7, False, 8, False, 2000}, {12, False, 14, False, 2048}, {12, False, 7, False, 2052}, {8, False, 14, False, 2095}, {4, False, 6, False, 2151}} Interleave[list1_, list2_] := Thread[List[list1, list2]] FormatResult[d_] := If[Length[d] == 0, "OK", d] PrintResult[s_, d_] := Module[{result}, WriteString["stdout", s]; result = Evaluate[d]; Print[": ", FormatResult[result]] ] SetAttributes[PrintResult, HoldRest] PrintDotAfter[exp_] := ( WriteString["stdout", "."]; exp ) TestFixedToWeekDay[] := PrintResult["fixed -> week day", Select[Interleave[TestFixedDates[], TestDaysOfWeek[]], (NameFromDayOfWeek[DayOfWeekFromFixed[ #[[1]] ]] =!= #[[2]])&]]; TestFixedToJulianDay[] := PrintResult["fixed -> Julian day", Select[Interleave[TestFixedDates[], TestJulianDays[]], (JDFromFixed[ #[[1]] ] =!= #[[2]])&]]; TestJulianDayToFixed[] := PrintResult["Julian day -> fixed", Select[Interleave[TestJulianDays[], TestFixedDates[]], (FixedFromJD[ #[[1]] ] =!= #[[2]])&]]; TestFixedToCalendar[calendar_] := PrintResult["fixed -> " <> ToString[calendar], Select[Interleave[TestFixedDates[], TestDates[calendar]], (PrintDotAfter[calendar[ #[[1]] ] =!= Apply[calendar, #[[2]] ]])&]]; TestCalendarToFixed[calendar_] := PrintResult[ToString[calendar] <> " -> fixed", Select[Interleave[TestDates[calendar], TestFixedDates[]], (PrintDotAfter[ToFixed[ Apply[calendar, #[[1]] ] ] =!= #[[2]]])&]]; TestCalendar[calendar_] := { TestFixedToCalendar[calendar], TestCalendarToFixed[calendar] } TestHoliday[holiday_] := PrintResult["fixed -> " <> ToString[holiday], Select[TestDates[holiday], (Gregorian[ holiday[ #[[3]] ] ] =!= Apply[Gregorian, #])&]]; TestSun[sun_] := PrintResult["fixed -> " <> ToString[sun], Select[Interleave[TestFixedDates[], TestDates[sun]], (SetAccuracy[sun[ #[[1]], 31.8, 35.2 ], 8] =!= SetAccuracy[ #[[2]], 8 ])&]]; CalendricaTest[] := ( Print["Basic support algorithms:"]; TestFixedToWeekDay[]; TestFixedToJulianDay[]; TestJulianDayToFixed[]; Print["\nArithmetic calendars:"]; Map[TestCalendar, {Gregorian, ISO, Julian, Coptic, Ethiopic, Islamic, Bahai, Persian, Hebrew, MayanLongCount, OldHinduSolar, OldHinduLunar}]; Map[TestFixedToCalendar, {MayanHaab, MayanTzolkin}]; Print["\nSome holidays:"]; Map[TestHoliday, {Easter, YomKippur, Passover, Purim}]; Print["\nAstronomical cycles:"]; Map[TestSun, {Sunrise, Sunset}]; Print["\nAstronomical calendars:"]; Map[TestCalendar, {French, ModifiedFrench, Chinese, HinduSolar, HinduLunar}]; Print["\nDone!"]; ) End[] EndPackage[]