R/Periodicy.R

In vreijs/ARCHAEOCOSMO: In-Line Documentation for ARCHAEOCOSMO

####################################################################
# (c) V. Reijs, 2006-201
# ARCHAEOCOSMO library (leaf from version 1.01)
# astronomical, geodetic and meteorological formula
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program (http://www.fsf.org/licensing/licenses/ );
#    if not, write to the Free Software
#    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
#
# The author (web.victor.reijs@gmail.com ) is interested in constructive feedback.
# http://archaeocosmology.org/eng/archaeocosmology.htm
#
# Explantion of the procdures can also be found at
# http://www.archaeocosmology.org/eng/archaeocosmoprocedures.htm
####################################################################

#Determines which algorimths are used
keuze <- "VR" #for DeltaT can be "VR" or "swisseph"
FormAstroRefrac <-
  "sinclair" #for Astronomical refraction can be "bennetth" or "sinclair"
Chapront <- TRUE

#info on LOD/deltaT formula, gained from Simplex method by V. Reijs, 2006
#using Stephenson 1997 data
#StartYear <- 1857.688 #[year]
#Average <- 1.7223697 #[msec/cy]
#Periodicy <- 1554.422 #[year]
#Amplitude <- 3.712246 #[msec]
#Phase <- -0.03652 #[deg]
#using Stephenson&Morrison, 2004 data
StartYear <- 1820 #[year]
Average <- 1.80546834626888 #[msec/cy]
Periodicy <- 1443.67123144531 #[year]
Amplitude <- 3.75606495492684 #[msec]
Phase <- 0 #[deg]

#epochs
#http://www.ephemeris.com/space-time.html
J1900 <- 2415020
J2000 <- 2451545

#Maya parameter
GMTCorrelation <-
  584283 #GMT Correlation ('Modified THompson 2': http://research.famsi.org/date_mayaLC.php
BaseMayaHour <- 0 #???
ClassicSM <-
  29.53333 #Glyph X of Maya Lunar Series, 1986, J.H. Linden, page 134
CopanSM <-
  29.5302 #Glyph X of Maya Lunar Series, 1986, J.H. Linden, page 134
PalenqueSM <-
  29.53086 #Glyph X of Maya Lunar Series, 1986, J.H. Linden, page 134
ModernMayanSM <- 29.53058573 #average value between 682CE and 878CE
ModernSM <- 29.5305805 #average value between 3114BCE and 878CE
BaseHaab <-
  8 #0.0.0.0.0 is equal to 8 : http://mathdl.maa.org/mathDL/46/?pa=content&sa=viewDocument&nodeId=3536&pf=1
BaseTzolkin <-
  4 #0.0.0.0.0 is equal to 4:Ajaw http: // mathdl.maa.org / mathDL / 46 /  ? pa = content &sa = viewDocument & nodeId = 3536 & pf = 1
BaseGlyphG <-
  0  #0.0.0.0.0 is equal to 0:G9: http://research.famsi.org/date_mayaLC.php
BaseGlyphZ <-
  3  #determined to get value 5 as stated in '2012: Science&Prophecy of ancient Maya', Mark van Stone, page 113
BaseCopanSM <-
  CopanSM - 22 #Glyph X of Maya Lunar Series, 1994, J.H. Linden, page 351
BaseCobaSM <- CopanSM - 23 #Glyph X of Coba Stela 1, Stone, page 44
BasePalenqueSM <- PalenqueSM - 24 #1986, J.H. Linden, page 127
BaseModernSM <-
  ModernSM - 11.4948 #based on new moon before basedate using SkyMap
BaseClassicSM <-
  ClassicSM - 2 #extrapolation of copan and palenque to be determined


###################################################################
S_DateConversion <- function (JDNDays, Yeartype, Epoch) {
  # JDNDays [Day]
  # Yeartype [julian,tropical]
  # Epoch [year]
  # DateCoversion [year]
  
  Yeartype <- tolower(Yeartype)
  if (Epoch == 1900) {
    EpochJDN = J1900
  }
  if (Epoch == 2000) {
    EpochJDN = J2000
  }
  DeltaJDN <- JDNDays - EpochJDN
  if (Yeartype == "julian") {
    Conversion = DeltaJDN * D2Y
  }
  if (Yeartype == "tropical") {
    Conversion = DeltaJDN / TY(DeltaJDN * D2Y / 100 / 2)
  }
  return(Conversion)
}

#' @export
###################################################################
JDutfromDate <- function (DateString,
                          hour = 0,
                          DeltaCorrelation = 0) {
  functionvector <-
    data.frame(DateString, hour, DeltaCorrelation, stringsAsFactors = FALSE)
  
  #    print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_JDutfromDate(
      functionvector$DateString[i],
      functionvector$hour[i],
      functionvector$DeltaCorrelation[i]
    )
  }
  return(ResultVector)
}

#' @export
###################################################################
DatefromJDut <- function (JDNDaysUT,
                          Argument = 0,
                          CalType = 0) {
  functionvector <-
    data.frame(JDNDaysUT,
               Argument,
               CalType)
  
  #    print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_DatefromJDut(
      functionvector$JDNDaysUT[i],
      functionvector$Argument[i],
      functionvector$CalType[i]
    )
  }
  return(ResultVector)
}


