// Indian Astronomy calculations
// Thanks to Toke Lindegaard Knudsen for the information.

// Much of this is based on his Ph.D dissertation, "The Siddhantasundara of
// Jnanaraja:  A Critical Edition of Select Chapters with English Translation
// and Commentary", Brown University, 2008.

// Some references are made to _Calendrical Calculations_ by Reingold and
// Dershowitz.
//
//  See also:
// http://www.new.dli.ernet.in/rawdataupload/upload/insa/INSA_1/20008275_17.pdf
//   WARNING:  The above is fundamentally flawed in that it states that Julian
// day begins at midnight, when it actually begins at noon.  This creates
// errors in *all* calculations in the paper.

// And a paper by Reingold and Dershowitz:
// http://emr.cs.iit.edu/~reingold/hindu-paper.pdf

use sun.frink

// There is some discrepancy in sources concerning the size of a sidereal year.
//
// Calendrical Calculations gives 1,577,917,500 civil days in a kalpa (see 9.3)
// (equivalent Frink definition is
//   saurayear := 1_577_917_500/4_320_000 days 
// )
// while Knudsen's dissertation (table 3.1, p. 102) gives 1,577,917,828 civil
// days.

// The modern value, using the best-known size of a sidereal year gives
// 1,577,907,477 civil days.  (In Frink notation
//  "4320000 siderealyear -> days" )
//
//   It's unclear which to use for high-accuracy
// calculations.

kalpa := 4_320_000_000 siderealyears
mahayuga := 1/1000 kalpa

saurayear := siderealyear

// A mahayuga is divided into 4 unequal yugas.
krtayuga    := 4/10 mahayuga
tretayuga   := 3/10 mahayuga
dvaparayuga := 2/10 mahayuga
kaliyuga    := 1/10 mahayuga

// At present, we are in a Kaliyuga, which began 4,567*432,000 = 1,972,944,000
// sidereal years after the beginning of the Kalpa.
// The time of the beginning of the
// present Kaliyuga is 6 a.m. on February 18, 3102 BCE (Julian) in Lanka (this
// corresponds to local sunrise in Lanka, which is a mythological city said
// to be located at the intersection of the terrestrial equator and the prime
// meridian of the Indian tradition, which goes through the city of Ujjain in
// Madhya Pradesh).

UjjainLongitude = DMS[75, 46, 6] East

// This is the time correction that has to be *added* to Ujjain time to get
// UTC.  (Note that it is a negative number.)  UTC is found by subtracting
// (5 hours + 3 min + 4.4 seconds) from Ujjain local time.
// Alternately, it is the time that has to be *subtracted* from UTC to get
// Ujjain time.  (Note that it's a negative number, so to get Ujjain local
// time, add (5 hours + 3 min + 4.4 seconds) to UTC.
UjjainMeanTimeCorrection = UjjainLongitude * (day/circle)

IndianMoonRotation := kalpa / 57_753_300_000
IndianMoonApogee   := kalpa / 488_105_858
IndianMoonNode     := kalpa / 232_311_168


// Converts a number which may be in the form 7;12,13 into Devanagari digits.
numberToDevanagari[str] :=
{
   str =~ %s/[,;]/\u0964/g
   str =~ %s/([0-9])/char[char[$1]-char["0"]+ 0x0966]/eg
   return str
}


// Converts a Devanagari number to the form popularized by Otto Neugebauer,
// for example, 14;23,12,14
DevanagariToNeugebauer[str] :=
{
   r = new array
   first = true
   for c = chars[str]
   {
      if (c >= 0x0966 and c <= 0x096f)
         c = c - 0x0966 + char["0"]
      if (c == 0x0964)
      {
         if (first)
            c = char[";"]
         else
            c = char[","]
         first = false
      }
      r.push[c]
   }
   return char[r]
}


// Parses a Devanagari number as an integer part in base 10, then an arbitrary
// number of parts in base 60 separated by the DEVANAGARI DANDA character
parseDevanagariFraction[str] :=
{
   parts = split["\u0964", str] // Break on danda character
   len = length[parts]
   s = parseDevanagariInteger[parts@0]
   base = 1
   for i = 1 to len-1
   {
      base = base * 60
      s = s + parseDevanagariInteger[parts@i]/base
   }
   return s
}


// Parses an integer either as Devanagari characters or Arabic numerals.
parseDevanagariInteger[str] :=
{
   zero = char["0"]
   nine = char["9"]

   s = 0
   for c=chars[str]
   {
      if (c >= 0x0966 and c <= 0x096f)
         c = c - 0x0966
      else
         if (c >= zero and c <= nine)
            c = c - zero
         else
            return "Invalid character in $str"
      s = s*10 + c
   }

   return s
}


// Converts a floating-point number to the sexagesimal form popularized by
// Otto Neugebauer, such as 14;12,13
sexagesimal[n, fracDigits=6] :=
{
   n = round[n, 60^-fracDigits]
   int = trunc[n]                 // Integer part
   frac = n - int

   str = "$int"

   d = 0

   while ((d < fracDigits) and (frac != 0))
   {
      if (d < fracDigits-1)
         part = trunc[frac * 60]
      else
         part = round[frac * 60]
      
      frac = frac*60 - part
      if (d==0)
         str = str + ";"
      else
         str = str + ","
      str = str + part
      d=d+1
   }

   return str
}


// Rata Die from Calendrical Calculations.
// This should be corrected for timezone, although the RD calculations are
// more of a toy calculation which assumes that every day is reckoned at noon
// and omits any consideration of timezone or time of day which are
// oversimplifications which makes it less than useful in the real world.
rd[date] := ((date - # 0001-01-03 00:00 UTC #) / day) + 1


// Calculate the ahargana based on offset from Julian day.  This should
// be corrected based on the time of day at which the ahargana began (e.g.
// midnight or 6 AM) and for the longitude in which the beginning of the
// ahargana was said to begin.
ahargana[date] := (JD[date] - JD[# BC 3102-02-18 06:00 Ujjain #]) / day