R/heliocentric.R

In Susarro/arqastwb: Archaeoastronomy algorithms

#' Cosine auxiliary function
#'
#' @param A Parameter A
#' @param B Parameter B
#' @param C Parameter C
#' @param t Parameter t
#' @return A * cos(B + C * t)
#' @examples
#' jd<-str2jd("1972-03-08)
#' t<-(jd - 2451545.0) / 365250
#' cos_component(54127,0,0,str2jd())
#'

cos_component<- function(A,B,C,t){
  return(A * cos(B + C * t))
}


#' Heliocentric longitude of a planet
#'
#' @param jd Julian day
#' @param celestial_body planet name as string: mercury, venus,earth, mars, jupyter, saturn
#' @return Heliocentric longitude in degrees
#' @examples
#' jd<-julian_day(13.19,11,2028)
#' heliocentric_longitude(jd,"venus")
#'

heliocentric_longitude <- function(jd, celestial_body) {
  t<-(jd - 2451545.0) / 365250

  if(celestial_body=="venus"){

    L0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(317614667,1353968,89892,5477,3456,2372,1664,1438,1317,1201,769,761,708,585,500,429,327,326,232,180,155,128,128,106),
                    B=c(0,5.5931332,5.30650,4.4163,2.6996,2.9938,4.2502,4.1575,5.1867,6.1536,0.816,1.950,1.065,3.998,4.123,3.586,5.677,4.591,3.163,4.653,5.570,4.226,0.962,1.537),
                    C=c(0,10213.2855462,20426.57109,7860.4194,11790.6291,3930.2097,1577.3435,9683.5946,26.2983,30639.8566,9437.763,529.691,775.523,191.448,15720.839,19367.189,5507.553,10404.734,9153.904,1109.379,19651.048,20.775,5661.332,801.821),
                    t=x)))
    })

    L1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1021352943053.0,95708,14445,213,174,152,82,70,52,38,30,25),
                    B=c(0,2.46424,0.51625,1.795,2.655,6.106,5.70,2.68,3.60,1.03,1.25,6.11),
                    C=c(0,10213.28555,20426.57109,30639.857,26.298,1577.344,191.45,9437.76,775.52,529.69,5507.55,10404.73),
                    t=x)))
    })
    L2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(54127,3891,1338,24,19,10,7,6),
                    B=c(0,0.3451,2.0201,2.05,3.54,3.97,1.52,1.00),
                    C=c(0,10213.2855,20426.5711,26.30,30639.86,775.52,1577.34,191.45),
                    t=x)))
    })

    L3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(136,78,26),
                    B=c(4.8004,3.67,0),
                    C=c(10213.286,20426.57,0),
                    t=x)))
    })

    L4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(114,3,2),
                    B=c(3.1416,5.21,2.51),
                    C=c(0,20426.57,10213.29),
                    t=x)))
    })

    L5 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1),
                    B=c(3.14),
                    C=c(0),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((L0 + L1 * t + L2 * t^2 + L3 * t^3 + L4 * t^4 + L5 * t^5) / 10^8)))
  }
  if(celestial_body=="mercury")
  {
    L0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(440250710,40989415,5046294,855347,165590,34562,7583,3560,1803,1726,1590,1365,1017,714,644,451,404,352,345,343,339,325,273,264,260,239,235,217,209,183,182,176,173,142,138,125,118,106),
                    B=c(0,1.48302034,4.4778549,1.165203,4.119692,0.77931,3.7135,1.5120,4.1033,0.3583,2.9951,4.5992,0.8803,1.541,5.303,6.050,3.282,5.242,2.792,5.765,5.863,1.337,2.495,3.917,0.987,0.113,0.267,0.660,2.092,2.629,2.434,4.536,2.452,3.360,0.291,3.721,2.781,4.206),
                    C=c(0,26087.90314157,52175.8062831,78263.709425,104351.612566,130439.51571,156527.4188,1109.3786,5661.3320,182615.322,25028.5212,27197.2817,31749.2352,24978.525,21535.950,51116.424,208703.225,20426.571,15874.618,955.600,25558.212,53285.185,529.691,57837.138,4551.953,1059.382,11322.664,13521.751,47623.853,27043.503,25661.305,51066.428,24498.830,37410.567,10213.286,39609.655,77204.327,19804.827),
                    t=x)))
    })

    L1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(2608814706223.0,1126008,303471,80538,21245,5592,1472,388,352,103,94,91,52,44,28,27),
                    B=c(0,6.2170397,3.055655,6.10455,2.83532,5.8268,2.5185,5.480,3.052,2.149,6.12,0.00,5.62,4.57,3.04,5.09),
                    C=c(0,26087.9031416,52175.806283,78263.70942,104351.61257,130439.5157,156527.4188,182615.322,1109.379,208703.225,27197.28,24978.52,5661.33,25028.52,51066.43,234791.13),
                    t=x)))
    })

    L2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(53050,16904,7397,3018,1107,378,123,39,15,12),
                    B=c(0,4.69072,1.3474,4.4564,1.2623,4.320,1.069,4.08,4.63,0.79),
                    C=c(0,26087.90314,52175.8063,78263.7094,104351.6126,130439.516,156527.419,182615.32,1109.8,208703.23),
                    t=x)))
    })

    L3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(188,142,97,44,35,18,7,3),
                    B=c(0.035,3.125,3.00,6.02,0,2.78,5.82,2.57),
                    C=c(52175.806,26087.903,78263.71,104351.61,0,130439.52,156527.42,182615.32),
                    t=x)))
    })

    L4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(114,3,2,2,1,1),
                    B=c(3.1416,2.03,1.42,4.50,4.50,1.27),
                    C=c(0,26087.90,78263.71,52175.81,104351.61,130439.52),
                    t=x)))
    })

    L5 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1),
                    B=c(3.14),
                    C=c(0),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((L0 + L1 * t + L2 * t^2 + L3 * t^3 + L4 * t^4 + L5 * t^5) / 10^8)))
  }
  if(celestial_body=="mars")
  {
    L0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(620347712,18656368,1108217,91798,27745,12316,10610,8927,8716,7775,6798,4161,3575,3075,2938,2628,2580,2389,1799,1546,1528,1286,1264,1025,892,859,833,833,749,724,713,655,636,553,550,472,426,415,312,307,302,299,293,284,281,274,274,239,236,231,221,204,193,189,179,174,172,160,144,140,138,131,128,128,117,113,110,105,100),
                    B=c(0,5.05037100,5.4009984,5.75479,5.97050,0.84956,2.93959,4.1570,6.1101,3.3397,0.3646,0.2281,1.6619,0.8570,6.0789,0.6481,0.0300,5.0390,0.6563,2.9158,1.1498,3.0680,3.6228,3.6933,0.183,2.401,4.495,2.464,3.822,0.675,3.663,0.489,2.922,4.475,3.810,3.625,0.554,0.497,0.999,0.381,4.486,2.783,4.221,5.769,5.882,0.542,0.134,5.372,5.755,1.282,3.505,2.821,3.357,1.491,1.006,2.414,0.439,3.949,1.419,3.326,4.301,4.045,2.208,1.807,3.128,3.701,1.052,0.785,3.243),
                    C=c(0,3340.61242670,6681.2248534,10021.83728,3.52312,2810.92146,2281.23050,0.0173,13362.4497,5621.8429,398.1490,2942.4634,2544.3144,191.4483,0.0673,3337.0893,3344.1355,796.2980,529.6910,1751.5395,6151.5339,2146.1654,5092.1520,8962.4553,16703.062,2914.014,3340.630,3340.595,155.420,3738.761,1059.382,3127.313,8432.764,1748.016,0.980,1194.447,6283.076,213.299,6677.702,6684.748,3532.061,6254.627,20.775,3149.164,1349.867,3340.545,3340.680,4136.910,3333.499,3870.303,382.897,1221.849,3.590,9492.146,951.718,553.569,5486.778,4562.461,135.065,2700.715,7.114,12303.068,1592.596,5088.629,7903.073,1589.073,242.729,8827.390,11773.377),
                    t=x)))
    })

    L1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(334085627474.0,1458227,164901,19963,3452,2485,842,538,521,433,430,382,314,283,206,169,158,134,134,118,117,114,114,91,85,83,81,80,73,73,71,68,65,65,62,57,48,48,47,41,40,40,33,28,27,27),
                    B=c(0,3.6042605,3.926313,4.26594,4.7321,4.6128,4.459,5.016,4.994,2.561,5.316,3.539,4.963,3.160,4.569,1.329,4.185,2.233,5.974,6.024,2.213,2.129,5.428,1.10,3.91,5.30,4.43,2.25,2.50,5.84,3.86,5.02,1.02,3.05,4.15,3.89,4.87,1.18,1.31,0.71,2.73,5.32,5.41,0.05,3.89,5.11),
                    C=c(0,3340.6124267,6681.224853,10021.83728,3.5231,13362.4497,2281.230,398.149,3344.136,191.448,155.420,796.298,16703.062,2544.314,2146.165,3337.089,1751.540,0.980,1748.016,6151.534,1059.382,1194.447,3738.761,1349.87,553.57,6684.75,529.69,8962.46,951.72,242.73,2914.01,382.90,3340.60,3340.63,3149.16,4136.91,213.30,3333.50,3185.19,1592.60,7.11,20043.67,6283.08,9492.15,1221.85,2700.72),
                    t=x)))
    })

    L2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(58016,54188,13908,2465,398,222,121,62,54,34,32,30,23,22,20,16,16,16,15,14,14,13,12,11,10,9,9,9,8,7,7,6,6),
                    B=c(2.04979,0,2.45742,2.8000,3.141,3.194,0.543,3.49,3.54,6.00,4.14,2.00,4.33,3.45,5.42,0.66,6.11,1.22,6.10,4.02,2.62,0.60,3.86,4.72,0.25,0.68,3.83,3.88,5.46,2.58,2.38,5.48,2.34),
                    C=c(3340.61243,0,6681.22485,10021.8373,13362.450,3.523,155.420,16703.06,3344.14,2281.23,191.45,796.30,242.73,398.15,553.57,0.98,2146.17,1748.02,3185.19,951.72,1349.87,1194.45,6684.75,2544.31,382.90,1059.38,20043.67,3738.76,1751.54,3149.16,4136.91,1592.60,3097.88),
                    t=x)))
    })

    L3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1482,662,188,41,26,23,10,8,5,4,3,3),
                    B=c(0.443,0.885,1.288,1.65,0,2.05,1.58,2.00,2.82,2.02,4.59,0.65),
                    C=c(3340.6124,6681.225,10021.837,13362.45,0,155.42,3.52,16703.06,242.73,3344.14,3185.19,553.57),
                    t=x)))
    })

    L4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(114,29,24,11,3,3,1,1),
                    B=c(3.1416,5.64,5.14,6.03,0.13,3.56,0.49,1.32),
                    C=c(0,6681.22,3340.61,10021.84,13362.45,155.42,16703.06,242.73),
                    t=x)))
    })

    L5 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1,1),
                    B=c(3.14,4.04),
                    C=c(0,6681.22),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((L0 + L1 * t + L2 * t^2 + L3 * t^3 + L4 * t^4 + L5 * t^5) / 10^8)))
  }
  if(celestial_body == "jupyter")
  {
    L0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(59954691,9695899,573610,306389,97178,72903,64264,39806,38858,27965,13590,8769,8246,7368,6263,6114,5305,5305,4905,4647,3045,2610,2028,1921,1765,1723,1633,1432,973,884,733,731,709,692,614,582,495,441,417,390,376,341,330,262,261,257,244,235,220,207,202,197,175,175,175,158,151,149,141,138,131,117,117,106),
                    B=c(0,5.0619179,1.444062,5.417347,4.14265,3.64043,3.41145,2.29377,1.27232,1.78455,5.77481,3.3600,3.5823,5.0810,0.0250,4.5132,4.1863,1.3067,1.3208,4.6996,4.3168,1.5667,1.0638,0.9717,2.1415,3.8804,3.5820,4.2968,4.098,2.437,6.85,3.806,1.293,6.134,4.109,4.540,3.756,2.958,1.036,4.897,4.703,5.715,4.740,1.877,0.820,3.724,5.220,1.227,1.651,1.855,1.807,5.293,3.730,3.226,5.910,4.365,3.906,4.377,3.136,1.318,4.169,2.500,3.389,4.554),
                    C=c(0,529.6909651,7.113547,1059.381930,632.78374,522.57742,103.09277,419.48464,316.39187,536.80451,1589.07290,949.1756,206.1855,735.8765,213.2991,1162.4747,1052.2684,14.2271,110.2063,3.9322,426.5982,846.0828,3.1814,639.8973,1066.4955,1265.5675,515.4639,625.6702,95.979,412.371,838.969,1581.959,742.990,2118.764,1478.867,309.278,323.505,454.909,2.448,1692.166,1368.660,533.623,0.048,0.963,380.128,199.072,728.763,909.819,543.918,525.759,1375.774,1155.361,942.062,1898.351,956.289,1795.258,74.782,1685.052,491.558,1169.588,1045.155,1596.186,0.521,526.510),
                    t=x)))
    })

    L1 <- sapply(t,function(x){
        sum(do.call(cos_component,
                    data.frame(
                      A=c(52993480757.0,489741,228919,27655,20721,12106,6068,5434,4238,2212,1746,1296,1173,1163,1099,1007,1004,848,827,816,725,568,474,413,345,336,234,234,199,195,187,184,171,131,115,115,108,80,72,70,67,66,65,59,58,57,57,55,52,52,50,47,47,40,34,33,32,29,29,29,25),
                      B=c(0,4.220667,6.026475,4.57266,5.45939,0.16986,4.4242,3.9848,5.8901,5.2677,4.9267,5.5513,5.8565,0.5145,5.3070,0.4648,3.1504,5.758,4.803,0.586,5.518,5.989,4.132,5.737,4.242,3.3732,4.035,6.243,1.505,2.219,6.086,6.280,5.417,0.626,0.680,5.286,4.493,5.82,5.34,5.97,5.73,0.13,6.09,0.59,0.99,5.97,1.41,5.43,5.73,0.23,6.08,3.63,0.51,4.16,0.10,5.04,5.37,5.42,3.36,0.76,1.61),
                      C=c(0,529.690965,7.113547,1059.38193,522.57742,536.80451,103.0928,419.4846,14.2271,206.1855,1589.0729,3.1814,1052.2684,3.9322,515.4639,735.8765,426.5982,110.206,213.299,1066.495,639.897,625.670,412.371,95.979,632.784,1162.475,949.176,309.278,838.969,323.505,742.990,543.918,199.072,728.763,846.083,2118.764,956.289,1045.15,942.06,532.87,21.34,526.51,1581.96,1155.36,1596.19,1169.59,533.62,10.29,117.32,1368.66,525.76,1478.87,1265.57,1692.17,302.16,220.41,508.35,1272.68,4.67,88.87,831.86),
                      t=x)))
    })

    L2 <- sapply(t,function(x){
          sum(do.call(cos_component,
                      data.frame(
                        A=c(47234,38966,30629,3189,2729,2723,1721,383,378,367,337,308,218,199,197,156,146,142,130,117,97,91,87,79,72,58,57,49,40,40,36,29,28,26,26,25,24,19,18,17,17,15,15,15,14,14,13,13,11,10,9,9,9,8,8,7,6),
                        B=c(4.32148,0,2.93021,1.0550,4.8455,3.4141,4.1873,5.768,0.760,6.055,3.786,0.694,3.814,5.340,2.484,1.406,3.814,1.634,5.837,1.414,4.03,1.11,2.52,4.64,2.22,0.83,3.12,1.67,4.02,0.62,2.33,3.61,3.24,4.50,2.51,1.22,3.01,4.29,0.81,4.20,1.83,5.81,0.68,4.00,5.95,1.80,2.52,4.37,4.44,1.72,2.18,3.29,3.32,5.76,2.71,2.18,0.50),
                        C=c(7.11355,0,529.69097,522.5774,536.8045,1059.3819,14.2271,419.485,515.464,103.093,3.181,206.186,1589.073,1066.495,3.932,1052.268,639.897,426.598,412.371,625.670,110.21,95.98,632.78,543.92,735.88,199.07,213.30,309.28,21.34,323.51,728.76,10.29,838.97,742.99,1162.47,1045.15,956.29,532.87,508.35,2118.76,526.51,1596.19,942.06,117.32,316.39,302.16,88.87,1169.59,525.76,1581.96,1155.36,220.41,831.86,846.06,533.62,1265.57,949.18),
                        t=x)))
    })

    L3 <- sapply(t,function(x){
            sum(do.call(cos_component,
                        data.frame(
                          A=c(6502,1357,471,417,353,155,87,44,34,28,24,23,20,20,19,17,17,16,16,13,13,13,9,9,7,7,6,5,5,4,4,4,3,3,3,3,3,2,2),
                          B=c(2.5986,1.3464,2.475,3.245,2.974,2.076,2.51,0,3.83,2.45,1.28,2.98,2.10,1.40,1.59,2.30,2.60,3.15,3.36,2.76,2.54,6.27,1.76,2.27,3.43,4.04,2.52,2.91,5.25,4.30,3.52,4.09,1.43,4.36,1.25,5.02,2.24,2.90,2.36),
                          C=c(7.1135,529.6910,14.227,536.805,522.577,1059.383,515.46,0,1066.50,206.19,412.37,543.92,639.90,419.48,103.09,21.34,1589.07,625.67,1052.27,95.98,199.07,426.60,10.29,110.21,309.28,728.76,508.35,1045.15,323.51,88.87,302.16,735.88,956.29,1596.19,213.30,838.97,117.32,742.99,942.06),
                          t=x)))
    })

    L4 <- sapply(t,function(x){
              sum(do.call(cos_component,
                          data.frame(
                            A=c(669,114,100,50,44,32,15,9,5,4,4,3,2,2,2,2,1,1,1),
                            B=c(0.853,3.142,0.743,1.65,5.82,4.86,4.29,0.71,1.30,2.32,0.48,3.00,0.40,4.26,4.91,4.26,5.26,4.72,1.29),
                            C=c(7.114,0,14.227,536.80,529.69,522.58,515.46,1059.38,543.92,1066.50,21.34,412.37,639.90,199.07,625.67,206.19,1052.27,95.98,1589.07),
                            t=x)))
    })

    L5 <- sapply(t,function(x){
                sum(do.call(cos_component,
                            data.frame(
                              A=c(50,16,4,2,1),
                              B=c(5.26,5.25,0.01,1.10,3.14),
                              C=c(7.11,14.23,536.80,522.58,0),
                              t=x)))
    })

    return (reduce_degrees(rad2degrees((L0 + L1 * t + L2 * t^2 + L3 * t^3 + L4 * t^4 + L5 * t^5) / 10^8)))
  }
  if(celestial_body=="saturn")
  {
    L0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(87401354,11107660,1414151,398379,350769,206816,79271,23990,16574,15820,15054,14907,14610,13160,13005,10725,6126,5863,5228,5020,4593,4006,3874,3269,2954,2461,1758,1640,1581,1391,1124,1087,1017,957,853,849,789,749,744,687,654,634,625,580,546,543,530,478,474,452,449,372,355,347,343,330,322,322,309,287,278,249,227,220,209,208,208,204,185,184,182,174,165,149,148,146,140,131,125,122,118,117,114,112,110,109,107,104,103,101),
                    B=c(0,3.96205090,4.5858152,0.521120,3.303299,0.246584,3.84007,4.66977,0.43719,0.93809,2.71670,5.76903,1.56519,4.44891,5.98119,3.12940,1.7633,0.2366,4.2078,3.1779,0.6198,2.2448,3.2228,0.7749,0.9828,2.0316,3.2658,5.5050,4.3727,4.0233,2.8373,4.1834,3.7170,0.507,3.421,3.191,5.007,2.144,5.253,1.747,1.599,2.299,0.970,3.093,2.127,1.518,4.449,2.965,5.475,1.044,1.290,2.278,3.013,1.539,0.246,0.247,0.961,2.572,3.495,2.370,0.400,1.470,4.910,4.204,1.345,0.483,1.283,6.011,3.503,0.973,5.491,1.863,0.440,5.736,1.535,6.231,4.295,4.068,6.277,1.976,5.341,2.679,5.594,1.105,0.166,3.438,4.012,2.192,1.197,4.965),
                    C=c(0,213.29909544,7.1135470,206.185548,426.598191,103.092774,220.41264,110.20632,419.48464,632.78374,639.89729,316.39187,3.93215,14.22709,11.04570,202.25340,277.0350,529.6910,3.1814,433.7117,199.0720,63.7359,138.5175,949.1756,95.9792,735.8765,522.5774,846.0828,309.2783,323.5054,415.5525,2.4477,227.5262,1265.567,175.166,209.367,0.963,853.196,224.345,1052.268,0.048,412.371,210.118,74.782,350.332,9.561,117.320,137.033,742.990,490.334,127.472,217.231,838.969,340.771,0.521,1581.959,203.738,647.011,216.480,351.817,211.815,1368.660,12.530,200.769,625.670,1162.475,39.357,265.989,149.563,4.193,2.921,0.751,5.417,52.690,5.629,195.140,21.341,10.295,1898.351,4.666,554.070,1155.361,1059.382,191.208,1.484,536.805,956.289,88.866,1685.052,269.921),
                    t=x)))
    })

    L1 <- sapply(t,function(x){
        sum(do.call(cos_component,
                    data.frame(
                      A=c(21354295596.0,1296855,564348,107679,98323,40255,19942,10512,6939,4803,4056,3769,3385,3302,3071,1953,1249,922,706,650,628,487,479,468,417,408,352,344,340,336,332,289,281,266,230,192,173,167,136,131,128,109,98,94,92,87,83,78,67,66,62,62,58,57,55,54,51,47,47,46,44,40,40,38,38,37,35,34,33,33,33,32,31,30,30,30,29,28,26),
                      B=c(0,1.8282054,2.885001,2.277699,1.08070,2.04128,1.27955,2.74880,0.4049,2.4419,2.9217,3.6497,2.4169,1.2626,2.3274,3.5639,2.6280,1.961,4.417,6.174,6.111,6.040,4.988,4.617,2.117,1.299,2.317,3.959,3.634,3.772,2.861,2.733,5.744,0.543,1.644,2.965,4.077,2.597,2.286,3.441,4.095,6.161,4.73,3.48,3.95,1.22,3.11,6.24,0.29,5.65,4.29,1.83,2.48,5.02,0.28,5.13,1.46,1.18,5.15,2.23,2.71,0.41,3.89,0.65,2.53,3.78,6.08,3.21,4.64,5.43,0.30,4.39,2.43,2.84,6.19,3.39,2.03,2.74,4.51),
                      C=c(0,213.2990954,7.113547,206.185548,426.59819,220.41264,103.09277,14.22709,639.8973,419.4846,110.2063,3.9322,3.1814,433.7117,199.0720,11.0457,95.9792,227.526,529.691,202.253,309.278,853.196,522.577,63.736,323.505,209.367,632.784,412.371,316.392,735.877,210.118,117.320,2.448,647.011,216.480,224.345,846.083,21.341,10.295,742.990,217.231,415.552,838.97,1052.27,88.87,440.83,625.67,302.16,4.67,9.56,127.47,195.14,191.96,137.03,74.78,490.33,536.80,149.56,515.46,956.29,5.42,269.92,728.76,422.67,12.53,2.92,5.63,1368.66,277.03,1066.50,351.82,1155.36,52.69,203.00,284.15,1059.38,330.62,265.99,340.77),
                      t=x)))
    })

    L2 <- sapply(t,function(x){
        sum(do.call(cos_component,
                    data.frame(
                      A=c(116441,91921,90592,15277,10631,10605,4265,1216,1165,1082,1045,1020,634,549,457,425,274,162,129,117,105,101,96,95,85,83,82,75,67,66,64,61,53,46,45,42,32,32,31,27,25,20,18,17,16,14,14,12,12,12,11,11,11,10,10,10,9,8,8,8,7,6,6),
                      B=c(1.179879,0.07425,0,4.06492,0.25778,5.40964,1.0460,2.9186,4.6094,5.6913,4.0421,0.6337,4.388,5.573,1.268,0.209,4.288,1.381,1.566,3.881,4.900,0.893,2.91,5.63,5.73,6.05,1.02,4.76,0.46,0.48,0.35,4.88,2.75,5.69,1.67,5.71,0.07,1.67,4.16,0.83,5.66,5.94,4.90,1.63,0.58,0.21,3.76,4.72,0.13,3.12,5.92,5.60,3.20,4.99,0.26,4.15,0.46,2.14,5.25,4.03,5.40,4.46,5.93),
                      C=c(7.113547,213.29910,0,206.18555,220.41264,426.59819,14.2271,103.0928,639.8973,433.7117,199.0720,3.1814,419.485,3.932,110.206,227.526,95.979,11.046,309.278,853.196,647.011,21.341,316.39,412.37,209.37,216.48,117.32,210.12,522.58,10.29,323.51,632.78,529.69,440.83,202.25,88.87,63.74,302.16,191.96,224.34,735.88,217.23,625.67,742.99,515.46,838.97,195.14,203.0,234.64,846.08,536.80,728.76,1066.50,422.67,330.62,860.31,956.29,269.92,429.78,9.56,1052.27,284.15,405.26),
                      t=x)))
    })

    L3 <- sapply(t,function(x){
          sum(do.call(cos_component,
                      data.frame(
                        A=c(10, 1067, 11, 11, 1162, 13, 131, 1466, 151, 16, 16, 16039, 166, 18, 18, 18, 18, 19, 1907, 2, 2, 2, 2, 2, 2, 2, 237, 239, 25, 28, 3, 3, 3, 3, 39, 4, 4, 40, 40, 4250, 5, 6, 6, 62, 63, 7, 8, 9),
                        B=c(5.73945,4.5854,4.7608,5.9133,5.6197,3.6082,3.861,5.768,5.116,2.736,4.743,0.23,4.74,5.47,5.96,5.83,3.01,0.99,1.92,4.97,1.03,4.20,3.32,3.90,5.62,1.18,5.58,5.93,3.95,3.39,4.88,0.38,2.25,1.06,4.64,3.14,2.31,2.20,0.59,4.93,0.42,4.77,3.35,3.20,1.19,1.35,4.16,3.07),
                        C=c(7.11355,213.2991,220.4126,206.1855,14.2271,426.5982,433.712,199.072,3.181,639.897,227.526,419.48,103.09,21.34,95.98,110.21,647.01,3.93,853.20,10.29,412.37,216.48,309.28,440.83,117.32,88.87,11.05,191.96,209.37,302.16,323.51,632.78,522.58,210.12,234.64,0,515.46,860.31,529.69,224.34,625.67,330.62,429.78,202.25,1066.50,405.26,223.59,654.12),
                        t=x)))
    })



    L4 <- sapply(t,function(x){
              sum(do.call(cos_component,
                          data.frame(
                            A=c(1662,257,236,149,114,110,68,40,38,31,15,9,6,6,4,4,3,3,3,3,3,2,2,2,2,2,1),
                            B=c(3.9983,2.984,3.902,2.741,3.142,1.515,1.72,2.05,1.24,3.01,0.83,3.71,2.42,1.16,1.45,2.12,4.09,2.77,3.01,0.00,0.39,3.78,2.83,5.08,2.24,5.19,1.55),
                            C=c(71135,220.413,14.227,213.299,0,206.186,426.60,433.71,199.07,227.53,639.90,21.34,419.48,647.01,95.98,440.83,110.21,412.37,88.87,853.20,103.09,117.32,234.64,309.28,216.48,302.16,191.96),
                            t=x)))
    })

    L5 <- sapply(t,function(x){
                sum(do.call(cos_component,
                            data.frame(
                              A=c(124,34,28,6,5,4,3,3,2,1,1,1),
                              B=c(2.259,2.16,1.20,1.22,0.24,6.23,2.97,4.29,6.25,5.28,0.24,3.14),
                              C=c(7.114,14.23,220.41,227.53,433.71,426.60,199.07,206.19,213.30,639.90,440.83,0),
                              t=x)))
    })

    return (reduce_degrees(rad2degrees((L0 + L1 * t + L2 * t^2 + L3 * t^3 + L4 * t^4 + L5 * t^5) / 10^8)))
  }
  if(celestial_body=="earth")
  {
    L0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(175347046,3341656,34894,3497,3418,3136,2676,2343,1324,1273,1199,990,902,857,780,753,505,492,357,317,284,271,243,206,205,202,156,132,126,115,103,102,102,99,98,86,85,85,80,79,75,74,74,70,62,61,57,56,56,52,52,51,49,41,41,39,37,37,36,36,33,30,30,25),
                    B=c(0,4.6692568,4.62610,2.7441,2.8289,3.6277,4.4181,6.1352,0.7425,2.0371,1.1096,5.233,2.045,3.508,1.179,2.533,4.583,4.205,2.920,5.849,1.899,0.315,0.345,4.806,1.869,2.458,0.833,3.411,1.083,0.645,0.636,0.976,4.267,6.21,0.68,5.98,1.30,3.67,1.81,3.04,1.76,3.50,4.68,0.83,3.98,1.82,2.78,4.39,3.47,0.19,1.33,0.28,0.49,5.37,2.40,6.17,6.04,2.57,1.71,1.78,0.59,0.44,2.74,3.16),
                    C=c(0,6283.0758500,12566.15170,5753.3849,3.5231,77713.7715,7860.4194,3930.2097,11506.7698,529.6910,1577.3435,5884.927,26.298,398.149,5223.694,5507.553,18849.228,775.523,0.067,11790.629,796.298,10977.079,5486.778,2544.314,5573.143,6069.777,213.299,2942.463,20.775,0.980,4694.003,15720.839,7.114,2146.17,155.42,161000.69,6275.96,71430.70,17260.15,12036.46,5088.63,3154.69,801.82,9437.76,8827.39,7084.90,6286.60,14143.50,6279.55,12139.55,1748.02,5856.48,1194.45,8429.24,19651.05,10447.39,10213.29,1059.38,2352.87,6812.77,17789.85,83996.85,1349.87,4690.48),
                    t=x)))
    })

    L1 <- sapply(t,function(x){
        sum(do.call(cos_component,
                    data.frame(
                      A=c(628331966747.0,206059,4303,425,119,109,93,72,68,67,59,56,45,36,29,21,19,19,17,16,16,15,12,12,12,12,11,10,10,9,9,8,6,6),
                      B=c(0,2.678235,2.6351,1.590,5.796,2.966,2.59,1.14,1.87,4.41,2.89,2.17,0.40,0.47,2.65,5.34,1.85,4.97,2.99,0.03,1.43,1.21,2.83,3.26,5.27,2.08,0.77,1.30,4.24,2.70,5.64,5.30,2.65,4.67),
                      C=c(0,6283.075850,12566.1517,3.523,26.298,1577.44,18849.23,529.69,398.15,5507.55,5223.69,155.42,796.30,775.52,7.11,0.98,5486.78,213.30,6275.96,2544.31,2146.17,10977.08,1748.02,5088.63,1194.45,4694.00,553.57,6286.60,1349.87,242.73,951.72,2352.87,9437.76,4690.48),
                      t=x)))
    })

    L2 <- sapply(t,function(x){
          sum(do.call(cos_component,
                      data.frame(
                        A=c(52919,8720,309,27,16,16,10,9,7,5,4,4,3,3,3,3,3,3,2,2),
                        B=c(0,1.0721,0.867,0.05,5.19,3.68,0.76,2.06,0.83,4.66,1.03,3.44,5.14,6.05,1.19,6.12,0.31,2.28,4.38,3.75),
                        C=c(0,6283.0758,12566.152,3.52,26.30,155.42,18849.23,77713.77,775.52,1577.34,7.11,5573.14,796.30,5507.55,242.73,529.69,398.15,553.57,5223.69,0.98),
                        t=x)))
    })

    L3 <- sapply(t,function(x){
            sum(do.call(cos_component,
                        data.frame(
                          A=c(289,35,17,3,1,1,1),
                          B=c(5.844,0,5.49,5.20,4.72,5.30,5.97),
                          C=c(6283.076,0,12566.15,155.42,3.52,18849.23,242.73),
                          t=x)))
    })

    L4 <- sapply(t,function(x){
              sum(do.call(cos_component,
                          data.frame(
                            A=c(114,8,1),
                            B=c(3.142,4.13,3.84),
                            C=c(0,6283.08,12566.15),
                            t=x)))
    })

    L5 <- sapply(t,function(x){
                sum(do.call(cos_component,
                            data.frame(
                              A=c(1),
                              B=c(3.14),
                              C=c(0),
                              t=x)))
    })

    return (reduce_degrees(rad2degrees((L0 + L1 * t + L2 * t^2 + L3 * t^3 + L4 * t^4 + L5 * t^5) / 10^8)))
  }
}