###################################################################
S_DatefromJDut <- function(JDNDaysUT,
                           Argument = 0,
                           CalType = 0) {
  # ' JDNDaysUT [-]
  # ' Argument []
  # ' CalType (=0: Greg/Julian 1582 CE, =1: Proleptic Gregorian always, =2: Proleptic Julian always)
  # ' DatefromJDut [yyy/mm/dd hh:mm:ss (xx. cal.)]
  
  # translated from BASIC program of P. Duffett-Smith, Astronomy
  # with your personal computer
  JDNDaysUTi <- JDNDaysUT - J1900
  cal <- "(Jul. cal.)"
  D <- JDNDaysUTi + 0.5
  i <- floor(D)
  fd <- D - i
  if (fd == 1) {
    fd = 0
    i = i + 1
  }
  if (CalType != 2) {
    if ((i > -115860) || (CalType == 1)) {
      A = floor((i / 36524.25) + 0.99835726) + 14
      i = i + 1 + A - floor(A / 4)
      cal = "(Greg. cal.)"
    }
  }
  bb <- floor((i / 365.25) + 0.802601)
  C <- i - floor((365.25 * bb) + 0.750001) + 416
  Gd <- floor(C / 30.6001)
  mn <- Gd - 1
  dy <- C - floor(30.6001 * Gd) + fd
  YR <- bb + 1899
  if (Gd > 13.5) {
    mn <- Gd - 13
  }
  if (mn < 2.5) {
    YR <- bb + 1900
  }
  if (YR < 1) {
    YR <- YR - 1
  }
  di <- floor(dy)
  fd <- fd * 24
  uur <- floor(fd)
  minu <- floor((fd - uur) * 60)
  sec <- floor(((fd - uur) * 60 - minu) * 60)
  if (Argument == -4) {
    Date <- cal
  }
  if (Argument == -3) {
    Date <- fd
  }
  if (Argument == -2) {Date <- fd/24} #TimeValue(Str$(uur) + ":" + Str$(minu) + ":" + Str$(sec))}
  if (Argument == -5) {
    Date <-
      paste(as.character(uur),
            ":",
            as.character(minu),
            ":",
            as.character(sec),
            sep = "")
  }
  if (YR > 0) {
    if (Argument == 0) {
      Date <-
        paste(
          as.character(YR),
          " CE /",
          as.character(mn),
          "/",
          as.character(di),
          " ",
          as.character(uur),
          ":",
          as.character(minu),
          ":",
          as.character(sec),
          " ",
          cal,
          sep = ""
        )
    }
    if (Argument == -1) {
      Date <-
        paste(as.character(YR),
              "/",
              as.character(mn),
              "/",
              as.character(di),
              sep = "")
    }
    if (Argument == 1) {
      Date <- YR
    }
  }
  else {
    if (Argument == 0) {
      Date <-
        paste(
          as.character(-YR),
          " BCE /",
          as.character(mn),
          "/",
          as.character(di),
          " ",
          as.character(uur),
          ":",
          as.character(minu),
          ":",
          as.character(sec),
          " ",
          cal,
          sep = ""
        )
    }
    if (Argument == -1) {
      Date <-
        paste(as.character(YR + 1),
              "/",
              as.character(mn),
              "/",
              as.character(di),
              sep = "")
    }
    if (Argument == 1) {
      Date <- YR + 1
    }
  }
  if (Argument == 2) {
    Date <- mn
  }
  if (Argument == 3) {
    Date <- di
  }
  if (Argument == 4) {
    Date <- uur
  }
  if (Argument == 5) {
    Date <- minu
  }
  if (Argument == 6) {
    Date <- sec
  }
  return(Date)
}


