picplanner/src/calculations.c

#include <calculations.h>

/* TODO
 * Azimuth and Elevation are wrong in the calculation of the position of the Milky Way */


/* Convert degree to radiant */
double calc_deg_rad (double deg) {
  return M_PI/180*deg;
}
/* Convert radiant to degree */
double calc_rad_deg (double rad) {
  return 180/M_PI*rad;
}

int times_to_zone (int day_utc, int hour_utc, int day_local, int hour_local){
  int time_zone;
  if (day_utc == day_local){
    time_zone = hour_local - hour_utc;
  }
  else{
    if (abs(day_utc-day_local)>1){
      if (day_utc < day_local){
        day_local = day_utc-1;
      }
      else{
        day_local = day_utc+1;
      }
    }
    time_zone = 24*(day_local-day_utc)+(hour_local-hour_utc);
  }
  return time_zone;
}

void utc_zone_to_time (int *time_utc, int *time_local){
  int sign = *(time_utc+5)/abs(*(time_utc+5));

  *time_local=*time_utc;
  *(time_local+1)=*(time_utc+1);
  *(time_local+2)=*(time_utc+2);
  *(time_local+3)=*(time_utc+3);
  *(time_local+4)=*(time_utc+4);
  *(time_local+5)=*(time_utc+5);

  *(time_local+3) = *(time_utc+3)+*(time_utc+5);
  if (*(time_local+3)>24 || *(time_local+3) < 0){ /* Test if hours are still inside the same day */
    *(time_local+3) = (*(time_local+3)+24)%24;
    *(time_local+2) = *(time_utc+2)+sign;

    if (*(time_local+2) > max_day_in_month (*(time_utc+1), *(time_utc)) || *(time_local+2) < 1){ /* Test if days are still inside the same month */
      *(time_local+2) = (*(time_local+2)-1)%max_day_in_month (*(time_utc+1)+sign, *(time_utc))+1; /* -1 and +1 necessary due to mod calculations */
      *(time_local+1) = *(time_utc+1)+sign;

      if (*(time_local+1) > 12 || *(time_local+1) < 1){ /* Test if months are still inside the same year */
        *(time_local+1) = (*(time_local+1)-1)%12+1;
        *(time_local) = *(time_utc)+sign;
      }
    }
  }
}

void time_to_utc (int *time_local, int *time_utc){
  int sign = *(time_utc+5)/abs(*(time_utc+5));

  *time_utc=*time_local;
  *(time_utc+1)=*(time_local+1);
  *(time_utc+2)=*(time_local+2);
  *(time_utc+3)=*(time_local+3);
  *(time_utc+4)=*(time_local+4);
  *(time_local+5)=*(time_utc+5);

  *(time_utc+3) = *(time_local+3)-*(time_utc+5);
  if (*(time_utc+3) > 24 || *(time_utc+3) < 0){
    *(time_utc+3) = (*(time_utc+3)+24)%24;
    *(time_utc+2) = *(time_utc+2)+sign;

    if (*(time_utc+2) > max_day_in_month (*(time_local+1), *(time_local)) || *(time_utc+2) < 1){
      *(time_utc+2) = (*(time_utc+2)-1)%max_day_in_month (*(time_local+1)+sign, *(time_local))+1;
      *(time_utc+1) = *(time_local+1)+sign;

      if (*(time_utc+1) > 12 || *(time_utc+1) < 1){
        *(time_utc+1) = (*(time_local+1)-1)%12+1;
        *(time_utc) = *(time_local)+sign;
      }
    }
  }
}

int max_day_in_month (int month, int year){
  int max_days = 30;
  /* check which month is set */
  switch (month) {
    case 1: case 3: case 5: case 7: case 8: case 10: case 12:
      max_days = 31;
      break;
    case 2:
      /* definition of a leap year: Every fourth year is a leap year expect if year%100 is 0 e.g. 1900
       * but every year where year%400 is 0 is a leap year again e.g. 2000 */
      if ((year%4==0 && year%100!=0) || year%400==0) {
        max_days = 29;
      }
      else {
        max_days = 28;
      }
      break;
  }
  return max_days;
}