#' Heliocentric latitude of a planet
#'
#' @param jd Julian day
#' @param celestial_body planet name as string: mercury, venus,earth, mars, jupyter, saturn
#' @return Heliocentric longitude  in degrees
#' @examples
#' jd<-julian_day(13.19,11,2028)
#' heliocentric_latitude(jd,"venus")
#'

heliocentric_latitude <- function(jd,celestial_body) {
  t<-(jd - 2451545.0) / 365250

  if(celestial_body=="venus")
  {
    B0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(5923638,40108,32815,1011,149,138,130,120,108),
                    B=c(0.2670278,1.14737,3.14159,1.0895,6.254,0.860,3.672,3.705,4.539),
                    C=c(10213.285462,20426.57109,0,30639.8566,18073.705,1577.344,9437.763,2352.866,22003.915),
                    t=x)))
    })

    B1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(513348,4380,199,197),
                    B=c(1.803643,3.3862,0,2.530),
                    C=c(10213.285546,20426.5711,0,30639.857),
                    t=x)))
    })

    B2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(22378,282,173,27),
                    B=c(3.38509,0,5.256,3.87),
                    C=c(10213.2855,0,20426.571,30639.86),
                    t=x)))
    })

    B3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(647,20,6,3),
                    B=c(4.992,3.14,0.77,5.44),
                    C=c(10213.286,0,20426.57,30639.86),
                    t=x)))
    })

    B4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(14),
                    B=c(0.32),
                    C=c(10213.29),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((B0 + B1 * t + B2 * t^2 + B3 * t^3 + B4 * t^4) / 10^8)))
  }
  if(celestial_body=="mercury")
  {
    B0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(11737529,2388077,1222840,543252,129779,31867,7963,2014,514,209,208,132,121,100),
                    B=c(1.98357499,5.0373896,3.1415927,1.796444,4.832325,1.58088,4.6097,1.3532,4.378,2.020,4.918,1.119,1.813,5.657),
                    C=c(26087.90314157,52175.8062831,0,78263.709425,104351.612566,130439.51571,156527.4188,182615.3220,208703.225,24978.525,27197.282,234791.128,53285.185,20426.571),
                    t=x)))
    })

    B1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(429151,146234,22675,10895,6353,2496,860,278,86,28,26),
                    B=c(3.501698,3.141593,0.01515,0.48540,3.4294,0.1605,3.185,6.210,2.95,0.29,5.98),
                    C=c(26087.903142,0,52175.80628,78263.70942,104351.6126,130439.5157,156527.419,182615.322,208703.23,27197.28,234791.13),
                    t=x)))
    })

    B2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(11831,1914,1045,266,170,96,45,18,7),
                    B=c(4.79066,0,1.2122,4.434,1.623,4.80,1.6,4.67,1.43),
                    C=c(26087.90314,0,52175.8063,78263.709,104351.613,130439.52,156527.42,182615.32,208703.23),
                    t=x)))
    })

    B3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(235,161,19,6,5,3,2),
                    B=c(0.354,0,4.36,2.51,6.14,3.12,6.27),
                    C=c(26087.903,0,52175.81,78263.71,104351.61,130439.52,156527.42),
                    t=x)))
    })

    B4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(4,1),
                    B=c(1.75,3.14),
                    C=c(26087.90,0),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((B0 + B1 * t + B2 * t^2 + B3 * t^3 + B4 * t^4) / 10^8)))
  }
  if(celestial_body=="mars")
  {
    B0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(3197135,298033,289105,31366,3484,443,443,399,293,182,163,160,149,143,143,139),
                    B=c(3.7683204,4.106170,0,4.44651,4.7881,5.026,5.652,5.131,3.793,6.136,4.264,2.232,2.165,1.182,3.213,2.418),
                    C=c(3340.6124267,6681.224853,0,10021.83728,13362.4497,3344.136,3337.089,16703.062,2281.230,6151.534,529.691,1059.382,5621.843,3340.595,3340.630,8962.455),
                    t=x)))
    })

    B1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(350069,14116,9671,1472,426,102,79,33,26),
                    B=c(5.368478,3.14159,5.4788,3.2021,3.408,0.776,3.72,3.46,2.48),
                    C=c(3340.612427,0,6681.2249,10021.8373,13362.450,3337.089,16703.06,5621.84,2281.23),
                    t=x)))
    })

    B2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(16727,4987,302,26,21,12,8),
                    B=c(0.60221,3.1416,3.559,1.90,0.92,2.24,2.25),
                    C=c(3340.61243,0,6681.225,13362.45,10021.84,3337.09,16703.06),
                    t=x)))
    })

    B3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(607,43,14,3),
                    B=c(1.981,0,1.80,3.45),
                    C=c(3340.612,0,6681.22,10021.84),
                    t=x)))
    })

    B4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(13,11,1),
                    B=c(0,3.46,0.50),
                    C=c(0,3340.61,6681.22),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((B0 + B1 * t + B2 * t^2 + B3 * t^3 + B4 * t^4) / 10^8)))
  }
  if(celestial_body=="jupyter")
  {
    B0<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(2268616,110090,109972,8101,6438,6044,1107,944,942,894,836,767,684,629,559,532,464,431,351,132,123,116,115,104,103,102),
                    B=c(3.5585261,0,3.908093,3.6051,0.3063,4.2588,2.9853,1.675,2.936,1.754,5.179,2.155,3.678,0.643,0.014,2.703,1.173,2.608,4.611,4.778,3.350,1.387,5.049,3.701,2.319,3.153),
                    C=c(529.6909651,0,1059.381930,522.5774,536.8045,1589.0729,1162.4747,426.598,1052.268,7.114,103.093,632.784,213.299,1066.495,846.083,110.206,949.176,419.485,2118.764,742.990,1692.166,323.505,316.392,515.464,1478.867,1581.959),
                    t=x)))
    })

    B1<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(177352,3230,3081,2212,1694,346,234,196,150,114,97,82,77,77,74,61,50,46,45,37,36,32),
                    B=c(5.701665,5.7794,5.4746,4.7348,3.1416,4.746,5.189,6.186,3.927,3.439,2.91,5.08,2.51,0.61,5.50,5.45,3.95,0.54,1.90,4.70,6.11,4.92),
                    C=c(529.690965,1059.3819,522.5774,536.8045,0,1052.268,1066.495,7.114,1589.073,632.784,949.18,1162.47,103.09,419.48,515.46,213.30,735.88,110.21,846.08,543.92,316.39,1581.96),
                    t=x)))
    })

    B2<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(8094,813,742,399,342,74,46,30,29,23,14,12,11,6),
                    B=c(1.4632,3.1416,0.957,2.899,1.447,0.41,3.48,1.93,0.99,4.27,2.92,5.22,4.88,6.21),
                    C=c(529.6910,0,522.577,536.805,1059.382,1052.27,1066.50,1589.07,515.46,7.11,543.92,632.78,949.18,1045.15),
                    t=x)))
    })

    B3<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(252,122,49,11,8,7,6,4,9),
                    B=c(3.381,2.733,1.04,2.31,2.77,4.25,1.78,1.13,3.14),
                    C=c(529.691,522.577,536.80,1052.27,515.46,1059.38,1066.50,543.92,0),
                    t=x)))
    })

    B4<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(15,5,4,3,2,1),
                    B=c(4.53,4.47,5.44,0,4.52,4.20),
                    C=c(522.58,529.69,536.80,0,515.46,1052.27),
                    t=x)))
    })

    B5<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1),
                    B=c(0.09),
                    C=c(522.58),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((B0 + B1 * t + B2 * t^2 + B3 * t^3 + B4 * t^4 + B5 * t^5) / 10^8)))
  }
  if(celestial_body=="saturn")
  {
    B0<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(4330678,240348,84746,34116,30863,14734,9917,6994,4808,4788,3432,1506,1060,969,942,708,552,400,319,316,314,284,236,215,209,207,179,141,139,139,135,122,116,114),
                    B=c(3.6028443,2.852385,0,0.57297,3.48442,2.11847,5.7900,4.7360,5.4331,4.9651,2.7326,6.0130,5.6310,5.204,1.396,3.803,5.131,3.359,3.626,1.997,0.465,4.886,2.139,5.950,2.120,0.730,2.954,0.644,4.595,1.998,5.245,3.115,3.109,0.963),
                    C=c(213.2990954,426.598191,0,206.18555,220.41264,639.89729,419.4846,7.1135,316.3919,110.2063,433.7117,103.0928,529.6910,632.784,853.196,323.505,202.253,227.526,209.367,647.011,217.231,224.345,11.046,846.083,415.552,199.072,63.736,490.334,14.227,735.877,742.990,522.577,216.480,210.118),
                    t=x)))
    })

    B1<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(397555,49479,18572,14801,9644,3757,2717,1455,1291,853,298,292,284,275,172,166,158,128,110,82,81,69,65,61,59,46,36,34,33,32,27,27),
                    B=c(5.332900,3.14159,6.09919,2.30586,1.6967,1.2543,5.9117,0.8516,2.9177,0.436,0.919,5.316,1.619,3.889,0.052,2.444,5.209,1.207,2.457,2.76,2.86,1.66,1.26,1.25,1.82,0.82,1.82,2.84,1.31,1.19,4.65,4.44),
                    C=c(213.299095,0,426.59819,206.18555,220.4126,419.4846,639.8973,433.7117,7.1135,316.392,632.784,853.196,227.526,103.093,647.011,199.072,110.206,529.691,217.231,210.12,14.23,202.25,216.48,209.37,323.51,440.83,224.34,117.32,412.37,846.08,1066.50,11.05),
                    t=x)))
    })

    B2<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(20630,3720,1627,1346,706,365,330,219,139,104,93,71,52,49,41,29,24,21,20,18,17,16,14,12,8,7,7,6,6),
                    B=c(0.50482,3.9983,6.1819,0,3.039,5.099,5.279,3.828,1.043,6.157,1.98,4.15,2.88,4.43,3.16,4.53,1.12,4.35,5.31,0.85,5.68,4.26,3.00,2.53,3.32,5.56,0.29,1.16,3.61),
                    C=c(213.29910,206.1855,220.4126,0,419.485,426.598,433.712,639.897,7.114,227.526,316.39,199.07,632.78,647.01,853.20,210.12,14.23,217.23,440.83,110.21,216.48,103.09,412.37,529.69,202.25,209.37,323.51,117.32,860.31),
                    t=x)))
    })

    B3<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(666,632,398,188,92,52,42,26,21,18,11,10,7,6,6,6,5,5,4,3,2),
                    B=c(1.990,5.698,0,4.338,4.84,3.42,2.38,4.40,5.85,1.99,5.37,2.55,3.46,4.80,0.02,3.52,5.64,1.22,4.71,0.63,3.72),
                    C=c(213.299,206.186,0,220.413,419.48,433.71,426.60,227.53,199.07,639.90,7.11,647.01,316.39,632.78,210.12,440.83,14.23,853.20,412.37,103.09,216.48),
                    t=x)))
    })

    B4<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(80,32,17,12,9,6,5,5,1,1,1,1),
                    B=c(1.12,3.12,2.48,3.14,0.38,1.56,2.63,1.28,1.43,0.67,1.72,6.18),
                    C=c(206.19,213.30,220.41,0,419.48,433.71,227.53,199.07,426.60,647.01,440.83,639.90),
                    t=x)))
    })

    B5<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(8,1),
                    B=c(2.82,0.51),
                    C=c(206.19,220.41),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((B0 + B1 * t + B2 * t^2 + B3 * t^3 + B4 * t^4 + B5 * t^5) / 10^8)))
  }
  if(celestial_body=="earth")
  {
    B0<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(280.0,102,80,44,32),
                    B=c(3.199,5.422,3.88,3.70,4.00),
                    C=c(84334.662,5507.553,5223.69,2352.87,1577.34),
                    t=x)))
    })

    B1<- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(9,6),
                    B=c(3.90,1.73),
                    C=c(5507.55,5223.69),
                    t=x)))
    })

    return (reduce_degrees(rad2degrees((B0 + B1 * t) / 10^8)))
  }
}