###################################################################
S_JDutfromDate <- function (DateString,
                            hour = 0,
                            DeltaCorrelation = 0) {
  # DateString [yyyy/mm/dd] or [0.0.0.0.0]
  # Hour [-]
  # DeltaCorrelation [Day]
  # JDutfromDate [Days]
  
  Correlation <- GMTCorrelation + DeltaCorrelation
  M <- 1
  D <- 0
  Maya <- 0
  ystring <- DateString
  yend <- unlist(gregexpr(pattern = "/", ystring))[1]
  yendMaya <- unlist(gregexpr(pattern = "\\.", ystring))[1]
  if (yend != -1) {
    Y <- as.numeric(substr(ystring, 1, yend - 1))
    mstring <- substr(ystring, yend + 1, nchar(ystring))
    mend <- unlist(gregexpr(pattern = "/", mstring))[1]
    if (mend != -1) {
      M <- as.numeric(substr(mstring, 1, mend - 1))
      dstring <- substr(mstring, mend + 1, nchar(mstring))
      D <- as.numeric(dstring)
    }
    else {
      M = as.numeric(mstring)
      D = 1
    }
  }
  else {
    if (yendMaya != -1) {
      baktun <- as.numeric(substr(ystring, 1, yendMaya - 1))
      mstring <- substr(ystring, yendMaya + 1, nchar(ystring))
      yendMaya <- unlist(gregexpr(pattern = "\\.", mstring))[1]
      katun <- as.numeric(substr(mstring, 1, yendMaya - 1))
      mstring <- substr(mstring, yendMaya + 1, nchar(mstring))
      yendMaya <- unlist(gregexpr(pattern = "\\.", mstring))[1]
      tun <- as.numeric(substr(mstring, 1, yendMaya - 1))
      mstring <- substr(mstring, yendMaya + 1, nchar(mstring))
      yendMaya <- unlist(gregexpr(pattern = "\\.", mstring))[1]
      winal <- as.numeric(substr(mstring, 1, yendMaya - 1))
      mstring <- substr(mstring, yendMaya + 1, nchar(mstring))
      kin <- as.numeric(mstring)
      Maya <- 1
    }
    else {
      Y <- as.numeric(ystring)
      M <- 1
      D <- 1
    }
  }
  ut <- hour
  if (Maya == 0) {
    stmonth <- sprintf("%02.0f", M)
    stday <- sprintf("%02.0f", D)
    yyear <-
      as.numeric(paste(sprintf("%04.0f", Y), stmonth, stday, sep = ""))
    # P. Bretagnon, Planetary Programs and tables from -4000 to +2800, 1986, page 6
    if (M > 2) {
      Ym <- Y
      mm <- M
    }
    else {
      Ym <- Y - 1
      mm <- M + 12
    }
    dd <- D
    cm <- 0
    if (Ym < 0) {
      cm = -0.75
    }
    dm <- ut / 24
    Bmonth <- 0
    # Gregorian calendar from 1582/10/15, before that Julian calendar
    if (yyear >= 15821015) {
      am <- Intdiffer(Ym / 100)
      Bmonth <- 2 - am + Intdiffer(am / 4)
    }
    JDut <-
      Intdiffer(365.25 * Ym + cm) + Intdiffer(30.6001 * (mm + 1)) + dd + dm + 1720994.5 + Bmonth
  }
  else {
    #GMT correlation
    #Maya date Long Count 0.0.0.0.0 equivalent with 3114 BCE Sept. 6 (Julian Calendar and GMt correlataion) http://en.wikipedia.org/wiki/Mesoamerican_Long_Count_calendar
    JDut <-
      baktun * 144000 + katun * 7200 + tun * 360 + winal * 20 + kin + (ut + BaseMayaHour) / 24 + (Correlation)
  }
  return(JDut)
}

#' @export
###################################################################
Intdiffer <- function (Number) {
  # number
  # Intdiffer (Intdiffer(-4.98)=-4)
  
  if (Number >= 0 || Number == floor(Number)) {
    newinteger = floor(Number)
  }
  else {
    newinteger = floor(Number) + 1
  }
  return(newinteger)
}

#' @export
###################################################################
Sunobliquity <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_Sunobliquity(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_Sunobliquity <- function (JDNDays) {
  # JDNDays [Day]
  # Sunobliquity [deg]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  Tbret <- S_DateConversion(JDNDays, "julian", 2000) / 1000
  if (Chapront) {
    # J Hilton et al, Celest.Mech.Dyn.Astron. 94 (2006), 351
    Tbret <- 10 * Tbret
    obl <-
      (84381.406 + (-46.836769 + (
        -0.0001831 + (0.0020034 + (-0.000000576 + (-0.0000000434) * Tbret) * Tbret) * Tbret
      ) * Tbret) * Tbret) / 3600
  }
  else {
    # P. Bretagnon, Planetary Programs and tables from -4000 to +2800, 1986, page 6
    obl <-
      23.4392911 - 0.130025833 * Tbret - 0.00000430556 * Tbret ^ 2 + 0.000555347 * Tbret ^ 3 - 0.00000142722 * Tbret ^ 4 - 0.000000693528 * Tbret ^ 5 - 0.0000000108472 * Tbret ^ 6 + 0.000000000197778 * Tbret ^ 7
  }
  return(obl)
}

#' @export
###################################################################
SolarDayOpt <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_SolarDayOpt(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_SolarDayOpt <- function (JDNDays) {
  # JDNDays [Day]
  # SolarDayOpt [Hour]
  
  # http://archaeocosmology.org/eng/moonfluct.htm#analysis
  OffSetYear <- (S_JDutfromDate(StartYear, 0) - JDNDays) / 365.25
  DayL <-
    -(OffSetYear / 100 * Average) + Amplitude * sin((Phase + 360 * OffSetYear / Periodicy) * Deg2Rad)
  DayL <- 24 + DayL / 1000 / 3600
  return(DayL)
}

#' @export
###################################################################
SiderealDay <- function (JDNDays, COD = 0) {
  functionvector <- data.frame(JDNDays, COD)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_SiderealDay(functionvector$JDNDays[i], functionvector$COD[i])
  }
  return(ResultVector)
}

###################################################################
S_SiderealDay <- function (JDNDays, COD = 0) {
  # JDNDays [Day]
  # COD [msec/century]
  # SiderealDay [Hour]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  DayL <-
    S_HarmonicSum(S_TropicalYear(JDNDays, "mean") * D2H,
                  S_SolarDay(JDNDays, COD),
                  -1)
  return(DayL)
}

#' @export
###################################################################
HarmonicSum <- function (PeriodA, PeriodB, PlusMin = 1) {
  functionvector <- data.frame(PeriodA, PeriodB, PlusMin)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_HarmonicSum(functionvector$PeriodA[i],
                                    functionvector$PeriodB[i],
                                    functionvector$PlusMin[i])
  }
  return(ResultVector)
}