/* Get current time and date and return as a pointer to an array */
void calc_time_utc (int *time_utc)
{
  time_t t;
  struct tm tmp_local;
  struct tm *tmp;

  t = time(&t);

  tmp = localtime(&t);
  memcpy(&tmp_local, tmp, sizeof(struct tm));
  tmp = gmtime(&t);

  int time_zone = times_to_zone (tmp->tm_mday, tmp->tm_hour, tmp_local.tm_mday, tmp_local.tm_hour);

  g_print ("UTC: %02d:%02d - %d.%d.%d, Time zone: %d\n", tmp->tm_hour, tmp->tm_min, tmp->tm_mday, tmp->tm_mon+1, tmp->tm_year+1900, time_zone);

  *(time_utc) = tmp->tm_year+1900;
  *(time_utc+1) = tmp->tm_mon+1;
  *(time_utc+2) = tmp->tm_mday;
  *(time_utc+3) = tmp->tm_hour;
  *(time_utc+4) = tmp->tm_min;
  *(time_utc+5) = time_zone;
}


/* Calculate the Julian Day */
void calc_jd (int *time_utc, double *time_jd)
{
  float year = *time_utc;
  float month = *(time_utc+1);
  double day = *(time_utc+2);
  float hour = *(time_utc+3);
  float min = *(time_utc+4);

  /* g_print ("Minute: %f, Hour: %f, Day: %f, Month: %f, Year: %f\n",min,hour,day,month,year); */
  /* calculate a float type day number which is direct proportional to the time went by */
  day += (hour + min/60.)/24.;
  /*g_print("Day JD: %f\n", day);*/

  if (month<=2){
    year --;
    month += 12;
  }

  *time_jd = (int)(365.25*(year+4716.)) + (int)(30.6001*(month+1.)) - (int)(year/100.) + (int)(year/400.) + day + 2. - 1524.5;
  /* g_print ("Julian Day: %f, Day JD: %f, hour JD %f\n", *time_jd, day, hour); */

  /*Explanation:
   * There are 5 main parts in this calculation.
   * 1. The first question is if our month is in January or February.
   *    This is necessary because if we have a leap year we should first add the extra day in March not allready in January or February.
   *    We than simply see January and February as the 13th or 14th month of the last year.
   *
   * 2. The secont main part is the expression "(int)(365.25*(year+4716))" and partly combined with the " - 1524.5".
   *    This mainly results from the definition of the Julian Day which is defined to be zero at the 1. January of -4712 at 12 UT
   *    Due to step 1 this would lead to a negativ year. To avoid this we define the 0 point 4 years earlier and substrate these 1461 days.
   *    To get then the right amount of days on the 1st of March 4712 at 00:00 we have to add these "lost" days which ar 59.5.
   *    It is 59.5 and not 58.5 because the year -4712 was also a leap year.
   *    The (int) definition is necessary because we only want to add a day every 4 years. The int gets rid of the .25 .50 or .75.
   *
   * 3. (int)(30.6001*(month+1)): This part is not as obvios as the rest of the calculation. This formular is an "estimation" function of the days passed until the given month.
   *    Due to step 1 we don't need to look at January and February so let's start with March. How many days have passed in March since the "beginning" of the year.
   *    Where the beginning is March it selfe. It is 0. How many days have passed on 1st of April? 31. The following months are
   *    61, 92, 122, 153, 184, 214, 245, 275, 306 and 337 on 1. February. The first of March is then again 0.
   *    We can receive these values also whith the equation (int)(30.6*(month+1)-122). Due to "errors" in floating point arithmetic 30.6 can not be represented in binary
   *    and can lead for example in april to 152.999998 which is in (int) 152 and not 153. Therefore we do not multiply by 30.6 and use 30.6001 instead.
   *
   * 4. The adding of the days (+ day) should be obvious. BUT we have to calculate (+ day - 1) because if we have the 3rd of Oktober for example in principle it is
   *    the 2nd of Oktober + a few hous.
   *
   * 5. 2 - (int)(year/100) + (int)(year/400): In the gregorian calendar we have every 4 years a leap year (step 2) but every 100th year we do NOT have a leap year
   *    but every 400th year it is still a leap year. This is to compensate the rotation around the Sun more precise. This is calculated by "- (int)(year/100) + (int)(year/400)"
   *    The 2 is an offset because at JD 0 the Gregorian calendar did not exist at this time. It was defined that the 5. Oktober 1582 is the 15. Oktober 1582. If we calculate
   *    the JD of this date with the Julian calendar we see a offset of two days to the same JD calculated with the Gregorian calendar. This is where the 2 comes from.
   *    Just keep in mind that the Julian Day is not the same as the Julian Date!
   *
   * --> combining all of this leads to the above formula.
   */
}