#' Heliocentric radius of a planet
#'
#' @param jd Julian day
#' @param celestial_body planet name as string: mercury, venus,earth, mars, jupyter, saturn
#' @return Heliocentric radius
#' @examples
#' jd<-julian_day(13.19,11,2028)
#' heliocentric_radius(jd,"venus")
#'

heliocentric_radius <- function(jd,celestial_body) {
  t<-(jd - 2451545.0) / 365250
  if(celestial_body=="venus")
  {
    R0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(72334821,489824,1658,1632,1378,498,374,264,237,222,126,119),
                    B=c(0,4.021518,4.9021,2.8455,1.1285,2.587,1.423,5.529,2.551,2.013,2.728,3.020),
                    C=c(0,10213.285546,20426.5711,7860.4194,11790.6291,9683.595,3930.210,9437.763,15720.839,19367.189,1577.344,10404.734),
                    t=x)))
    })

    R1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(34551,234,234),
                    B=c(0.89199,1.772,3.142),
                    C=c(10213.28555,20426.571,0),
                    t=x)))
    })

    R2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1407,16,13),
                    B=c(5.06337,5.47,0),
                    C=c(10213.2855,20426.57,0),
                    t=x)))
    })

    R3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(50),
                    B=c(3.22),
                    C=c(10213.29),
                    t=x)))
    })

    R4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1),
                    B=c(0.92),
                    C=c(10213.29),
                    t=x)))
    })

    return ((R0 + R1 * t + R2 * t^2 + R3 * t^3 + R4 * t^4) / 10^8)
  }
  if(celestial_body=="mercury")
  {
    R0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(39528272,7834132,795526,121282,21922,4354,918,290,260,202,201,142,100),
                    B=c(0,6.1923372,2.959897,6.010642,2.77820,5.8289,2.597,1.424,3.028,5.647,5.592,6.253,3.734),
                    C=c(0,26087.9031416,52175.806283,78263.709425,104351.61257,130439.5157,156527.419,25028.521,27197.282,182615.322,31749.235,24978.525,24535.950),
                    t=x)))
    })

    R1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(217348,44142,10094,2433,1624,604,153,39),
                    B=c(4.656172,1.42386,4.47466,1.2423,0,4.293,1.061,4.11),
                    C=c(26087.903142,52175.80628,78263.70942,104351.6126,0,130439.516,156527.419,182615.32),
                    t=x)))
    })

    R2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(3118,1245,425,136,42,22,13),
                    B=c(3.0823,6.1518,2.926,5.980,2.75,3.14,5.80),
                    C=c(26087.9031,52175.8063,78263.709,104351.613,130439.52,0,156527.42),
                    t=x)))
    })

    R3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(33,24,12,5,2),
                    B=c(1.68,4.63,1.39,4.44,1.21),
                    C=c(26087.90,52175.81,78263.71,104351.61,130439.52),
                    t=x)))
    })


    return ((R0 + R1 * t + R2 * t^2 + R3 * t^3) / 10^8)
  }
  if(celestial_body=="mars")
  {
    R0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(153033488,14184953,660776,46179,8110,7485,5523,3825,2484,2307,1999,1960,1167,1103,992,899,807,798,741,726,692,633,633,630,574,526,473,348,284,280,276,275,270,239,234,228,223,219,208,208,186,183,179,176,164),
                    B=c(0,3.47971284,3.817834,4.15595,5.5596,1.7724,1.3644,4.4941,4.9255,0.0908,5.3606,4.7425,2.1126,5.0091,5.839,4.408,2.102,3.448,1.499,1.245,2.134,0.894,2.924,1.287,0.829,5.383,5.199,4.832,2.907,5.257,1.218,2.908,3.764,2.037,5.105,3.255,4.199,5.583,5.255,4.846,5.699,5.081,4.184,5.953,3.799),
                    C=c(0,3340.61242670,6681.224853,10021.83728,2810.9215,5621.8429,2281.2305,13362.4497,2942.4634,2544.3144,3337.0893,3344.1355,5092.1520,398.1490,6151.534,529.691,1059.382,796.298,2146.165,8432.764,8962.455,3340.595,3340.630,1751.540,2914.014,3738.761,3127.313,16703.062,3532.061,6283.076,6254.627,1748.016,5884.927,1194.447,5486.778,6872.673,3149.164,191.448,3340.545,3340.680,6677.702,6684.748,3333.499,3870.303,4136.910),
                    t=x)))
    })

    R1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1107433,103176,12877,10816,1195,439,396,183,136,128,128,127,118,88,83,76,72,67,66,58,54,51,49,49,48,48,39),
                    B=c(2.0325052,2.370718,0,2.70888,3.0470,2.888,3.423,1.584,3.385,6.043,0.630,1.954,2.998,3.42,3.86,4.45,2.76,2.55,4.41,0.54,0.68,3.73,5.73,1.48,2.58,2.29,2.32),
                    C=c(3340.6124267,6681.224853,0,10021.83728,13362.4497,2281.230,3344.136,2544.314,16703.062,3337.089,1059.382,796.298,2146.165,398.15,3738.76,6151.53,529.69,1751.54,1748.02,1194.45,8962.46,6684.75,3340.60,3340.63,3149.16,2914.01,4136.91),
                    t=x)))
    })

    R2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(44242,8138,1275,187,52,41,27,18,12,10,10),
                    B=c(0.47931,0.8700,1.2259,1.573,3.14,1.97,1.92,4.43,4.53,5.39,0.42),
                    C=c(3340.61243,6681.2249,10021.8373,13362.450,0,3344.14,16703.06,2281.23,3185.19,1059.38,796.30),
                    t=x)))
    })

    R3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1113,424,100,20,5,3),
                    B=c(5.1499,5.613,5.997,0.08,3.14,0.43),
                    C=c(3340.6124,6681.225,10021.837,13362.45,0,16703.06),
                    t=x)))
    })

    R4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(20,16,6,2),
                    B=c(3.58,4.05,4.46,4.84),
                    C=c(3340.61,6681.22,10021.84,13362.45),
                    t=x)))
    })

    return ((R0 + R1 * t + R2 * t^2 + R3 * t^3 + R4 * t^4) / 10^8)
  }
  if(celestial_body=="jupyter")
  {
    R0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(520887429,25209327,610600,282029,187647,86793,72063,65517,30135,29135,23947,23453,22284,13033,12749,9703,9161,7895,7058,6138,5477,4170,4137,3503,2617,2500,2128,1912,1611,1479,1231,1217,1015,999,961,886,821,812,777,727,655,654,621,615,562,542),
                    B=c(0,3.49108640,3.841154,2.574199,2.075904,0.71001,0.21466,5.97996,2.16132,1.67759,0.27458,3.54023,4.19363,2.96043,2.71550,1.9067,4.4135,2.4791,2.1818,6.2642,5.6573,2.0161,2.7222,0.5653,2.0099,4.5518,6.1275,0.8562,3.0887,2.6803,1.8904,1.8017,1.3867,2.872,4.549,4.148,1.593,5.941,3.677,3.988,2.791,3.382,4.823,2.276,0.081,0.284),
                    C=c(0,529.69096509,1059.381930,632.783739,522.577418,419.48464,536.80451,316.39187,949.17561,103.09277,7.11355,735.87651,1589.07290,1162.47470,1052.26838,206.1855,213.2991,426.5982,1265.5675,846.0828,639.8973,515.4639,625.6702,1066.4955,1581.9593,838.9693,742.9901,412.3711,1368.6603,1478.8666,323.5054,110.2063,454.9094,309.278,2118.764,533.623,1898.351,909.819,728.763,1155.361,1685.052,1692.166,956.289,942.062,543.918,525.759),
                    t=x)))
    })

    R1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1271802,61662,53444,41390,31185,11847,9166,3404,3203,3176,2806,2677,2600,2412,2101,1646,1641,1050,1025,806,741,677,567,485,469,445,416,402,347,338,261,247,220,203,200,197,197,184,180,170,146,133,132),
                    B=c(2.6493751,3.0076,3.89718,0,4.88277,2.41330,4.7598,3.3469,5.2108,2.7930,3.7422,4.3305,3.6344,1.4695,3.9276,5.3095,4.4163,3.1611,2.5543,2.678,2.171,6.250,4.577,2.469,4.710,0.403,5.368,4.605,4.681,3.168,5.343,3.923,4.842,5.600,4.439,3.706,3.759,4.265,4.402,4.846,6.130,1.322,4.512),
                    C=c(529.6909651,1059.38193,522.57742,0,536.80451,419.48464,7.1135,1589.0729,735.8765,103.0928,515.4639,1052.2684,206.1855,426.5982,639.8973,1066.4955,625.6702,213.2991,412.3711,632.784,1162.475,838.969,742.990,949.176,543.918,323.505,728.763,309.278,14.227,956.289,846.083,942.062,1368.660,1155.361,1045.155,2118.764,199.072,95.979,532.872,526.510,533.623,110.206,525.759),
                    t=x)))
    })

    R2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(79645,8252,7030,5314,1861,964,836,498,427,406,377,363,342,339,333,280,257,230,201,200,139,114,95,86,83,80,75,70,67,62,56,52,50,45,44,40),
                    B=c(1.35866,5.7777,3.2748,1.8384,2.9768,5.480,4.199,3.142,2.228,3.783,2.242,5.368,6.099,6.127,0.003,4.262,0.963,0.705,3.069,4.429,2.932,0.787,1.70,5.14,0.06,2.98,1.60,1.51,5.47,6.10,0.96,5.58,2.72,5.52,0.27,5.95),
                    C=c(529.69097,522.5774,536.8045,1059.3819,7.1135,515.646,419.485,0,639.897,1066.495,1589.073,206.186,1052.268,625.670,426.598,412.371,632.784,735.877,543.918,103.093,14.227,728.763,838.97,323.51,309.28,742.99,956.29,213.30,199.07,1045.15,1162.47,942.06,532.87,508.35,526.51,95.98),
                    t=x)))
    })

    R3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(3519,1073,916,342,255,222,90,69,58,58,51,47,43,37,34,34,31,30,21,15,14,13,12,12,11,11,10,9),
                    B=c(6.058,1.6732,1.413,0.523,1.196,0.952,3.14,2.27,1.41,0.53,5.98,1.58,6.12,1.18,1.67,0.85,1.04,4.63,2.50,0.89,0.96,1.50,2.61,3.56,1.79,6.28,6.26,3.45),
                    C=c(529.6910,536.8045,522.577,1059.382,7.114,515.464,0,1066.50,543.92,639.90,412.37,625.67,419.48,14.23,1052.27,206.19,1589.07,426.60,728.76,199.07,508.35,1045.15,735.88,323.51,309.28,956.29,103.09,838.97),
                    t=x)))
    })

    R4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(129,113,83,38,27,18,13,9,8,7,6,5,3,3,3),
                    B=c(0.084,4.249,3.30,2.73,5.69,5.40,6.02,0.77,5.68,1.43,5.12,3.34,3.40,4.16,2.90),
                    C=c(536.805,529.691,522.58,515.46,7.11,1059.38,543.92,1066.50,14.23,412.37,639.90,625.67,1052.27,728.76,426.60),
                    t=x)))
    })

    R5 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(11,4,2,2,2,2,2),
                    B=c(4.75,5.92,5.57,4.30,3.69,4.13,5.49),
                    C=c(536.80,522.58,515.46,543.92,7.11,1059.38,1066.50),
                    t=x)))
    })

    return ((R0 + R1 * t + R2 * t^2 + R3 * t^3 + R4 * t^4 + R5 * t^5) / 10^8)
  }
  if(celestial_body=="saturn")
  {
    R0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(955758136,52921382,1873680,1464664,821891,547507,371684,361778,140618,108975,69007,61053,48913,34144,32402,20937,20839,20747,15298,14296,12884,11993,11380,9796,7753,6771,6466,5850,5307,4696,4044,3688,3461,3420,3401,3376,2976,2885,2881,2508,2448,2406,2174,2024),
                    B=c(0,2.39226220,5.2354961,1.6476305,5.935200,5.015326,2.271148,3.139043,5.704067,3.293136,5.94100,0.94038,1.55733,0.19519,5.47085,0.46349,1.52103,5.33256,3.05944,2.60434,1.64892,5.98051,1.73106,5.2048,5.8519,3.0043,0.1773,1.4552,0.5974,2.1492,1.6401,0.7802,1.8509,4.9455,0.5539,3.6953,5.6847,1.3876,0.1796,3.5385,6.1841,2.9656,0.0151,5.0541),
                    C=c(0,213.29909544,206.1855484,426.5981909,316.391870,103.092774,220.412642,7.113547,632.783739,110.206321,419.48464,639.89729,202.25340,277.03499,949.17561,735.87651,433.71174,199.07200,529.69097,323.50542,138.51750,846.08283,522.57742,1265.5675,95.9792,14.2271,1052.2684,415.5525,63.7359,227.5262,209.3669,412.3711,175.1661,1581.9593,350.3321,224.3448,210.1177,838.9693,853.1964,742.9901,1368.6603,117.3199,340.7709,11.0457),
                    t=x)))
    })

    R1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(6182981,506578,341394,188491,186262,143891,49621,20928,19953,18840,13877,12893,5397,4869,4247,3252,3081,2909,2856,1988,1941,1581,1340,1316,1203,1091,966,954,898,882,874,785,740,658,650,613,599,503),
                    B=c(0.2584352,0.711147,5.796358,0.472157,3.141593,1.407449,6.01744,5.09246,1.17560,1.60820,0.75886,5.94330,1.2885,0.8679,0.3930,1.2585,3.4366,4.6068,2.1673,2.4505,6.0239,1.2919,4.3080,1.2530,1.8665,0.0753,0.480,5.152,0.983,1.885,1.402,3.064,1.382,4.144,1.725,3.033,2.549,2.130),
                    C=c(213.2990954,206.185548,426.598191,220.412642,0,7.113547,103.09277,639.89729,419.48464,110.20632,199.07200,433.71174,14.2271,323.5054,227.5262,95.9792,522.5774,202.2534,735.8765,412.3711,209.3669,210.1177,853.1964,117.3199,316.3919,216.4805,632.784,647.011,529.691,1052.268,224.345,838.969,625.670,309.278,742.990,63.736,217.231,3.932),
                    t=x)))
    })

    R2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(436902,71923,49767,43221,29646,4721,4142,3789,2964,2556,2327,2208,2188,1957,924,706,546,431,405,391,374,361,356,326,207,204,180,178,154,148,133,132),
                    B=c(4.786717,2.50070,4.97168,3.86940,5.96310,2.4753,4.1067,3.0977,1.3721,2.8507,0,6.2759,5.8555,4.9245,5.464,2.971,4.129,5.178,4.173,4.481,5.834,3.277,3.192,2.269,4.022,0.088,3.597,4.097,3.135,0.136,2.594,5.933),
                    C=c(213.299095,206.18555,220.41264,426.59819,7.11355,199.0720,433.7117,639.8973,103.09228,419.4846,0,110.2063,14.2271,227.5262,323.505,95.979,412.371,522.577,209.367,216.480,117.320,647.011,210.118,853.196,735.877,202.253,632.784,440.825,625.670,302.165,191.958,309.278),
                    t=x)))
    })

    R3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(20315,8924,6909,4087,3879,1071,907,606,597,483,393,229,188,150,121,102,101,93,84,73,62,55,50,45,41,40,38,32),
                    B=c(3.02187,3.1914,4.3517,4.2241,2.0106,4.2036,2.283,3.175,4.135,1.173,0,4.698,4.590,3.202,3.768,4.710,5.819,1.44,2.63,4.15,2.31,0.31,2.39,4.37,0.69,1.84,5.94,4.01),
                    C=c(213.29910,220.4126,206.1855,7.1135,426.5982,199.0720,433.712,227.526,14.227,639.897,0,419.485,110.206,103.093,323.505,95.979,412.371,647.01,216.48,117.32,440.83,853.20,209.37,191.96,522.58,302.16,88.87,21.34),
                    t=x)))
    })

    R4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(1202,708,516,427,268,170,150,145,121,47,19,17,16,15,14,13,11,11,10,9,9,9,8),
                    B=c(1.4150,1.162,6.240,2.469,0.187,5.959,0.480,1.442,2.405,5.57,5.86,0.53,2.90,0.30,1.30,2.09,0.22,2.46,3.14,1.56,2.28,0.68,1.27),
                    C=c(220.4126,213.299,206.186,7.114,426.598,199.072,433.712,227.526,14.227,639.90,647.01,440.83,110.21,419.48,412.37,323.51,95.98,117.32,0,88.87,21.34,216.48,234.64),
                    t=x)))
    })

    R5 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(129,32,27,20,20,14,14,13,7,5,4,3,3,3,3,2,2,2),
                    B=c(5.913,0.69,5.91,4.95,0.67,2.67,1.46,4.59,4.630,3.61,4.90,4.07,4.66,0.49,3.18,3.70,3.32,0.56),
                    C=c(220.413,7.11,227.53,433.71,14.23,206.19,199.07,426.60,213.30,639.90,440.83,647.01,191.96,323.51,419.48,88.87,95.98,117.32),
                    t=x)))
    })

    return ((R0 + R1 * t + R2 * t^2 + R3 * t^3 + R4 * t^4 + R5 * t^5) / 10^8)
  }
  if(celestial_body=="earth")
  {
    R0 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(100013989,1670700,13956,3084,1628,1576,925,542,472,346,329,307,243,212,186,175,110,98,86,86,65,63,57,56,49,47,45,43,39,38,37,37,36,35,33,32,32,28,28,26),
                    B=c(0,3.0984635,3.05525,5.1985,1.1739,2.8469,5.453,4.564,3.661,0.964,5.900,0.299,4.273,5.847,5.022,3.012,5.055,0.89,5.69,1.27,0.27,0.92,2.01,5.24,3.25,2.58,5.54,6.01,5.36,2.39,0.83,4.90,1.67,1.84,0.24,0.18,1.78,1.21,1.90,4.59),
                    C=c(0,6283.0758500,12566.15170,77713.7715,5753.3849,7860.4194,11506.770,3930.210,5884.927,5507.553,5223.694,5573.143,11790.629,1577.344,10977.079,18849.228,5486.778,6069.78,15720.84,161000.69,17260.15,529.69,83996.85,71430.70,2544.31,775.52,9437.76,6275.96,4694.00,8827.39,19651.05,12139.55,12036.46,2942.46,7084.90,5088.63,398.15,6286.60,6279.55,10447.39),
                    t=x)))
    })

    R1 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(103019,1721,702,32,31,25,18,10,9,9),
                    B=c(1.107490,1.0644,3.142,1.02,2.84,1.32,1.42,5.91,1.42,0.27),
                    C=c(6283.075850,12566.1517,0,18849.23,5507.55,5223.69,1577.34,10977.08,6275.96,5486.76),
                    t=x)))
    })

    R2 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(4359,124,12,9,6,3),
                    B=c(5.7846,5.579,3.14,3.63,1.87,5.47),
                    C=c(6283.0758,12566.152,0,77713.77,5573.14,18849.23),
                    t=x)))
    })

    R3 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(145,7),
                    B=c(4.273,3.92),
                    C=c(6283.076,12566.15),
                    t=x)))
    })

    R4 <- sapply(t,function(x){
      sum(do.call(cos_component,
                  data.frame(
                    A=c(4),
                    B=c(2.56),
                    C=c(6283.08),
                    t=x)))
    })

    return ((R0 + R1 * t + R2 * t^2 + R3 * t^3 + R4 * t^4) / 10^8)
  }
}