###################################################################
S_HarmonicSum <- function (PeriodA, PeriodB, PlusMin = 1) {
  # PeriodA [x]
  # PeriodB [x]
  # PlusMin [-1: minus/same direction, 1: plus/opposite direction]
  # HarmonicSum [x]
  
  PlusMin <- -PlusMin
  HSum <- PeriodA * PeriodB / (PlusMin * PeriodA + PeriodB)
  if (HSum <= 0) {
    HSum = PeriodA * PeriodB / (PeriodA + PlusMin * PeriodB)
  }
  return(HSum)
}

##################################################################
S_HD <- function(PeriodA, PeriodB) {

  HSum <- S_HS(PeriodA, -PeriodB)
  return(HSum)
}

###################################################################
S_HS <- function(PeriodA, PeriodB) {
  HSum <- PeriodA * PeriodB / (PeriodB + PeriodA)
  return(HSum)
}


#' @export
###################################################################
TropicalYear <- function (JDNDays, TropType = "mean") {
  functionvector <- data.frame(JDNDays, TropType)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_TropicalYear(functionvector$JDNDays[i], functionvector$TropType[i])
  }
  return(ResultVector)
}

###################################################################
S_TropicalYear <- function (JDNDays, TropType = "mean") {
  # JDNDays [Date]
  # TropType [mean,vsop,spring,autumn,summer,winter]
  # TropicalYear [Day]
  
  TropType = tolower(TropType)
  if (TropType == "mean") {
    Epoch <- 2000
    Julcent <- S_DateConversion(JDNDays, "julian", Epoch) / 100
    YearL <- TY(Julcent)
  }
  if (TropType == "vsop") {
    Epoch <- 2000
    Julmil <- S_DateConversion(JDNDays, "julian", Epoch) / 1000
    YearL <- TYVSOP(Julmil)
  }
  if (TropType == "spring") {
    Day1 <- DayfromAngle(-90, JDNDays)
    Day2 <- DayfromAngle(270, JDNDays)
    YearL <- Day2 - Day1
  }
  if (TropType == "summer") {
    Day1 <- DayfromAngle(0, JDNDays)
    Day2 <- DayfromAngle(360, JDNDays)
    YearL <- Day2 - Day1
  }
  if (TropType == "autumn") {
    Day1 <- DayfromAngle(90, JDNDays)
    Day2 <- DayfromAngle(450, JDNDays)
    YearL = Day2 - Day1
  }
  if (TropType == "winter") {
    Day1 <- DayfromAngle(180, JDNDays)
    Day2 <- DayfromAngle(540, JDNDays)
    YearL <- Day2 - Day1
  }
  return (YearL)
}

#' @export
###################################################################
TY <- function (Julcent) {
  # Julcent [Julian Century]
  # TY [Day]
  
  if (Chapront) {
    YearL = 365.2421896698 - 0.00000615359 * Julcent - 0.000000000729 * Julcent ^ 2 + 0.000000000264 * Julcent ^ 3
  }
  else {
    # http://www.treasure-troves.com/astro/TropicalYear.html
    YearL = 365.2421896698 - 0.00000615359 * Julcent - 0.000000000729 * Julcent ^ 2 + 0.000000000264 * Julcent ^ 3
  }
  return(YearL)
}

#' @export
###################################################################
TYVSOP <- function (Julmil) {
  # Julmil [Julian Millenium]
  # TYVSOP [Day]
  
  #The history of the tropical year, Meeus J. and Savoie, D.
  #British Astronomical association, 102, 1, 1992
  YearL <-
    365.242189623 - 0.000061522 * Julmil - 0.0000000609 * Julmil ^ 2 + 0.00000026525 * Julmil ^ 3
  return(YearL)
}

#' @export
###################################################################
SolarDay <- function (JDNDays, COD = 0) {
  functionvector <- data.frame(JDNDays, COD)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_SolarDay(functionvector$JDNDays[i], functionvector$COD[i])
  }
  return(ResultVector)
}

###################################################################
S_SolarDay <- function (JDNDays, COD = 0) {
  # JDNDays [Day]
  # COD [msec/century]
  # SolarDay [Hour]
  
  if (COD == 0) {
    # determine more complex formula (using Simplex optimization)
    DayL <- S_SolarDayOpt(JDNDays)
  }
  else {
    # Simple lineair formula
    CODi = COD / 1000 * S2H
    # http://archaeocosmology.org/eng/moonfluct.htm
    Julcent <-
      (S_JDttfromJDut(S_JDutfromDate(StartYear, 0), COD) - JDNDays) / 365.25 / 100
    DayL <- 24 - Julcent * CODi
  }
  return(DayL)
}

###################################################################
S_JDttfromJDut <- function (JDNDaysUT, COD = 0) {
  # JDNDaysUT [yyyy/mm/dd]
  # COD [msec/cy]
  # JDttfromJDut [Days]
  
  JDtt <- JDNDaysUT + S_DeltaT(JDNDaysUT, COD) / D2S
  return(JDtt)
}

#' @export
###################################################################
DeltaT <- function (JDNDays, COD = 0) {
  functionvector <- data.frame(JDNDays, COD)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_DeltaTVR(functionvector$JDNDays[i], functionvector$COD[i])
  }
  return(ResultVector)
}