/* Calculate the Sidereal time in degree (German: Sternzeit) */
void calc_sidereal_time(float longitude, double *time_jd, double *sidereal_time) {
  float hours_passed;
  double jd_sidereal;

  /* Julian date at 0 UT at the current date */
  if (*time_jd - (int)*time_jd >= 0.5){
    jd_sidereal = (int)(*time_jd) + 0.5;
  }
  else{
    jd_sidereal = (int)(*time_jd) - 0.5;
  }

  hours_passed = 24*(*time_jd-jd_sidereal);
  /*g_print("hours_passed %f\n", hours_passed);*/
  double T = (jd_sidereal - 2451545)/36525;

  /* This formular calculatest the amount of 100 years since the fundamental epoche J2000 (1.1.2000)
   * This formular is NOT allowed to get julian days which end on .5 (00:00), because the formular below only considers the rotation around the sun NOT its' own prcession.
   * So we only want to know in which "direktion" we look every day at 00:00 in Greenwich in comparison to the "fixed" stars. */

  *sidereal_time = 100.46061837 + T*(36000.770053608 + T*(0.000387933 - T/38710000.0)) + longitude + hours_passed*1.00273790935*15;
  *sidereal_time = remainder(*sidereal_time, 360);

  /* g_print("time_jd: %f, sidereal_time %f, hours_passed: %f\n", *time_jd, *sidereal_time, hours_passed); */


  /* This formula for the sidereal time is not that easy to explain and strongly based on numerical calculations which do not allow analytical arguments.
   * For further details I recommend reading the reports of the International Atronomical Union (IAU) and the publications from P. K. Seidelmann from 1981.
   *
   * Short: The first part: "100.46061837 + T*(36000.770053608 + T*(0.000387933 - T/38710000))" calculates the startime at 00:00 in Greenwich.
   *        The first offset at 100.46061837 degree is due to its definition to the 1.1.2000 which was the offset at this spesific date (fundamental epoch J2000).
   *
   *        The first linear part T*36000.770053608 is also quite easily to anderstand. At 00:00 every day the earth had to spin about 361 degree to be at the same
   *        position towards the sun due to the rotation around the sun. This means every day at 00:00 we look about 1 degree differently to a fixed star then the day
   *        bevore at the same time. In one year this is 360 degree and in 100 years this is 36000 degree.
   *        The quadratic and cubic parts are due to different effects like the precission and nutation of the ecliptic.
   *        --> more detailed explanations where these terms come from may follow. If someone knows the answer PLEASE tell me!
   *
   *        The adding of the longitude is to get the correct sidereal time at the location of interest. This is necessary because the sidereal time
   *        is not the same around the globe.
   *
   *        The adding of the time multiplied by "1.00273790935" is necessary because a classical day is more than a 360 degree rotation of the earth.
   *        This can be explained because we can see the same star at the skye one time more often than the sun due to the earth rotating around the sun.
   *        The precise time can not be took into account in the calculation of T as mentioned above.
   *        It has to be multiplied by 15 to receive an angle (360/24=15).*/
}

void calc_coordinates_sun(double time_jd, float *coordinates_sun){
  float right_ascension;
  float declination;
  double T;
  double mean_longitude;
  double mean_anomaly;
  double equation_of_center;
  double ecliptic_longitude;
  double ecliptic;

  T = (time_jd-2451545.0)/36525;
  mean_longitude = 280.46645 + (36000.76983 + 0.0003032*T)*T;
  mean_anomaly = 357.52910 + (35999.05030 - (0.0001559 - 0.00000048*T)*T)*T;
  equation_of_center = (1.914600 - (0.004817 - 0.000014*T)*T)*sin(calc_deg_rad (mean_anomaly))
                     + (0.019993 - 0.000101*T)*sin(2*calc_deg_rad (mean_anomaly))
                     + (0.00029)*sin(3*calc_deg_rad (mean_anomaly));
  ecliptic_longitude = mean_longitude + equation_of_center;
  ecliptic = 23.43928 + 0.01301*T;

  ecliptic_longitude = calc_deg_rad (ecliptic_longitude);
  ecliptic = calc_deg_rad (ecliptic);

  right_ascension = atan2 (cos(ecliptic)*sin(ecliptic_longitude), cos(ecliptic_longitude));
  declination = asin (sin(ecliptic)*sin(ecliptic_longitude));

  *coordinates_sun = calc_rad_deg (right_ascension);
  *(coordinates_sun+1) = calc_rad_deg (declination);
}

