kopia lustrzana https://gitlab.com/Zwarf/picplanner
335 wiersze
13 KiB
C
335 wiersze
13 KiB
C
/*
|
|
* calculations_transformations.c
|
|
* Copyright (C) 2021 Zwarf <zwarf@mail.de>
|
|
*
|
|
* 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 3 of the License, or
|
|
* (at your option) 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. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
|
|
#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.
|
|
*/
|
|
|
|
g_date_time_unref (date_time_utc);
|
|
|
|
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 refraction;
|
|
double azimuth, elevation;
|
|
double right_ascension, declination;
|
|
double *coordinates_horizontal = 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);
|
|
|
|
azimuth = atan2 (-y, -x);
|
|
|
|
elevation = asin (sin (latitude) * sin (declination)
|
|
+ cos (latitude) * cos (declination) * cos (time_sidereal - right_ascension));
|
|
|
|
|
|
/* TODO
|
|
* explanation missing!!
|
|
* */
|
|
|
|
azimuth = calc_rad_to_deg (azimuth);
|
|
elevation = calc_rad_to_deg (elevation);
|
|
|
|
|
|
/*
|
|
* Add atmospheric refraction
|
|
* Formula by Saemundsson (University of Iceland)
|
|
* Refraction in angle minutes.
|
|
* The formular is only an approximation. Higher precisions can only be
|
|
* achieved with the additional knowledge of the temperature and air pressure.
|
|
* However, in real situtations higher precisions are useless because one has
|
|
* nearly never a perfect horizon (e.g. mountains, buildings, trees ...).
|
|
* For photography reasons the achieved precision should be enough.
|
|
*/
|
|
|
|
refraction = 1.02 / tan ( calc_deg_to_rad (0 + 10.3 / (0 + 5.11) ) );
|
|
elevation += refraction / 60.;
|
|
|
|
|
|
coordinates_horizontal [0] = azimuth;
|
|
coordinates_horizontal [1] = elevation;
|
|
|
|
return coordinates_horizontal;
|
|
}
|
|
|
|
|
|
int
|
|
*picplanner_get_index_rise_upper_set_lower (double *coordinates_array)
|
|
{
|
|
int *rise_upper_set_lower_array = malloc (sizeof (int) * 4);
|
|
|
|
double elevation;
|
|
double upper_max = coordinates_array[0+1];
|
|
double lower_min = coordinates_array[0+1];
|
|
double last_elevation = coordinates_array[0+1];
|
|
|
|
for (int i=0; i<4; i++)
|
|
rise_upper_set_lower_array[i] = -1;
|
|
|
|
for (int i=0; i<NUM_DATA_POINTS; i++)
|
|
{
|
|
elevation = coordinates_array[i*2+1];
|
|
|
|
if (elevation >= upper_max)
|
|
{
|
|
rise_upper_set_lower_array[1] = i;
|
|
upper_max = elevation;
|
|
}
|
|
|
|
if (elevation <= lower_min)
|
|
{
|
|
rise_upper_set_lower_array[3] = i;
|
|
lower_min = elevation;
|
|
}
|
|
|
|
|
|
if (last_elevation < 0 && elevation > 0)
|
|
{
|
|
rise_upper_set_lower_array[0] = i;
|
|
}
|
|
else if (last_elevation > 0 && elevation < 0)
|
|
{
|
|
rise_upper_set_lower_array[2] = i;
|
|
}
|
|
|
|
last_elevation = elevation;
|
|
}
|
|
|
|
return rise_upper_set_lower_array;
|
|
}
|
|
|
|
|
|
int
|
|
*picplanner_get_index_dark_blue_golden (double *coordinates_array)
|
|
{
|
|
int *dark_blue_golden_array = malloc (sizeof (int) * 8);
|
|
|
|
double elevation;
|
|
double last_elevation = coordinates_array[0+1];
|
|
|
|
for (int i=0; i<8; i++)
|
|
dark_blue_golden_array[i] = -1;
|
|
|
|
for (int i=0; i<NUM_DATA_POINTS; i++)
|
|
{
|
|
elevation = coordinates_array[i*2+1];
|
|
|
|
/* Morning */
|
|
/* End of the dark night */
|
|
if (last_elevation < -18 && elevation > -18)
|
|
dark_blue_golden_array[0] = i;
|
|
/* Beginning of the blue hour */
|
|
else if (last_elevation < -8 && elevation > -8)
|
|
dark_blue_golden_array[1] = i;
|
|
/* Ending of the blue hour and beginning of the golden hour */
|
|
else if (last_elevation < -4 && elevation > -4)
|
|
dark_blue_golden_array[2] = i;
|
|
/* Ennding of the golden hour */
|
|
else if (last_elevation < 6 && elevation > 6)
|
|
dark_blue_golden_array[3] = i;
|
|
|
|
/* Evening */
|
|
/* Beginning of the golden hour */
|
|
else if (last_elevation > 6 && elevation < 6)
|
|
dark_blue_golden_array[4] = i;
|
|
/* Ending of the golden hour and beginning of the blue hour */
|
|
else if (last_elevation > -4 && elevation < -4)
|
|
dark_blue_golden_array[5] = i;
|
|
/* Ending of the blue hour */
|
|
else if (last_elevation > -8 && elevation < -8)
|
|
dark_blue_golden_array[6] = i;
|
|
/* Beginning of the dark night */
|
|
else if (last_elevation > -18 && elevation < -18)
|
|
dark_blue_golden_array[7] = i;
|
|
|
|
last_elevation = elevation;
|
|
}
|
|
|
|
return dark_blue_golden_array;
|
|
}
|
|
|
|
void
|
|
picplanner_get_date_time_from_index (GDateTime **date_time_index,
|
|
GDateTime **date_time_noon,
|
|
int index_time)
|
|
{
|
|
*date_time_index = g_date_time_add_minutes (*date_time_noon, index_time*24*60/NUM_DATA_POINTS-12*60);
|
|
}
|
|
|