###################################################################
S_DeltaT <- function (JDNDays, COD = 0) {
  # JDNDays [Days]
  # COD [msec/cy]
  # DeltaT [Sec]
  
  if (keuze == "swisseph") {
    DT <- S_DeltaTSE(JDNDays, COD)
  }
  if (keuze == "VR") {
    DT <- S_DeltaTVR(JDNDays, COD)
  }
  return(DT)
}

###################################################################
S_DeltaTVR <- function (JDNDays, COD = 0) {
  # JDNDays [Day]
  # COD [msec/cy]
  # DeltaTVR [Sec]
  
  # Determined by V. Reijs
  OffSetYear <- (S_JDutfromDate(StartYear, 0) - JDNDays) / 365.25
  if (COD == 0) {
    DT <-
      (OffSetYear ^ 2 / 100 / 2 * Average + Periodicy / 2 / pi * Amplitude * (cos((
        2 * pi * OffSetYear / Periodicy
      )) - 1)) * Y2D
  }
  else {
    DT <- OffSetYear ^ 2 / 100 / 2 * COD * Y2D
  }
  DT = DT / 1000
  return(DT)
}

#' @export
###################################################################
EarthRotation <- function (JDNDays, COD = 0) {
  functionvector <- data.frame(JDNDays, COD)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_EarthRotation(functionvector$JDNDays[i], functionvector$COD[i])
  }
  return(ResultVector)
}

###################################################################
S_EarthRotation <- function(JDNDays, COD = 0) {
  # JDNDays [Day]
  # COD [msec/century]
  # EarthRotation [Hour]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  Period = S_HarmonicSum(S_LunarSolarPrecession(JDNDays) * Y2D * D2H,
                         S_SiderealDay(JDNDays, COD),
                         1)
  return(Period)
}

#' @export
###################################################################
LunarSolarPrecession <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  # print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_LunarSolarPrecession(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_LunarSolarPrecession <- function(JDNDays) {
  # JDNDays [Day]
  # LunarSolarPrecession [Year]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  if (Chapront) {
    Tbret <- S_DateConversion(JDNDays, "julian", 2000) / 100
    #PO3rev01
    #    crev1 = (5038.481507 + Tbret * ((-1.0790069 * 2 + Tbret * (-0.00114045 * 3 + Tbret * (0.000132851 * 4 - Tbret * 0.0000000951 * 5))))) / 3600 / 100
    #    Lrev1 = 360 / crev1
    #PO3rev02
    #    crev2 = (5038.48209 + Tbret * ((-1.0789921 * 2 + Tbret * (-0.0011404 * 3 + Tbret * (0.000132851 * 4 - Tbret * 0.0000000951 * 5))))) / 3600 / 100
    #    Lrev2 = 360 / crev2
    ctp <-
      (5028.7946 + Tbret * (2.224066 + Tbret * (0.0002319 - Tbret * 0.00009412))) / 3600 / 100
    Ltp <- 360 / ctp
    #P03
    cP03pa <-
      (5028.796195 + Tbret * ((
        1.1054348 * 2 + Tbret * (
          0.00007964 * 3 + Tbret * (-0.000023857 * 4 - Tbret * 0.0000000383 * 5)
        )
      ))) / 3600 / 100
    LP03pa <- 360 / cP03pa
    Period <- LP03pa
  }
  #        dif = Ltp - LP03pa
  else {
    #derived from P. Bretagnon, Planetary Programs and tables from -4000 to +2800, 1986, page 6
    Tbret <- S_DateConversion(JDNDays, "julian", 2000) / 1000
    Period <-
      360 / (
        13.96971278 + 0.030888083 * 2 * Tbret + 0.0000214778 * 3 * Tbret ^ 2 - 0.0000653656 * 4 * Tbret ^ 3 - 0.000000501528 * 5 * Tbret ^ 4 + 0.000000048475 * 6 * Tbret ^ 5 + 0.0000000036375 * 7 * Tbret ^ 6
      ) * 1000
  }
  return(Period)
}

#' @export
###################################################################
SiderealYear <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_SiderealYear(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_SiderealYear <- function(JDNDays) {
  # JDNDays [Day]
  # SiderealYear [Day]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  if (Chapront) {
    Tbret = S_DateConversion(JDNDays, "julian", 2000) / 100
    YearL = 365.256362953 + Tbret * (0.0000001139 + Tbret * (-0.000000000076 - 0.00000000000169 * Tbret))
  }
  else {
    YearL = S_HarmonicSum(S_TropicalYear(JDNDays, "mean"),
                          S_LunarSolarPrecession(JDNDays) * Y2D,
                          1)
  }
  return(YearL)
}

