const unsigned int st[51] = {
  69754, //19:22:34
   4526, //01:15:26
  25493, //07:04:53
  45761, //12:42:41
  66961, //18:36:01
   1390, //00:23:10
  22134, //06:08:54
  43491, //12:04:51
  64073, //17:47:53
  85173, //23:39:33
  19869, //05:31:09
  40363, //11:12:43
  61927, //17:12:07
  83048, //23:04:08
  17344, //04:49:04
  38718, //10:45:18
  59345, //16:29:05
  80633, //22:23:53
  15634, //04:20:34
  36209, //10:03:29
  57626, //16:00:26
  78561, //21:49:21
  12510, //03:28:30
  33702, //09:21:42
  54253, //15:04:13
  75082, //20:51:22
   9797, //02:43:17
  30048, //08:20:48
  51314, //14:15:14
  72646, //20:10:46
   7001, //01:56:41
  28624, //07:57:04
  49628, //13:47:08
  70508, //19:35:08
   5521, //01:32:01
  26041, //07:14:01
  47333, //13:08:53
  68608, //19:03:28
   2505, //00:41:45
  23640, //06:34:00
  44370, //12:19:30
  65114, //18:05:14
    145, //00:02:25
  20687, //05:44:47
  41779, //11:36:19
  62982, //17:29:42
  83308, //23:08:28
  18210, //05:03:30
  39206, //10:53:26
  60000  //16:40:00
};

double MeanObliquityOfEcliptic(double JD)		/* number of centuries since J2000 */
{
    //static double OblPoly[11] = { 84381.448, -4680.93, -1.55, 1999.25, -51.38, -249.67, -39.05, 7.12, 27.87, 5.79, 2.45};
    double T = (JD - 2451545) / 36525;
    double u = T / 100;
    val = (0.130218283951 - u*(0.00722365740741  + u*(0.239197530864197e-05 + u*(-0.003085262345679 + u*(0.00007929012345679 + u*(0.00038529321 + u*(0.000060262345679
                          - u*(0.000010987654321 + u*(0.00004300925925925925 + u*(0.0000089351851855185 + 0.000003780864197531*u))))))))));
    return val;
}

double EarthMeanLongitude(double JD)
{
  double T = (JD - 2451545) / 36525;
  /* Solar mean anomaly, (47.3)  */
  double M = M_PI*(1.98627282889 + T*(199.9947238389 - 0.853333333e-06*T));
  /* Solar mean longitude, (25.2) */
  double L0 = (1.558147 + 200.004276839*T);
  double F = (0.51817830555556 + T*(2684.45565290722222 - T*(0.0002029944445 + 0.0000000157559*T)));
  double C = 0.03341543385*sin(M) + 0.000348972895596125*sin(2*M) + 0.00000998187972*sin(M_PI*F);
  double solar_lon = L0 + C;
  return (solar_lon - lrint(solar_lon)) - 0.5;
}

double NorthernSolstice(int Year)
{
  //calculate the approximate date
  double Y = (Year - 2000) / 100.0;
  double Ysquared = Y*Y;
  return 2451717.06767 + 36524.162603*Y + 0.0000325*Ysquared + Ysquared*(0.00000889*Y - 0.000000030*Ysquared);
}

double SouthernSolstice(int Year)
{
  //calculate the approximate date
  double Y = (Year - 2000) / 100.0;
  double Y2 = Y*Y;
  return 2451900.55952 + 36524.274049*Y - 0.006223*Y2 + Y2*(-0.00000823*Y + 0.000000032*Y2);
}