kopia lustrzana https://gitlab.com/Zwarf/picplanner
225 wiersze
9.4 KiB
C
225 wiersze
9.4 KiB
C
#include "calculations_transformations.h"
|
|
|
|
|
|
/* Convert degree to radiant */
|
|
double
|
|
calc_deg_to_rad (double deg)
|
|
{
|
|
return M_PI/180*deg;
|
|
}
|
|
|
|
|
|
/* Convert radiant to degree */
|
|
double
|
|
calc_rad_to_deg (double rad)
|
|
{
|
|
return 180/M_PI*rad;
|
|
}
|
|
|
|
|
|
/* Calculate the Julian Day in UTC*/
|
|
double
|
|
calc_jd (GDateTime *date_time)
|
|
{
|
|
GDateTime *date_time_utc;
|
|
|
|
double time_jd;
|
|
double year, month, day, hour, min;
|
|
|
|
date_time_utc = g_date_time_to_utc (date_time);
|
|
|
|
year = g_date_time_get_year (date_time_utc);
|
|
month = g_date_time_get_month (date_time_utc);
|
|
day = g_date_time_get_day_of_month (date_time_utc);
|
|
hour = g_date_time_get_hour (date_time_utc);
|
|
min = g_date_time_get_minute (date_time_utc);
|
|
|
|
/*
|
|
* Calculate a float type day number which is direct proportional to the time went by
|
|
*/
|
|
day += (hour + min/60.)/24.;
|
|
|
|
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;
|
|
|
|
/*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.
|
|
*/
|
|
|
|
return time_jd;
|
|
}
|
|
|
|
|
|
|
|
/* Calculate the Sidereal time in degree (German: Sternzeit) */
|
|
double
|
|
time_jd_to_sidereal_time(double longitude,
|
|
GDateTime *date_time)
|
|
{
|
|
double time_jd;
|
|
double hours_passed;
|
|
double T;
|
|
double jd_sidereal;
|
|
double sidereal_time;
|
|
|
|
time_jd = calc_jd (date_time);
|
|
|
|
/* 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);*/
|
|
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).*/
|
|
|
|
return sidereal_time;
|
|
}
|
|
|
|
|
|
|
|
/* Convert between the rotation coordinate system and the horizontal coordinate system */
|
|
double
|
|
*picplanner_transform_rotational_to_horizontal (double *coordinates_rot,
|
|
double latitude,
|
|
double time_sidereal)
|
|
{
|
|
double x, y;
|
|
double azimuth, elevation;
|
|
double right_ascension, declination;
|
|
double *coordinates_hor = malloc (sizeof (double) * 2);
|
|
|
|
latitude = calc_deg_to_rad (latitude);
|
|
time_sidereal = calc_deg_to_rad (time_sidereal);
|
|
right_ascension = calc_deg_to_rad (coordinates_rot[0]);
|
|
declination = calc_deg_to_rad (coordinates_rot[1]);
|
|
|
|
x = -cos (latitude) * sin (declination)
|
|
+ sin (latitude) * cos (declination) * cos (time_sidereal - right_ascension);
|
|
|
|
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;
|
|
|
|
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_to_deg (azimuth);
|
|
elevation = calc_rad_to_deg (elevation);
|
|
|
|
coordinates_hor [0] = azimuth;
|
|
coordinates_hor [1] = elevation;
|
|
|
|
/*g_print("Azimuth: %f, Elevation: %f\n", azimuth, elevation);*/
|
|
|
|
return coordinates_hor;
|
|
}
|
|
|
|
|
|
int
|
|
*picplanner_get_rise_upper_set_index (double *coordinates_array)
|
|
{
|
|
int rise = -1;
|
|
int upper = -1;
|
|
int set = -1;
|
|
int *rise_upper_set_array = malloc (sizeof (int) * 2);
|
|
|
|
for (int i=0; i<NUM_DATA_POINTS; i++)
|
|
{
|
|
|
|
}
|
|
|
|
rise_upper_set_array[0] = rise;
|
|
rise_upper_set_array[1] = upper;
|
|
rise_upper_set_array[1] = set;
|
|
|
|
return rise_upper_set_array;
|
|
}
|