#' @export
###################################################################
EclipticYear <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_EclipticYear(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_EclipticYear <- function(JDNDays) {
  # JDNDays [Day]
  # EclipticYear [Day]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  # using my own methodology of simple periods
  YearL <-
    S_HarmonicSum(S_TropicalYear(JDNDays, "mean"),
                  S_LunarNodalCycle(JDNDays) * Y2D,-1)
  return(YearL)
}

#' @export
###################################################################
LunarNodalCycle <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_LunarNodalCycle(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_LunarNodalCycle <- function(JDNDays) {
  # JDNDays [Day]
  # LunarNodalCycle [Year]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  Period <-
    S_HD(S_ICRSLunarNodalCycle(JDNDays),
         S_LunarSolarPrecession(JDNDays))
  return(Period)
}

#' @export
###################################################################
ICRSLunarNodalCycle <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_ICRSLunarNodalCycle(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_ICRSLunarNodalCycle <- function(JDNDays) {
  # JDNDays [Day]
  # ICRSLunarNodalCycle [Year]
  
  # (derived by T. Peters from Chapront [2002], page 704)
  # http://archaeocosmology.org/eng/moonfluct.htm
  tret <- S_DateConversion(JDNDays, "julian", 2000) / 100
  Period = (6793.476501 + tret * (0.0124002 + tret * (0.000022325 - tret * 0.00000013985))) * D2Y
  return(Period)
}

#' @export
###################################################################
EarthRotationperSiderealYear <- function (JDNDays, COD = 0) {
  functionvector <- data.frame(JDNDays, COD)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_EarthRotationperSiderealYear(functionvector$JDNDays[i], functionvector$COD[i])
  }
  return(ResultVector)
}

###################################################################
S_EarthRotationperSiderealYear <- function (JDNDays, COD = 0) {
  # JDNDays [Day]
  # COD [msec/century]
  # EarthRotationperSiderealYear [-]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  Period = S_SiderealYear(JDNDays) * D2H / S_EarthRotation(JDNDays, COD)
  return(Period)
}

#' @export
###################################################################
EquatorPrecession <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_EquatorPrecession(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_EquatorPrecession <- function(JDNDays) {
  # JDNDays [Day]
  # EquatorPrecession [Year]
  
  #is same as LunarSolarPrecession
  Period <- S_LunarSolarPrecession(JDNDays)
  return (Period)
}

#' @export
###################################################################
AnomalisticYear <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_AnomalisticYear(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_AnomalisticYear <- function(JDNDays) {
  # JDNDays [Day]
  # AnomalisticYear [Day]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  # using my own methodology of simple periods
  YearL <-
    S_HarmonicSum(S_TropicalYear(JDNDays, "mean"),
                  S_ClimaticPrecession(JDNDays) * Y2D,
                  1)
  return(YearL)
}

#' @export
###################################################################
ClimaticPrecession <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_ClimaticPrecession(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_ClimaticPrecession <- function(JDNDays) {
  # JDNDays [Day]
  # ClimaticPrecession [Year]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  dd <- S_DateConversion(JDNDays, "julian", 1900) / 10000 * Y2D
  # drived from Expl. Suppl, page 98
  Period <-
    (360 / (
      0.0000470684 + 0.0000339 * 2 / 10000 * dd + 0.00000007 * 3 / 10000 * dd ^ 2
    ) * D2Y)
  return(Period)
}

#' @export
###################################################################
AnomalisticMonth <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_AnomalisticMonth(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_AnomalisticMonth <- function(JDNDays) {
  # JDNDays [Day]
  # AnomalisticMonth [Day]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  # using my own methodology of simple periods
  Period <-
    HarmonicSum(SiderealMonth(JDNDays), LunarApseCycle(JDNDays) * Y2D, 1)
  return(Period)
}

#' @export
TropicalMonth <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_TropicalMonth(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_TropicalMonth <- function(JDNDays) {
  # JDNDays [Day]
  # TropicalMonth [Day]
  
  if (Chapront) {
    MonthL <-
      HarmonicSum(LunarSolarPrecession(JDNDays) * Y2D,-SiderealMonth(JDNDays),
                  1)
  }
  else {
    # http://archaeocosmology.org/eng/moonfluct.htm
    dd <- S_DateConversion(JDNDays, "julian", 1900) / 10000 * Y2D
    # drived from Expl. Suppl, page 107
    MonthL <-
      360 / (13.1763965268 - 0.000085 / 10000 * 2 * dd + 0.000000039 / 10000 * 3 * dd ^ 2)
  }
  return(MonthL)
}

#' @export
SynodicMonth <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_SynodicMonth(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_SynodicMonth <- function(JDNDays) {
  # JDNDays [Day]
  # SynodicMonth [Day]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  # using my own methodology of simple periods
  MonthL <-
    S_HarmonicSum(S_TropicalYear(JDNDays, "mean"),
                  S_TropicalMonth(JDNDays),
                  1)
  return(MonthL)
}

#' @export
SiderealMonth <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_SiderealMonth(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_SiderealMonth <- function(JDNDays) {
  # JDNDays [Day]
  # SiderealMonth [Day]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  # using my own methodology of simple periods
  if (Chapront) {
    # derived by T. Peters from Chapront [2002], page 704
    tret <- S_DateConversion(JDNDays, "julian", 2000) / 100
    MonthL <-
      27.32166155356 + tret * (0.000000216673 + tret * (-0.00000000031243 + tret * 1.9989E-12))
  }
  else {
    # using my own methodology of simple periods
    MonthL <-
      S_HarmonicSum(S_LunarSolarPrecession(JDNDays) * Y2D,
                    S_TropicalMonth(JDNDays),
                    1)
  }
  return(MonthL)
}

#' @export
DraconicMonth <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_DraconicMonth(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_DraconicMonth <- function(JDNDays) {
  # JDNDays [Day]
  # DraconicMonth [Day]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  # using my own methodology of simple periods
  MonthL = S_HarmonicSum(S_LunarNodalCycle(JDNDays) * Y2D,
                         S_TropicalMonth(JDNDays),-1)
  return(MonthL)
}

#' @export
Eccentricity <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_Eccentricity(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_Eccentricity <- function(JDNDays) {
  # JDNDays [Day]
  # Eccentricity [-]
  
  # http://archaeocosmology.org/eng/moonfluct.htm
  # from Expl. Suppl, page 98
  Julcent <- S_DateConversion(JDNDays, "julian", 1900) / 100
  Period <-
    0.01675104 - 0.0000418 * Julcent - 0.000000126 * Julcent ^ 2
  return(Period)
}

#' @export
PerihelionNumber <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_PerihelionNumber(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_PerihelionNumber <- function(JDNDays) {
  # JDNDays [Day]
  # PerihelionNumber [Day]
  
  # Determined by V. Reijs (using SkyMap)
  # hthttp://archaeocosmology.org/eng/moonfluct.htm
  # http://archaeocosmology.org/eng/season.htm
  Julcent <- S_DateConversion(JDNDays, "julian", 2000)
  Number <- 197.26 + Julcent * 0.017808333
  return(Number)
}

#' @export
ICRSLunarApseCycle <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_ICRSLunarApseCycle(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_ICRSLunarApseCycle <- function(JDNDays) {
  # JDNDays [Day]
  # ICRSLunarApseCycle [Year]
  
  # (Chapront [2002], page 704)
  # hhttp://archaeocosmology.org/eng/moonfluct.htm
  tret <- S_DateConversion(JDNDays, "julian", 2000) / 100
  Period <-
    (3232.60542496 + tret * (0.0168939 + tret * (0.000029833 - tret * 0.00000018809))) * D2Y
  return(Period)
}

#' @export
LunarApseCycle <- function (JDNDays) {
  functionvector <- data.frame(JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_LunarApseCycle(functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_LunarApseCycle <- function(JDNDays) {
  # JDNDays [Day]
  # LunarApseCycle [Year]
  
  # Determined by V. Reijs
  # http://archaeocosmology.org/eng/moonfluct.htm
  Period <-
    S_HS(S_LunarSolarPrecession(JDNDays),
         S_ICRSLunarApseCycle(JDNDays))
  return(Period)
}

#' @export
###################################################################
DayfromAngle <- function (PathAngle, JDNDays) {
  functionvector <- data.frame(PathAngle, JDNDays)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_DayfromAngle(functionvector$PathAngle[i], functionvector$JDNDays[i])
  }
  return(ResultVector)
}

###################################################################
S_DayfromAngle <- function(PathAngle, JDNDays) {
  # ' PathAngle [Deg]
  # ' JDNDays [Day]
  # ' DayfromAngle [-]
  
  TropYear <- S_TropicalYear(JDNDays, "mean")
  AnoYear <- S_AnomalisticYear(JDNDays)
  Ecc <- S_Eccentricity(JDNDays)
  Perihelion <- S_PerihelionNumber(JDNDays)
  maxerror <- 0.0000001
  richting <- 1
  stap <- 1
  # V. Reijs, http://archaeocosmology.org//eng/season.htm
  dayoud <- (PathAngle / 90) * 365.2424 / 4
  angleoud <-
    S_AngleinSunsPathfromDay(dayoud, TropYear, AnoYear, Ecc, Perihelion)
  difoud <- angleoud - PathAngle
  if (difoud < 0) {
    signoud <- -1
  } else {
    signoud <- 1
  }
  while (abs(difoud) > maxerror) {
    daynew <- dayoud + richting * stap
    anglenew <-
      S_AngleinSunsPathfromDay(daynew, TropYear, AnoYear, Ecc, Perihelion)
    difnew <- anglenew - PathAngle
    if (difnew < 0) {
      signnew <- -1
    } else {
      signnew <- 1
    }
    if (signnew == signoud) {
      if (abs(difnew) > abs(difoud)) {
        richting <- -richting
      }
    }
    else {
      richting <- -richting
      stap <- stap / 2
    }
    difoud <- difnew
    signoud <- signnew
    dayoud <- daynew
  }
  return(dayoud)
}

#' @export
###################################################################
AngleinSunsPathfromDay <-
  function (DaysSummer,
            TropYear,
            AnoYear,
            Ecc,
            Perihelion) {
    functionvector <-
      data.frame(DaysSummer, TropYear, AnoYear, Ecc, Perihelion)
    #  print(functionvector)
    ResultVector <- c(0)
    for (i in 1:nrow(functionvector))
    {
      ResultVector[i] = S_AngleinSunsPathfromDay(
        functionvector$DaysSummer[i],
        functionvector$TropYear[i],
        functionvector$AnoYear[i],
        functionvector$Ecc[i],
        functionvector$Perihelion[i]
      )
    }
    return(ResultVector)
  }

###################################################################
S_AngleinSunsPathfromDay <-
  function(DaysSummer,
           TropYear,
           AnoYear,
           Ecc,
           Perihelion) {
    # ' DaysSummer [-]
    # ' TropYear [Day]
    # ' AnoYear [Day]
    # ' Ecc [-]
    # ' Perihelion [Day]
    # ' AngleinSunsPathfromDay [Deg]
    
    # V. Reijs, 2004, hhttp://archaeocosmology.org/eng/season.htm
    Angle <-
      360 / TropYear * DaysSummer + Ecc * 2 * Rad2Deg * sin(360 / AnoYear * (DaysSummer - Perihelion) * Deg2Rad)
    return(Angle)
  }

#' @export
###################################################################
MayaDatefromJDut <- function (JDNDays,
                              MayaDateType = 0,
                              DeltaCorrelation = 0) {
  functionvector <- data.frame(JDNDays, MayaDateType,
                               DeltaCorrelation)
  #  print(functionvector)
  ResultVector <- c(0)
  for (i in 1:nrow(functionvector))
  {
    ResultVector[i] = S_MayaDatefromJDut(
      functionvector$JDNDays[i],
      functionvector$MayaDateType[i],
      functionvector$DeltaCorrelation[i]
    )
  }
  return(ResultVector)
}

###################################################################
S_MayaDatefromJDut <-
  function(JDNDaysUT,
           MayaDateType = 0,
           DeltaCorrelation = 0) {
    # ' JDNDaysUT [-
    # ' MayaDateType [-1=all, 0=Long count, 1=Tzolkin,2=Haab, 3=Night, 4=Earth ]
    # ' DeltaCorrelation [Day]
    # ' DatefromJDut [baktun.katun.tun.winal.kin Tzolkin Haab]
    
    Correlation <- GMTCorrelation + DeltaCorrelation
    DaysfromStart <- JDNDaysUT - Correlation
    baktun <- floor(DaysfromStart / 144000)
    DaysfromStart <- DaysfromStart - baktun * 144000
    katun <- floor(DaysfromStart / 7200)
    DaysfromStart <- DaysfromStart - katun * 7200
    tun <- floor(DaysfromStart / 360)
    DaysfromStart <- DaysfromStart - tun * 360
    winal <- floor(DaysfromStart / 20)
    DaysfromStart <- DaysfromStart - winal * 20
    kin <- floor(DaysfromStart)
    ut <- DaysfromStart - kin
    MayaDays <-
      baktun * 144000 + katun * 7200 + tun * 360 + winal * 20 + kin
    #method for Tzolkin and Haab from http://mathdl.maa.org/mathDL/46/?pa=content&sa=viewDocument&nodeId=3536&pf=1
    RoundCalendarRest <-
      MayaDays - 18980 * floor((MayaDays) / 18980)
    NightGod <-
      MayaDays + BaseGlyphG - 9 * floor((MayaDays + BaseGlyphG) / 9)
    EarthGod <-
      MayaDays + BaseGlyphZ - 7 * floor((MayaDays + BaseGlyphZ) / 7)
    if (NightGod == 0) {
      NightGod <- 9
    }
    if (EarthGod == 0) {
      EarthGod <- 7
    }
    TzolkinDay <-
      RoundCalendarRest - 13 * floor(RoundCalendarRest / 13) + BaseTzolkin
    TzolkinGod = RoundCalendarRest - 20 * floor(RoundCalendarRest / 20)
    if (TzolkinDay > 13) {
      TzolkinDay <- TzolkinDay - 13
    }
    if (TzolkinGod == 0) {
      TzolkinGod <- 20
    }
    Haab <- RoundCalendarRest - 365 * floor(RoundCalendarRest / 365)
    if (Haab <= 11) {
      HaabGod <- 18
      HaabDay <- Haab + BaseHaab
    }
    else {
      if (Haab <= 16) {
        HaabGod <- 19
        HaabDay <- Haab - 12
      }
      else {
        HaabGod <- floor((Haab - 17) / 20) + 1
        HaabDay <- (Haab - 17) - (HaabGod - 1) * 20
      }
    }
    if (MayaDateType == -1) {
      Date <-
        paste(
          baktun,
          ".",
          katun,
          ".",
          tun,
          ".",
          winal,
          ".",
          kin,
          " ",
          TzolkinDay,
          ":",
          TzolkinGod,
          " ",
          HaabDay,
          ".",
          HaabGod,
          " ",
          NightGod,
          sep = ""
        )
    }
    if (MayaDateType == 0) {
      Date = paste(baktun, ".", katun, ".", tun, ".", winal, ".", kin, sep = "")
    }
    if (MayaDateType == 1) {
      Date <- paste(TzolkinDay, ":", TzolkinGod, sep = "")
    }
    if (MayaDateType == 2) {
      Date <- paste(HaabDay, ".", HaabGod, sep = "")
    }
    if (MayaDateType == 3) {
      Date <- NightGod
    }
    if (MayaDateType == 4) {
      Date <- EarthGod
    }
    if (MayaDateType == 5) {
      Date <- HaabDay
    }
    if (MayaDateType == 6) {
      Date <- HaabGod
    }
    if (MayaDateType == 7) {
      Date <- TzolkinDay
    }
    if (MayaDateType == 8) {
      Date <- TzolkinGod
    }
    if (MayaDateType == -2) {
      Date <- ut
    }
    return(Date)
  }