void calc_coordinates_moon(double time_jd, float *coordinates_moon){
  float right_ascension;
  float declination;
  double T;
  double mean_anomaly;
  double mean_longitude; /* ecliptic longitude? */
  double mean_distance;
  double longitude_moon;
  double latitude_moon;
  double ecliptic;

  T = (time_jd-2451545.0)/36525;

  mean_longitude = 218.316 + (13.176396*36525) * T;
  mean_anomaly = 134.963 + (13.064993*36525) * T;
  mean_distance = 93.272 + (13.229350*36525) * T;

  longitude_moon = mean_longitude + 6.289 * sin(calc_deg_rad (mean_anomaly));
  latitude_moon = 5.128 * sin(calc_deg_rad(mean_distance));
  /*double distance_moon;
  distance_moon = 385001. - 20905. * cos(calc_deg_rad(mean_anomaly)); Useful for super moon calculation?*/

  ecliptic = calc_deg_rad(23.43928 + 0.01301*T);
  longitude_moon = calc_deg_rad (longitude_moon);
  latitude_moon = calc_deg_rad (latitude_moon);

  right_ascension = atan2(cos(ecliptic)*sin(longitude_moon)-sin(ecliptic)*tan(latitude_moon), cos(longitude_moon));
  declination = asin(cos(ecliptic)*sin(latitude_moon)+sin(ecliptic)*cos(latitude_moon)*sin(longitude_moon));

  /* g_print("time_jd: %f, ecliptic: %f, longitude_moon: %f, latitude_moon: %f, right_ascension %f\n", time_jd, ecliptic, longitude_moon, latitude_moon, right_ascension); */

  *coordinates_moon = calc_rad_deg (right_ascension);
  *(coordinates_moon+1) = calc_rad_deg (declination);
}

/* Convert between the rotation coordinate system and the horizontal coordinate system */
void calc_convert_rotation_horizontal(float right_ascension, float declination, float latitude, double time_sidereal, float *coordinates)
{
  latitude = calc_deg_rad (latitude);
  right_ascension = calc_deg_rad (right_ascension);
  declination = calc_deg_rad(declination);
  time_sidereal = calc_deg_rad (time_sidereal);

  float azimuth, elevation;
  float x = -cos(latitude)*sin(declination)+sin(latitude)*cos(declination)*cos(time_sidereal-right_ascension);
  float y = cos(declination)*sin(time_sidereal-right_ascension);

  if (x < 0 && y <= 0){
    azimuth = atan(y/x);
  }
  else if (x < 0 && y > 0){
    azimuth = atan(y/x)+2*M_PI;
  }
  else if (x > 0){
    azimuth = atan(y/x)+M_PI;
  }
  else if (x == 0 && y < 0){
    azimuth = M_PI/2;
  }
  else if (x == 0 && y > 0){
    azimuth = -3*M_PI/2;
  }
  else{
    azimuth = 0;
    g_print("Error calculating azimuth!");
  }

  elevation = asin( sin(latitude)*sin(declination) + cos(latitude)*cos(declination)*cos(time_sidereal-right_ascension) );

  /* TODO
   * explanation missing!!
   * */
  /* g_print("Latitude %f, RA %f, dec %f, time_sidereal %f, azimuth %f, elevation %f\n", latitude, right_ascension, declination, time_sidereal, azimuth, elevation); */

  azimuth = calc_rad_deg (azimuth);
  elevation = calc_rad_deg (elevation);

  *coordinates = azimuth;
  *(coordinates+1) = elevation;

  /*g_print("Azimuth: %f, Elevation: %f\n", azimuth, elevation);*/
}