Ik ben op zoek naar een formule waarmee ik de Positie (compass richting vanaf het Noorden in graden) en de Hoogte (in graden) van de Maan kan berekenen op een bepaalde dag en tijd.
Input: Tijd, Datum, latitude en longitude
Na dagen zoeken en proberen kom ik tot de conclusie dat dit te hoog gegrepen is voor mij, wie kan mij hieraan helpen ?
Maan Positie
Re: Maan Positie
Zie bijvoorbeeld hoofdstuk 17 uit https://ytliu0.github.io/starCharts/doc ... charts.pdf
Re: Maan Positie
Ik begrijp niet wat hier staat (heb er een plaatje van gemaakt maar weet niet hoe ik dat hier kan plaatsen)
Op Blz. 39, bij (115) het teken met daarboven de 3 en eronder 1=0
Ik ben zo'n formules gewend:
d := Math.degtorad(0.396372 - 22.91327 * cos(g) + 4.02543 * sin(g)
- 0.387205 * cos(2 * g) + 0.051967 * sin(2 * g)
- 0.154527 * cos(3 * g) + 0.084798 * sin(3 * g));
Op Blz. 39, bij (115) het teken met daarboven de 3 en eronder 1=0
Ik ben zo'n formules gewend:
d := Math.degtorad(0.396372 - 22.91327 * cos(g) + 4.02543 * sin(g)
- 0.387205 * cos(2 * g) + 0.051967 * sin(2 * g)
- 0.154527 * cos(3 * g) + 0.084798 * sin(3 * g));
Re: Maan Positie
De griekse sigma wordt gebruikt voor sommatie, zie bv. https://nl.wikipedia.org/wiki/Sommatie.
Voor formule (115):
\(V(T)=W(T)+\displaystyle \sum_{i=0}^3 T^i\sum_jA_{ij}^{(V)}\sin \left[\phi_{ij}^{(V)}(T)\right]\)
betekent dit dat index i loopt van 0 t/m 3, wat een sommatie van 4 termen oplevert:
\(V(T)=W(T)+T^0\displaystyle\sum_jA_{0j}^{(V)}\sin \left[\phi_{0j}^{(V)}(T)\right]+T^1 \sum_jA_{1j}^{(V)}\sin \left[\phi_{1j}^{(V)}(T)\right]\)
\(+T^2 \displaystyle\sum_jA_{2j}^{(V)}\sin \left[\phi_{2j}^{(V)}(T)\right]+T^3\sum_jA_{3j}^{(V)}\sin \left[\phi_{3j}^{(V)}(T)\right]\)
vervolgens staat in elk van die 4 termen een sommatie over j:
\(V(T) = W(T) + T^0 \left(A_{0,0}^{(V)}\sin \left[\phi_{0,0}^{(V)}(T)\right]+A_{0,1}^{(V)}\sin \left[\phi_{0,1}^{(V)}(T)\right]+A_{0,2}^{(V)}\sin \left[\phi_{0,2}^{(V)}(T)\right]+... \right)\)
\(+T^1 \left(A_{1,0}^{(V)}\sin \left[\phi_{1,0}^{(V)}(T)\right]+A_{1,1}^{(V)}\sin \left[\phi_{1,1}^{(V)}(T)\right]+A_{1,2}^{(V)}\sin \left[\phi_{1,2}^{(V)}(T)\right]+... \right)\)
\(+T^2 \left(A_{2,0}^{(V)}\sin \left[\phi_{2,0}^{(V)}(T)\right]+A_{2,1}^{(V)}\sin \left[\phi_{2,1}^{(V)}(T)\right]+A_{2,2}^{(V)}\sin \left[\phi_{2,2}^{(V)}(T)\right]+... \right)\)
\(+T^3 \left(A_{3,0}^{(V)}\sin \left[\phi_{3,0}^{(V)}(T)\right]+A_{3,1}^{(V)}\sin \left[\phi_{3,1}^{(V)}(T)\right]+A_{3,2}^{(V)}\sin \left[\phi_{3,2}^{(V)}(T)\right]+... \right)\)
In dit geval zijn er voor j geen gefixeerde waarden, maar gebruiken ze voor elk van deze 4 sommaties net zo veel termen totdat \(A_{ij}\) onder een drempelwaarde zakt = totdat de termen zodanig klein worden dat je ze kan verwaarlozen.
Onder aan pag. 39 staat een GitHub link, waar je de waarden voor alle constanten vindt
(zie bv. https://github.com/ytliu0/ElpMpp02/blob ... pMpp02.pdf
of https://github.com/ytliu0/ElpMpp02/blob ... ElpMpp02.h)
Voor formule (115):
\(V(T)=W(T)+\displaystyle \sum_{i=0}^3 T^i\sum_jA_{ij}^{(V)}\sin \left[\phi_{ij}^{(V)}(T)\right]\)
betekent dit dat index i loopt van 0 t/m 3, wat een sommatie van 4 termen oplevert:
\(V(T)=W(T)+T^0\displaystyle\sum_jA_{0j}^{(V)}\sin \left[\phi_{0j}^{(V)}(T)\right]+T^1 \sum_jA_{1j}^{(V)}\sin \left[\phi_{1j}^{(V)}(T)\right]\)
\(+T^2 \displaystyle\sum_jA_{2j}^{(V)}\sin \left[\phi_{2j}^{(V)}(T)\right]+T^3\sum_jA_{3j}^{(V)}\sin \left[\phi_{3j}^{(V)}(T)\right]\)
vervolgens staat in elk van die 4 termen een sommatie over j:
\(V(T) = W(T) + T^0 \left(A_{0,0}^{(V)}\sin \left[\phi_{0,0}^{(V)}(T)\right]+A_{0,1}^{(V)}\sin \left[\phi_{0,1}^{(V)}(T)\right]+A_{0,2}^{(V)}\sin \left[\phi_{0,2}^{(V)}(T)\right]+... \right)\)
\(+T^1 \left(A_{1,0}^{(V)}\sin \left[\phi_{1,0}^{(V)}(T)\right]+A_{1,1}^{(V)}\sin \left[\phi_{1,1}^{(V)}(T)\right]+A_{1,2}^{(V)}\sin \left[\phi_{1,2}^{(V)}(T)\right]+... \right)\)
\(+T^2 \left(A_{2,0}^{(V)}\sin \left[\phi_{2,0}^{(V)}(T)\right]+A_{2,1}^{(V)}\sin \left[\phi_{2,1}^{(V)}(T)\right]+A_{2,2}^{(V)}\sin \left[\phi_{2,2}^{(V)}(T)\right]+... \right)\)
\(+T^3 \left(A_{3,0}^{(V)}\sin \left[\phi_{3,0}^{(V)}(T)\right]+A_{3,1}^{(V)}\sin \left[\phi_{3,1}^{(V)}(T)\right]+A_{3,2}^{(V)}\sin \left[\phi_{3,2}^{(V)}(T)\right]+... \right)\)
In dit geval zijn er voor j geen gefixeerde waarden, maar gebruiken ze voor elk van deze 4 sommaties net zo veel termen totdat \(A_{ij}\) onder een drempelwaarde zakt = totdat de termen zodanig klein worden dat je ze kan verwaarlozen.
Onder aan pag. 39 staat een GitHub link, waar je de waarden voor alle constanten vindt
(zie bv. https://github.com/ytliu0/ElpMpp02/blob ... pMpp02.pdf
of https://github.com/ytliu0/ElpMpp02/blob ... ElpMpp02.h)
Re: Maan Positie
Mocht er iemand de kennis, tijd en zin hebben om de formule(s) om te zetten in een pascal achtige code, dan zou ik daar erg blij mee zijn.
Re: Maan Positie
Hieronder een C-versie van de ELP/MPP02L benadering beschreven onderaan pag. 39 van
https://ytliu0.github.io/starCharts/doc ... charts.pdf
(nauwkeurigheid hoeken gemiddeld binnen 1 boogminuut en en afstand aarde-maan gemiddeld binnen 100 km).
Het aantal constanten is daardoor zeer beperkt, ik heb ze allemaal opgenomen in de code.
C-code:
Output:
Dezelfde tabel met hoge nauwkeurigheid vind je onderaan op pag. 6 van
https://github.com/ytliu0/ElpMpp02/blob ... pMpp02.pdf
Hier nog wat toelichting bij de taal C:
Resteert nog de omzetting van de geocentrische coordinaten (x, y, z) naar azimuth en altitude,
zie daarvoor bv. paragraaf 17.2 (pag. 40) van https://ytliu0.github.io/starCharts/doc ... charts.pdf
Kom je hiermee verder?
https://ytliu0.github.io/starCharts/doc ... charts.pdf
(nauwkeurigheid hoeken gemiddeld binnen 1 boogminuut en en afstand aarde-maan gemiddeld binnen 100 km).
Het aantal constanten is daardoor zeer beperkt, ik heb ze allemaal opgenomen in de code.
C-code:
Code: Selecteer alles
//---------------------------------------------------------------------------
// Lunar Ephemeris
//---------------------------------------------------------------------------
// Source:
// https://ytliu0.github.io/starCharts/docs/star_charts.pdf
// chapter 17
// Based on truncated ELP/MPP02 series using the following parameters:
// corr = 1 = Parameters fitted to JPL ephemeris DE405/DE406
// AthU = 1.454441e-04 radians = 30.000000 arcsecs
// AthV = 1.454441e-04 radians = 30.000000 arcsecs
// AthR = 100.000000 km
// tau = 50.000000
// Total number of terms in the original series: 35901
// Total number of terms in the truncated series: 42
// Statistical evaluation:
// Estimated error of the truncated series by comparing it to the original series
// at 10000 randomly chosen T between -50.000000 and 10.000000:
// Maximum deviation found in ecliptic longitude: 213.496798 arcsecs
// Root mean square deviation in ecliptic longitude: 47.665830 arcsecs
// Maximum deviation found in ecliptic latitude: 158.214315 arcsecs
// Root mean square deviation in ecliptic latitude: 34.778046 arcsecs
// Maximum deviation found in geocentric distance: 336.995129 km
// Root mean square deviation in geocentric distance: 85.856277 km
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
// standard C libs:
#include <stdio.h>
#include <math.h>
//---------------------------------------------------------------------------
// constants:
const double PI = 3.14159265358979323846;
const double TPI = 2.0*PI;
//---------------------------------------------------------------------------
// Arguments for the ELP/MPP02 series:
double W1, D, F, L, Lp;
// coefficients for the ELP/MPP02 main series:
int n_main_long=19;
int i_main_long[19][4]={
{0,2,0,0},{0,-2,1,0},{0,0,1,0},{0,2,1,0},{0,0,2,0},{0,0,3,0},{0,0,-1,1},{0,0,0,1},
{0,0,1,1},{1,0,0,0},{2,0,-1,-1},{2,0,0,-1},{2,0,-2,0},{2,0,-1,0},{2,-2,0,0},{2,0,0,0},
{2,0,1,0},{4,0,-2,0},{4,0,-1,0}};
double A_main_long[19]={
-1.99547249801829270800e-03, 1.91662856547817478400e-04, 1.09759808626191612800e-01, -2.18649080896658580300e-04,
3.72834182278018950000e-03, 1.75133187917569992000e-04, -7.14234214377495904400e-04, -3.23088338567132336800e-03,
-5.30290959258585535000e-04, -6.05959497615181906000e-04, 9.95981590587821435600e-04, 7.98626872936074750500e-04,
1.02613505456892778600e-03, 2.22356802322449947900e-02, 2.67505934260747291300e-04, 1.14896669562630695400e-02,
9.30629894589579480800e-04, 1.49189651057941429400e-04, 1.86313093175213899100e-04};
int n_main_lat=10;
int i_main_lat[10][4]={
{0,1,0,0},{0,-1,1,0},{0,1,1,0},{0,-1,2,0},{0,1,2,0},{2,-1,-1,0},{2,1,-1,0},{2,-1,0,0},{2,1,0,0},{2,-1,1,0}};
double A_main_lat[10]={
8.95026190647627839500e-02, 4.84665164815963230000e-03, 4.89742857492255542000e-03, 1.53975229251259538900e-04,
3.00157605100467080500e-04, 8.07574106044906285200e-04, 9.67124565462755482600e-04, 3.02355257133907857600e-03,
5.68495899772767385500e-04, 1.61720283648969048700e-04};
int n_main_dist=12;
int i_main_dist[12][4]={
{0,0,0,0},{0,0,1,0},{0,0,2,0},{0,0,-1,1},{0,0,1,1},{1,0,0,0},
{2,0,-1,-1},{2,0,0,-1},{2,0,-2,0},{2,0,-1,0},{2,0,0,0},{2,0,1,0}};
double A_main_dist[12]={
3.85000529032284044600e+05, -2.09053549352050904400e+04, -5.69925115335094687900e+02, -1.29620222172050574700e+02,
1.04755292674270279460e+02, 1.08742701272463051500e+02, -1.52137809425890736800e+02, -2.04586116460873086000e+02,
2.46158474853557947900e+02, -3.69911089639807642000e+03, -2.95596755797241485200e+03, -1.70733077124770574100e+02};
double main_long, main_lat, main_dist;
double pert_longT1;
//---------------------------------------------------------------------------
// restrict x to [-pi,pi>
double mod2pi(double x)
{
return(x - TPI*floor((x + PI)/TPI));
}
// Compute the lunar and planetary arguments used in the ELP/MPP02 series:
// set Elp arguments: W1, D, F, L, Lp
void compute_Elp_arguments(double T)
{
double T2 = T*T;
double T3 = T*T2;
double T4 = T2*T2;
W1 = -2.4728412163487064 + mod2pi(8399.6847300719273*T) + mod2pi(-3.3189520425500949e-05*T2)
+ mod2pi(3.1102494491060616e-08*T3) + mod2pi(-2.0328237648922845e-10*T4);
double W2 = 1.4547895404440774 + mod2pi(70.993305448479248*T) + mod2pi(-0.00018548192818782712*T2)
+ mod2pi(-2.1961637966359414e-07*T3) + mod2pi(1.0327016221314225e-09*T4);
double W3 = 2.1824388474237688 + mod2pi(-33.781427419672326*T) + mod2pi(3.0816644950982666e-05*T2)
+ mod2pi(3.6447710769397583e-08*T3) + mod2pi(-1.738541860458796e-10*T4);
double Ea = 1.7534699452640699 + mod2pi(628.30758508103168*T) + mod2pi(-9.7932363584126268e-08*T2)
+ mod2pi(4.3633231299858238e-11*T3) + mod2pi(7.2722052166430391e-13*T4);
double pomp = 1.7965956331206001 + mod2pi(0.0056298669711442699*T) + mod2pi(2.5659491293243858e-06*T2)
+ mod2pi(-5.7275888286280579e-10*T3) + mod2pi(5.5166948773454105e-11*T4);
// Mean longitude of the Moon:
W1 = mod2pi(W1);
// Arguments of Delaunay:
D = mod2pi(W1 - Ea + PI);
F = mod2pi(W1 - W3);
L = mod2pi(W1 - W2);
Lp = mod2pi(Ea - pomp);
}
//---------------------------------------------------------------------------
// Sum the ELP/MPP02 main and perturbation series:
void compute_Elp_sums(void)
{
int i;
double phase;
// main sums (long and lat: sin-sums; dist: cos-sum):
main_long=0.0;
for(i=0;i<n_main_long;i++){
phase=i_main_long[i][0]*D + i_main_long[i][1]*F + i_main_long[i][2]*L + i_main_long[i][3]*Lp;
main_long += A_main_long[i]*sin(phase);
}
main_lat=0.0;
for(i=0;i<n_main_lat;i++){
phase=i_main_lat[i][0]*D + i_main_lat[i][1]*F + i_main_lat[i][2]*L + i_main_lat[i][3]*Lp;
main_lat += A_main_lat[i]*sin(phase);
}
main_dist=0.0;
for(i=0;i<n_main_dist;i++){
phase=i_main_dist[i][0]*D + i_main_dist[i][1]*F + i_main_dist[i][2]*L + i_main_dist[i][3]*Lp;
main_dist += A_main_dist[i]*cos(phase);
}
// perturbations sums: pert_longT1 is the only non-zero for our approximation:
pert_longT1 = 8.1293558048447e-06 * sin(Lp);
}
//---------------------------------------------------------------------------
// Calculate the Moon's geocentric X,Y,Z coordinates with respect to
// J2000.0 mean ecliptic and equinox.
void getX2000(double JD, double &X, double &Y, double &Z)
{
// T is the TDB Julian century from J2000.0 = (JD - 2451545)/36525
double T=(JD - 2451545.0)/36525.0;
double T2 = T*T;
double T3 = T*T2;
double T4 = T2*T2;
double T5 = T2*T3;
compute_Elp_arguments(T);
compute_Elp_sums();
// Moon's longitude, latitude and distance
double longM = W1 + main_long + mod2pi(pert_longT1*T);
double latM = main_lat;
// ra0 = a0(DE405)/a0(ELP) = 384747.961370173/384747.980674318 = 0.99999994982652041.
double r = 0.99999994982652041*main_dist;
double x0 = r*cos(longM)*cos(latM);
double y0 = r*sin(longM)*cos(latM);
double z0 = r*sin(latM);
// Precession matrix
double P = 0.10180391e-4*T + 0.47020439e-6*T2 - 0.5417367e-9*T3 - 0.2507948e-11*T4 + 0.463486e-14*T5;
double Q = -0.113469002e-3*T + 0.12372674e-6*T2 + 0.12654170e-8*T3 - 0.1371808e-11*T4 - 0.320334e-14*T5;
double sq = sqrt(1 - P*P - Q*Q);
double p11 = 1 - 2*P*P;
double p12 = 2*P*Q;
double p13 = 2*P*sq;
double p21 = 2*P*Q;
double p22 = 1-2*Q*Q;
double p23 = -2*Q*sq;
double p31 = -2*P*sq;
double p32 = 2*Q*sq;
double p33 = 1 - 2*P*P - 2*Q*Q;
// Finally, calculate the return values X, Y and Z:
// (= the components of the position vector; wrt J2000.0 mean ecliptic and equinox)
X = p11*x0 + p12*y0 + p13*z0;
Y = p21*x0 + p22*y0 + p23*z0;
Z = p31*x0 + p32*y0 + p33*z0;
}
//---------------------------------------------------------------------------
// Julian date:
//---------------------------------------------------------------------------
// truncate x = drop fractional part of x:
double trunc(double x)
{
double sign = 1.0;
if(x<0.0){
sign = -1.0;
x = -x;
}
return(sign*floor(x));
}
// calculate Julian Date for date D-M-Y time h:m:s
double JulianDate(int D, int M, int Y, double h, double m, double s)
{
double day=(double)D + (h+m/60.0+s/3600.0)/24.0;
double year, month, tmp, a, b, c;
if(M<=2){
year = Y-1;
month = M+12;
}
else{
year = Y;
month = M;
}
// 15-10-1582 = start Gregorian calendar:
if((Y<1582) || (Y==1582 && M<10) || (Y==1582 && M==10 && D<15)) // before Gregorian cal:
a = 0;
else{ // Gregorian cal:
tmp = trunc(year/100.0);
a = 2 - tmp + floor(tmp/4.0);
}
if(year<0)
b = trunc(365.25*year-0.75);
else
b = trunc(365.25*year);
c = trunc(30.6001 * (month+1));
return(a + b + c + day + 1720994.5);
}
//---------------------------------------------------------------------------
// Demo:
//---------------------------------------------------------------------------
// Demonstration output for our ELP/MPP02 series approximation
// Compare with the high precision version on page 6 in
// https://github.com/ytliu0/ElpMpp02/blob/master/docs/ElpMpp02.pdf
void showTableHeader(void)
{
printf("\n");
printf("Lunar geocentric rectangular coordinates (in km), fitted to DE405/DE406\n");
printf("-------------------------------------------------------------------------------\n");
printf("Jul. Date Date Time X Y Z\n");
printf("-------------------------------------------------------------------------------\n");
}
void showMoonPosition(int D, int M, int Y, int h, int m, double s)
{
double JD=JulianDate(D,M,Y,h,m,s); // calculate Julian Date
double x,y,z; // lunar coordinates
getX2000(JD,x,y,z); // calculate lunar coordinates on JD
printf("%.2lf %5d-%02d-%02d %02d:%02d:%02.0lf %10.0lf %10.0lf %10.0lf\n",JD,Y,M,D,h,m,s,x,y,z);
}
// main program:
int main(int argc, char* argv[])
{
showTableHeader();
showMoonPosition(13,6,2192,4,4,48);
showMoonPosition(7,12,1490,19,55,12);
showMoonPosition(16,6,789,11,45,36);
showMoonPosition(25,12,87,3,36,0);
showMoonPosition(3,7,-614,19,26,24);
return 0; // exit without errors
}
Output:
Code: Selecteer alles
Lunar geocentric rectangular coordinates (in km), fitted to DE405/DE406
-------------------------------------------------------------------------------
Jul. Date Date Time X Y Z
-------------------------------------------------------------------------------
2521835.67 2192-06-13 04:04:48 -184072 346042 30400
2265621.33 1490-12-07 19:55:12 -297887 -213732 -23223
2009406.99 789-06-16 11:45:36 350156 -201162 1638
1753192.65 87-12-25 03:36:00 90182 351827 13573
1496978.31 -614-07-03 19:26:24 -403080 -2692 -28476
https://github.com/ytliu0/ElpMpp02/blob ... pMpp02.pdf
Hier nog wat toelichting bij de taal C:
Code: Selecteer alles
// dit is commentaar
double = double precision floating number
int = integer number
int i_main_long[19][4] is een matrix (=array) van 19 bij 4, genaamd i_main_long, bestaande uit integers
floor(x) = het grootste gehele getal dat kleiner is dan x
a+=b is hetzelfde als a=a+b
if((Y<1582) || (Y==1582 && M<10))
|| = OF = or
&& = EN = and
Y==1582 test of Y gelijk is aan 1582 (return: true als dat zo is, en anders false)
Y=1582 maakt Y gelijk aan 1582
printf("%.2lf %04d %5d\n", x, a, b)
schrijf double getal x met 2 cijfers achter de komma
schrijf spatie
schrijf integer a 4 karakters breed, eventueel voorafgegaan door nullen (0023)
schrijf spatie
schrijf integer b 5 karakters breed, eventueel voorafgegaan door spaties ( 23)
\n = nieuwe regel
zie daarvoor bv. paragraaf 17.2 (pag. 40) van https://ytliu0.github.io/starCharts/doc ... charts.pdf
Kom je hiermee verder?
Re: Maan Positie
Ik heb de code omgezet naar Delphi (Pascal):
Maar daarin zit blijkbaar een fout.
Datum/Tijd: 13,6,2192,4,4,48
Geeft: 2521835,67 356672,731550951 -0,00331534775033703 -15,2150931367566
Moet zijn: 2521835.67 -184072 346042 30400
Zou je een paar tussen waardes willen geven, zodat ik kan zoeken waar het mis gaat ?
Code: Selecteer alles
const
PI = 3.14159265358979323846;
TPI = 2.0 * PI;
var
W1, D, F, L, Lp: double;
n_main_long: Integer = 19;
i_main_long: array[0..18, 0..3] of Integer = (
(0, 2, 0, 0), (0, -2, 1, 0), (0, 0, 1, 0), (0, 2, 1, 0), (0, 0, 2, 0),
(0, 0, 3, 0), (0, 0, -1, 1), (0, 0, 0, 1), (0, 0, 1, 1), (1, 0, 0, 0),
(2, 0, -1, -1), (2, 0, 0, -1), (2, 0, -2, 0), (2, 0, -1, 0), (2, -2, 0, 0),
(2, 0, 0, 0), (2, 0, 1, 0), (4, 0, -2, 0), (4, 0, -1, 0));
A_main_long: array[0..18] of Double = (
-1.99547249801829270800e-03, 1.91662856547817478400e-04, 1.09759808626191612800e-01, -2.18649080896658580300e-04,
3.72834182278018950000e-03, 1.75133187917569992000e-04, -7.14234214377495904400e-04, -3.23088338567132336800e-03,
-5.30290959258585535000e-04, -6.05959497615181906000e-04, 9.95981590587821435600e-04, 7.98626872936074750500e-04,
1.02613505456892778600e-03, 2.22356802322449947900e-02, 2.67505934260747291300e-04, 1.14896669562630695400e-02,
9.30629894589579480800e-04, 1.49189651057941429400e-04, 1.86313093175213899100e-04);
n_main_lat: Integer = 10;
i_main_lat: array[0..9, 0..3] of Integer = (
(0,1,0,0),(0,-1,1,0),(0,1,1,0),(0,-1,2,0),(0,1,2,0),(2,-1,-1,0),(2,1,-1,0),(2,-1,0,0),(2,1,0,0),(2,-1,1,0));
A_main_lat: array[0..9] of Double = (
8.95026190647627839500e-02, 4.84665164815963230000e-03, 4.89742857492255542000e-03, 1.53975229251259538900e-04,
3.00157605100467080500e-04, 8.07574106044906285200e-04, 9.67124565462755482600e-04, 3.02355257133907857600e-03,
5.68495899772767385500e-04, 1.61720283648969048700e-04);
n_main_dist: Integer = 12;
i_main_dist: array[0..11, 0..3] of Integer = (
(0,0,0,0),(0,0,1,0),(0,0,2,0),(0,0,-1,1),(0,0,1,1),(1,0,0,0),
(2,0,-1,-1),(2,0,0,-1),(2,0,-2,0),(2,0,-1,0),(2,0,0,0),(2,0,1,0));
A_main_dist: array[0..11] of Double = (
3.85000529032284044600e+05, -2.09053549352050904400e+04, -5.69925115335094687900e+02, -1.29620222172050574700e+02,
1.04755292674270279460e+02, 1.08742701272463051500e+02, -1.52137809425890736800e+02, -2.04586116460873086000e+02,
2.46158474853557947900e+02, -3.69911089639807642000e+03, -2.95596755797241485200e+03, -1.70733077124770574100e+02);
main_long, main_lat, main_dist: Double;
pert_longT1: Double;
Geocentrische_Maan_Positie_X: Double;
Geocentrische_Maan_Positie_Y: Double;
Geocentrische_Maan_Positie_Z: Double;
//------------------------------------------------------------------------------
procedure TForm1.Button1Click(Sender: TObject);
var
JD: Double;
x,y,z: double;
begin
Memo1.Lines.Add( 'Lunar geocentric rectangular coordinates (in km), fitted to DE405/DE406' );
Memo1.Lines.Add( '-------------------------------------------------------------------------------' );
Memo1.Lines.Add( 'Jul. Date Date Time X Y Z' );
Memo1.Lines.Add( '-------------------------------------------------------------------------------' );
JD := JulianDate(13,6,2192,4,4,48);
// (13,6,2192,4,4,48);
// (7,12,1490,19,55,12);
// (16,6,789,11,45,36);
// (25,12,87,3,36,0);
// (3,7,-614,19,26,24);
getX2000(JD); // calculate lunar coordinates on JD
//Deze variable zijn nu berekend:
//Geocentrische_Maan_Positie_X
//Geocentrische_Maan_Positie_Y
//Geocentrische_Maan_Positie_Z
Memo1.Lines.Add( FloatToStr( JD ) + ' ' + FloatToStr( Geocentrische_Maan_Positie_X ) + ' ' + FloatToStr( Geocentrische_Maan_Positie_Y ) + ' ' + FloatToStr( Geocentrische_Maan_Positie_Z ) );
end;
//------------------------------------------------------------------------------
// restrict x to [-pi,pi>
function mod2pi(x: double): Double;
begin
result := (x - TPI*floor((x + PI)/TPI));
end;
//------------------------------------------------------------------------------
// Compute the lunar and planetary arguments used in the ELP/MPP02 series:
// set Elp arguments: W1, D, F, L, Lp
procedure TForm1.compute_Elp_arguments(T: double);
var
T2, T3, T4, W1, W2, W3: Double;
Ea: Double ;
pomp: Double;
D, F, L, Lp: Double;
begin
T2 := T * T;
T3 := T * T2;
T4 := T2 * T2;
W1 := -2.4728412163487064 + mod2pi(8399.6847300719273*T) + mod2pi(-3.3189520425500949e-05*T2)
+ mod2pi(3.1102494491060616e-08*T3) + mod2pi(-2.0328237648922845e-10*T4);
W2 := 1.4547895404440774 + mod2pi(70.993305448479248*T) + mod2pi(-0.00018548192818782712*T2)
+ mod2pi(-2.1961637966359414e-07*T3) + mod2pi(1.0327016221314225e-09*T4);
W3 := 2.1824388474237688 + mod2pi(-33.781427419672326*T) + mod2pi(3.0816644950982666e-05*T2)
+ mod2pi(3.6447710769397583e-08*T3) + mod2pi(-1.738541860458796e-10*T4);
Ea := 1.7534699452640699 + mod2pi(628.30758508103168*T) + mod2pi(-9.7932363584126268e-08*T2)
+ mod2pi(4.3633231299858238e-11*T3) + mod2pi(7.2722052166430391e-13*T4);
pomp := 1.7965956331206001 + mod2pi(0.0056298669711442699*T) + mod2pi(2.5659491293243858e-06*T2)
+ mod2pi(-5.7275888286280579e-10*T3) + mod2pi(5.5166948773454105e-11*T4);
// Mean longitude of the Moon:
W1 := mod2pi(W1);
// Arguments of Delaunay:
D := mod2pi(W1 - Ea + PI);
F := mod2pi(W1 - W3);
L := mod2pi(W1 - W2);
Lp := mod2pi(Ea - pomp);
end;
//---------------------------------------------------------------------------
// Sum the ELP/MPP02 main and perturbation series:
procedure TForm1.compute_Elp_sums();
var
i: Integer;
phase: Double;
main_long: Double;
main_lat: Double;
begin
// main sums (long and lat: sin-sums; dist: cos-sum):
main_long := 0.0;
for i := 0 to n_main_long - 1 do
begin
phase := i_main_long[i][0] * D + i_main_long[i][1] * F + i_main_long[i][2] * L + i_main_long[i][3] * Lp;
main_long := main_long + A_main_long[i] * Sin(phase);
end;
main_lat := 0.0;
for i := 0 to n_main_lat - 1 do
begin
phase := i_main_lat[i][0] * D + i_main_lat[i][1] * F + i_main_lat[i][2] * L + i_main_lat[i][3] * Lp;
main_lat := main_lat + A_main_lat[i] * sin(phase);
end;
main_dist := 0.0;
for i := 0 to n_main_dist - 1 do
begin
phase := i_main_dist[i][0]*D + i_main_dist[i][1]*F + i_main_dist[i][2]*L + i_main_dist[i][3]*Lp;
main_dist := main_dist + A_main_dist[i] * cos(phase);
end;
// perturbations sums: pert_longT1 is the only non-zero for our approximation:
pert_longT1 := 8.1293558048447e-06 * sin(Lp);
end;
//---------------------------------------------------------------------------
// Calculate the Moon's geocentric X,Y,Z coordinates with respect to
// J2000.0 mean ecliptic and equinox.
procedure TForm1.getX2000(JD: double);
// T is the TDB Julian century from J2000.0 = (JD - 2451545)/36525
var
T, T2, T3, T4, T5: Double;
longM, latM: Double;
r, x0, y0, z0: Double;
P, Q: Double;
sq: Double;
p11, p12, p13, p21, p22, p23, p31, p32, p33: Double;
begin
T := (JD - 2451545.0) / 36525.0;
T2 := T * T;
T3 := T * T2;
T4 := T2 * T2;
T5 := T2 * T3;
compute_Elp_arguments(T);
compute_Elp_sums();
// Moon's longitude, latitude and distance
longM := W1 + main_long + mod2pi(pert_longT1 * T);
latM := main_lat;
// ra0 = a0(DE405)/a0(ELP) = 384747.961370173/384747.980674318 = 0.99999994982652041.
r := 0.99999994982652041 * main_dist;
x0 := r * cos(longM) * cos(latM);
y0 := r * sin(longM) * cos(latM);
z0 := r * sin(latM);
// Precession matrix
P := 0.10180391e-4 * T + 0.47020439e-6 * T2 - 0.5417367e-9 * T3 - 0.2507948e-11 * T4 + 0.463486e-14 * T5;
Q := -0.113469002e-3 * T + 0.12372674e-6 * T2 + 0.12654170e-8 * T3 - 0.1371808e-11 * T4 - 0.320334e-14 * T5;
sq := sqrt(1 - P * P - Q * Q);
p11 := 1 - 2 * P * P;
p12 := 2 * P * Q;
p13 := 2 * P * sq;
p21 := 2 * P * Q;
p22 := 1 - 2 * Q * Q;
p23 := -2 * Q * sq;
p31 := -2 * P * sq;
p32 := 2 * Q * sq;
p33 := 1 - 2 * P * P - 2 * Q * Q;
// Finally, calculate the return values X, Y and Z:
// (= the components of the position vector; wrt J2000.0 mean ecliptic and equinox)
Geocentrische_Maan_Positie_X := p11 * x0 + p12 * y0 + p13 * z0;
Geocentrische_Maan_Positie_Y := p21 * x0 + p22 * y0 + p23 * z0;
Geocentrische_Maan_Positie_Z := p31 * x0 + p32 * y0 + p33 * z0;
end;
//------------------------------------------------------------------------------
// truncate x = drop fractional part of x:
function trunc(x: double): double;
var
sign: double;
begin
sign := 1.0;
if(x < 0.0) then
begin
sign := -1.0;
x := -x;
end;
result := (sign * floor(x));
end;
//------------------------------------------------------------------------------
// calculate Julian Date for date D-M-Y time h:m:s
//double JulianDate(D, M, Y: Integer; , double h, double m, double s)
function JulianDate(D, Maand, Y: Integer; h, m, s: double): Double;
var
day, year, month, tmp, a, b, c: double ;
begin
day := D + (h + m / 60.0 + s / 3600.0) / 24.0;
if( Maand <= 2) then
begin
year := Y - 1;
month := Maand + 12;
end
else
begin
year := Y;
month := Maand;
end;
// 15-10-1582 = start Gregorian calendar:
if ( (Y < 1582) or ( (Y = 1582) and ( Maand < 10 ) ) or ( ( Y = 1582) and (Maand = 10) and (D < 15) ) ) then // before Gregorian cal:
begin
a := 0;
end
else
begin // Gregorian cal:
tmp := trunc(year / 100.0);
a := 2 - tmp + floor(tmp / 4.0);
end;
if(year < 0) then
b := trunc(365.25 * year - 0.75)
else
b := trunc(365.25 * year);
c := trunc(30.6001 * (month + 1));
result := (a + b + c + day + 1720994.5);
end;
//------------------------------------------------------------------------------
Datum/Tijd: 13,6,2192,4,4,48
Geeft: 2521835,67 356672,731550951 -0,00331534775033703 -15,2150931367566
Moet zijn: 2521835.67 -184072 346042 30400
Zou je een paar tussen waardes willen geven, zodat ik kan zoeken waar het mis gaat ?
Re: Maan Positie
Code: Selecteer alles
getX2000 -> compute_Elp_arguments():
T = 1.92445366
W1 = 1.97846340e+00
W2 = -1.52649066e-01
W3 = 3.61460556e-03
Ea = 4.53072350e+00
pomp = 1.80743955e+00
D = 5.89332560e-01
F = 1.97484880e+00
L = 2.13111247e+00
Lp = 2.72328395e+00
getX2000 -> compute_Elp_sums():
main_long = 8.12128459e-02
main_lat = 7.77709587e-02
main_dist = 3.93130604e+05
pert_longT1 = 3.30227090e-06
getX2000:
longM = 2.05968261 latM = 0.07777096
r=393130.584080 x0=-184072.92 y0=346028.78 z0=30543.33
P=2.1329e-05 Q=-2.1790e-04 sq=9.99999976e-01
Precession matrix:
1.000000e+00 -9.295209e-09 4.265841e-05
-9.295209e-09 9.999999e-01 4.357972e-04
-4.265841e-05 -4.357972e-04 9.999999e-01
Output:
2521835.67 2192-06-13 04:04:48 -184072 346042 30400
Re: Maan Positie
Problemen gevonden en opgelost.
De Geocentrische_Maan_Positie (x, Y en Z) zijn nu gelijk aan je voorbeeld waardes.
De volgende stap in het berekenen van de azimuth en altitude, en ben paragraaf 17.2 (pag. 40) aan het lezen.
Helaas snap ik hier niet veel van.
De Geocentrische_Maan_Positie (x, Y en Z) zijn nu gelijk aan je voorbeeld waardes.
De volgende stap in het berekenen van de azimuth en altitude, en ben paragraaf 17.2 (pag. 40) aan het lezen.
Helaas snap ik hier niet veel van.
Re: Maan Positie
We weten nu de positie van de maan ten opzichte van de aarde als functie van tijd.
Hier de posities voor november 2023:
Dit is de tabel in de vorm van mijn eerdere post, uitgebreid met:
R = \(\sqrt{X^2 + Y^2 + Z^2}\)
RSG = de afstand aarde-maan volgens de Sterrengids 2023
R/RSG
|R-RSG| = het absolute verschil tussen R en RSG
Onze benadering lijkt behoorlijik goed te kloppen.
Als ik de X en Y waarden plot met de aarde in de oorsprong, krijg ik dit:
Hierin is de omlooptijd van de maan ongeveer 27.3 dagen = een siderische maand.
De aarde staat wel in het centrum, maar de assen zijn niet gefixeerd aan de aarde (waarschijnlijk op een of andere manier gefixeerd aan het ecliptisch vlak).
Ik zal in https://ytliu0.github.io/starCharts/doc ... charts.pdf duiken om te zien wat er precies gebeurt, maar dat zal nog even tijd kosten (het is een nogal omvangrijk boekwerk).
Hier de posities voor november 2023:
Code: Selecteer alles
Jul. Date Date Time X Y Z R RSG R/RSG |R-RSG|
-----------------------------------------------------------------------------------------------------------------
2460249.50 2023-11-01 00:00:00 79627 375267 28352 384669 384796 0.999669 127
2460250.50 2023-11-02 00:00:00 -6746 388596 32780 390034 390147 0.999711 113
2460251.50 2023-11-03 00:00:00 -92615 382293 35525 394952 395057 0.999735 105
2460252.50 2023-11-04 00:00:00 -173827 357357 36518 399066 399170 0.999740 104
2460253.50 2023-11-05 00:00:00 -246717 315507 35779 402112 402220 0.999732 108
2460254.50 2023-11-06 00:00:00 -308148 259019 33389 403932 404038 0.999737 106
2460255.50 2023-11-07 00:00:00 -355515 190622 29482 404471 404565 0.999768 94
2460256.50 2023-11-08 00:00:00 -386755 113426 24234 403772 403838 0.999837 66
2460257.50 2023-11-09 00:00:00 -400373 30883 17868 401960 401989 0.999927 29
2460258.50 2023-11-10 00:00:00 -395513 -53271 10658 399227 399224 1.000007 3
2460259.50 2023-11-11 00:00:00 -372043 -135086 2935 395819 395798 1.000054 21
2460260.50 2023-11-12 00:00:00 -330649 -210526 -4929 392013 391993 1.000051 20
2460261.50 2023-11-13 00:00:00 -272913 -275633 -12533 388088 388088 1.000000 0
2460262.50 2023-11-14 00:00:00 -201340 -326755 -19476 384299 384325 0.999932 26
2460263.50 2023-11-15 00:00:00 -119320 -360780 -25377 380846 380889 0.999887 43
2460264.50 2023-11-16 00:00:00 -31015 -375394 -29897 377858 377897 0.999897 39
2460265.50 2023-11-17 00:00:00 58840 -369296 -32759 375386 375397 0.999971 11
2460266.50 2023-11-18 00:00:00 145206 -342370 -33766 373419 373384 1.000095 35
2460267.50 2023-11-19 00:00:00 223068 -295774 -32820 371912 371828 1.000226 84
2460268.50 2023-11-20 00:00:00 287785 -231931 -29941 370822 370705 1.000315 117
2460269.50 2023-11-21 00:00:00 335445 -154405 -25276 370140 370021 1.000320 119
2460270.50 2023-11-22 00:00:00 363163 -67683 -19096 369909 369822 1.000236 87
2460271.50 2023-11-23 00:00:00 369313 23140 -11778 370224 370195 1.000079 29
2460272.50 2023-11-24 00:00:00 353652 112728 -3768 371203 371241 0.999897 38
2460273.50 2023-11-25 00:00:00 317325 195897 4453 372948 373049 0.999730 101
2460274.50 2023-11-26 00:00:00 262748 267988 12413 375511 375654 0.999619 143
2460275.50 2023-11-27 00:00:00 193383 325182 19672 378850 379013 0.999571 163
2460276.50 2023-11-28 00:00:00 113433 364719 25854 382825 382984 0.999586 159
2460277.50 2023-11-29 00:00:00 27494 385001 30656 387197 387338 0.999636 141
2460278.50 2023-11-30 00:00:00 -59774 385584 33867 391657 391778 0.999691 121
R = \(\sqrt{X^2 + Y^2 + Z^2}\)
RSG = de afstand aarde-maan volgens de Sterrengids 2023
R/RSG
|R-RSG| = het absolute verschil tussen R en RSG
Onze benadering lijkt behoorlijik goed te kloppen.
Als ik de X en Y waarden plot met de aarde in de oorsprong, krijg ik dit:
Hierin is de omlooptijd van de maan ongeveer 27.3 dagen = een siderische maand.
De aarde staat wel in het centrum, maar de assen zijn niet gefixeerd aan de aarde (waarschijnlijk op een of andere manier gefixeerd aan het ecliptisch vlak).
Ik zal in https://ytliu0.github.io/starCharts/doc ... charts.pdf duiken om te zien wat er precies gebeurt, maar dat zal nog even tijd kosten (het is een nogal omvangrijk boekwerk).
Re: Maan Positie
Er bestaan diverse benaderingsalgorithmen voor de verschillende stappen in de berekening.
Ik zal in de Kerstvakantie kijken in hoeverre de nauwkeurigheid is op te voeren.
Hier wat voorlopige resultaten:
DateTime.h = datum en tijd berekeningen
MoonPosition.h = Yuk Tung Liu's Lunar Ephemeris
Moon.c = hoofdprogramma
OUTPUT:
Ter vergelijking: (tussen rechte haken steeds de waarden uit bovenstaande tabel):
https://hemel.waarnemen.com/applets/maan.cgi
Utrecht, 9-12-2023, 10:20 uur wintertijd = 9:20 UT:
alt = 22.9° [23.00]
Az = 192.8° [194.56]
R = 391818 km [391887]
opkomst: 4:23 uur, Az = 110° [110.31]
doorgang: 9:29 uur, alt = 24.8° [24.08]
(noot: de doorgang op deze site wijkt af van de Sterrengids, zie hieronder)
ondergang: 14:21 uur, Az = 246° [246.47]
Sterrengids 2023:
opkomst: Az = 109.8° [110.31]
doorgang: alt = 24.0° [24.08]
ondergang: Az = 246.3° [246.47]
0h UT:
R = 393614 km [393682]
RA = 13:43.9 h [13:43]
Dec = -11° 11' [-11 49']
Ik zal in de Kerstvakantie kijken in hoeverre de nauwkeurigheid is op te voeren.
Hier wat voorlopige resultaten:
DateTime.h = datum en tijd berekeningen
Code: Selecteer alles
//---------------------------------------------------------------------------
// JULIAN DATE:
//---------------------------------------------------------------------------
// truncate x = drop fractional part of x:
double trunc(double x)
{
double sign = 1.0;
if(x<0.0){
sign = -1.0;
x = -x;
}
return(sign*floor(x));
}
// calculate Julian Date for date D-M-Y time h:m:s
double JulianDate(int D, int M, int Y, double h, double m, double s)
{
double day=(double)D + (h+m/60.0+s/3600.0)/24.0;
double year, month, tmp, a, b, c;
if(M<=2){
year = Y-1;
month = M+12;
}
else{
year = Y;
month = M;
}
// 15-10-1582 = start Gregorian calendar:
if((Y<1582) || (Y==1582 && M<10) || (Y==1582 && M==10 && D<15)) // before Gregorian cal:
a = 0;
else{ // Gregorian cal:
tmp = trunc(year/100.0);
a = 2 - tmp + floor(tmp/4.0);
}
if(year<0)
b = trunc(365.25*year-0.75);
else
b = trunc(365.25*year);
c = trunc(30.6001 * (month+1));
return(a + b + c + day + 1720994.5);
}
//---------------------------------------------------------------------------
// SIDEREAL TIME:
// (Fundamental Astronomy - 5th edition - 2007 - H. Karttunen e.a.)
//---------------------------------------------------------------------------
// gmst at 0 hour UT:
// in: J, such that fraction J = 0.5
// out: GMST at 0:00:00 hour
double gmstZeroUT(double J)
{
double t=J-2451545.0; // start of J2000
double T=t/36525.0; // Julian century
double T2=T*T;
double T3=T2*T;
double gmst = (24110.54841 + 8640184.812866*T + 0.093104*T2 - 0.0000062*T3) / 3600.0;
gmst=fmod(gmst,24.0);
if(gmst<0.0) gmst+=24.0;
return(gmst);
}
// get sidereal time at latitude lambda, date D-M-YYYY at h:m:s UT:
double siderealTime(double lambda, int D, int M, int Y, int h, int m, double s)
{
double J = JulianDate(D,M,Y,0,0,0);
double gmst=gmstZeroUT(J);
double H = (double)h + (double)m/60.0 + s/3600.0; // time in hours UT
double DH = 0.002743055555555*H;
double ST = gmst + DH + H + lambda/15.0;
ST=fmod(ST, 24.0);
if(ST<0.0) ST+=24.0;
return(ST);
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
Code: Selecteer alles
//---------------------------------------------------------------------------
// Lunar Ephemeris
//---------------------------------------------------------------------------
// Source:
// https://ytliu0.github.io/starCharts/docs/star_charts.pdf
// chapter 17
// Based on truncated ELP/MPP02 series using the following parameters:
// corr = 1 = Parameters fitted to JPL ephemeris DE405/DE406
// AthU = 1.454441e-04 radians = 30.000000 arcsecs
// AthV = 1.454441e-04 radians = 30.000000 arcsecs
// AthR = 100.000000 km
// tau = 50.000000
// Total number of terms in the original series: 35901
// Total number of terms in the truncated series: 42
// Statistical evaluation:
// Estimated error of the truncated series by comparing it to the original series
// at 10000 randomly chosen T between -50.000000 and 10.000000:
// Maximum deviation found in ecliptic longitude: 213.496798 arcsecs
// Root mean square deviation in ecliptic longitude: 47.665830 arcsecs
// Maximum deviation found in ecliptic latitude: 158.214315 arcsecs
// Root mean square deviation in ecliptic latitude: 34.778046 arcsecs
// Maximum deviation found in geocentric distance: 336.995129 km
// Root mean square deviation in geocentric distance: 85.856277 km
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
// standard C libs:
#include <stdio.h>
#include <math.h>
//---------------------------------------------------------------------------
// constants:
extern const double PI;
const double TPI = 2.0*PI;
//---------------------------------------------------------------------------
// Arguments for the ELP/MPP02 series:
double W1, D, F, L, Lp;
// coefficients for the ELP/MPP02 main series:
int n_main_long=19;
int i_main_long[19][4]={
{0,2,0,0},{0,-2,1,0},{0,0,1,0},{0,2,1,0},{0,0,2,0},{0,0,3,0},{0,0,-1,1},{0,0,0,1},
{0,0,1,1},{1,0,0,0},{2,0,-1,-1},{2,0,0,-1},{2,0,-2,0},{2,0,-1,0},{2,-2,0,0},{2,0,0,0},
{2,0,1,0},{4,0,-2,0},{4,0,-1,0}};
double A_main_long[19]={
-1.99547249801829270800e-03, 1.91662856547817478400e-04, 1.09759808626191612800e-01, -2.18649080896658580300e-04,
3.72834182278018950000e-03, 1.75133187917569992000e-04, -7.14234214377495904400e-04, -3.23088338567132336800e-03,
-5.30290959258585535000e-04, -6.05959497615181906000e-04, 9.95981590587821435600e-04, 7.98626872936074750500e-04,
1.02613505456892778600e-03, 2.22356802322449947900e-02, 2.67505934260747291300e-04, 1.14896669562630695400e-02,
9.30629894589579480800e-04, 1.49189651057941429400e-04, 1.86313093175213899100e-04};
int n_main_lat=10;
int i_main_lat[10][4]={
{0,1,0,0},{0,-1,1,0},{0,1,1,0},{0,-1,2,0},{0,1,2,0},{2,-1,-1,0},{2,1,-1,0},{2,-1,0,0},{2,1,0,0},{2,-1,1,0}};
double A_main_lat[10]={
8.95026190647627839500e-02, 4.84665164815963230000e-03, 4.89742857492255542000e-03, 1.53975229251259538900e-04,
3.00157605100467080500e-04, 8.07574106044906285200e-04, 9.67124565462755482600e-04, 3.02355257133907857600e-03,
5.68495899772767385500e-04, 1.61720283648969048700e-04};
int n_main_dist=12;
int i_main_dist[12][4]={
{0,0,0,0},{0,0,1,0},{0,0,2,0},{0,0,-1,1},{0,0,1,1},{1,0,0,0},
{2,0,-1,-1},{2,0,0,-1},{2,0,-2,0},{2,0,-1,0},{2,0,0,0},{2,0,1,0}};
double A_main_dist[12]={
3.85000529032284044600e+05, -2.09053549352050904400e+04, -5.69925115335094687900e+02, -1.29620222172050574700e+02,
1.04755292674270279460e+02, 1.08742701272463051500e+02, -1.52137809425890736800e+02, -2.04586116460873086000e+02,
2.46158474853557947900e+02, -3.69911089639807642000e+03, -2.95596755797241485200e+03, -1.70733077124770574100e+02};
double main_long, main_lat, main_dist;
double pert_longT1;
//---------------------------------------------------------------------------
// restrict x to [-pi,pi>
double mod2pi(double x)
{
return(x - TPI*floor((x + PI)/TPI));
}
// Compute the lunar and planetary arguments used in the ELP/MPP02 series:
// set Elp arguments: W1, D, F, L, Lp
void compute_Elp_arguments(double T)
{
double T2 = T*T;
double T3 = T*T2;
double T4 = T2*T2;
W1 = -2.4728412163487064 + mod2pi(8399.6847300719273*T) + mod2pi(-3.3189520425500949e-05*T2)
+ mod2pi(3.1102494491060616e-08*T3) + mod2pi(-2.0328237648922845e-10*T4);
double W2 = 1.4547895404440774 + mod2pi(70.993305448479248*T) + mod2pi(-0.00018548192818782712*T2)
+ mod2pi(-2.1961637966359414e-07*T3) + mod2pi(1.0327016221314225e-09*T4);
double W3 = 2.1824388474237688 + mod2pi(-33.781427419672326*T) + mod2pi(3.0816644950982666e-05*T2)
+ mod2pi(3.6447710769397583e-08*T3) + mod2pi(-1.738541860458796e-10*T4);
double Ea = 1.7534699452640699 + mod2pi(628.30758508103168*T) + mod2pi(-9.7932363584126268e-08*T2)
+ mod2pi(4.3633231299858238e-11*T3) + mod2pi(7.2722052166430391e-13*T4);
double pomp = 1.7965956331206001 + mod2pi(0.0056298669711442699*T) + mod2pi(2.5659491293243858e-06*T2)
+ mod2pi(-5.7275888286280579e-10*T3) + mod2pi(5.5166948773454105e-11*T4);
// Mean longitude of the Moon:
W1 = mod2pi(W1);
// Arguments of Delaunay:
D = mod2pi(W1 - Ea + PI);
F = mod2pi(W1 - W3);
L = mod2pi(W1 - W2);
Lp = mod2pi(Ea - pomp);
}
//---------------------------------------------------------------------------
// Sum the ELP/MPP02 main and perturbation series:
void compute_Elp_sums(void)
{
int i;
double phase;
// main sums (long and lat: sin-sums; dist: cos-sum):
main_long=0.0;
for(i=0;i<n_main_long;i++){
phase=i_main_long[i][0]*D + i_main_long[i][1]*F + i_main_long[i][2]*L + i_main_long[i][3]*Lp;
main_long += A_main_long[i]*sin(phase);
}
main_lat=0.0;
for(i=0;i<n_main_lat;i++){
phase=i_main_lat[i][0]*D + i_main_lat[i][1]*F + i_main_lat[i][2]*L + i_main_lat[i][3]*Lp;
main_lat += A_main_lat[i]*sin(phase);
}
main_dist=0.0;
for(i=0;i<n_main_dist;i++){
phase=i_main_dist[i][0]*D + i_main_dist[i][1]*F + i_main_dist[i][2]*L + i_main_dist[i][3]*Lp;
main_dist += A_main_dist[i]*cos(phase);
}
// perturbations sums: pert_longT1 is the only non-zero for our approximation:
pert_longT1 = 8.1293558048447e-06 * sin(Lp);
}
//---------------------------------------------------------------------------
// Calculate the Moon's geocentric X,Y,Z coordinates with respect to
// J2000.0 mean ecliptic and equinox.
void getX2000(double JD, double &X, double &Y, double &Z)
{
// T is the TDB Julian century from J2000.0 = (JD - 2451545)/36525
double T=(JD - 2451545.0)/36525.0;
double T2 = T*T;
double T3 = T*T2;
double T4 = T2*T2;
double T5 = T2*T3;
compute_Elp_arguments(T);
compute_Elp_sums();
// Moon's longitude, latitude and distance
double longM = W1 + main_long + mod2pi(pert_longT1*T);
double latM = main_lat;
// ra0 = a0(DE405)/a0(ELP) = 384747.961370173/384747.980674318 = 0.99999994982652041.
double r = 0.99999994982652041*main_dist;
double x0 = r*cos(longM)*cos(latM);
double y0 = r*sin(longM)*cos(latM);
double z0 = r*sin(latM);
// Precession matrix
double P = 0.10180391e-4*T + 0.47020439e-6*T2 - 0.5417367e-9*T3 - 0.2507948e-11*T4 + 0.463486e-14*T5;
double Q = -0.113469002e-3*T + 0.12372674e-6*T2 + 0.12654170e-8*T3 - 0.1371808e-11*T4 - 0.320334e-14*T5;
double sq = sqrt(1 - P*P - Q*Q);
double p11 = 1 - 2*P*P;
double p12 = 2*P*Q;
double p13 = 2*P*sq;
double p21 = 2*P*Q;
double p22 = 1-2*Q*Q;
double p23 = -2*Q*sq;
double p31 = -2*P*sq;
double p32 = 2*Q*sq;
double p33 = 1 - 2*P*P - 2*Q*Q;
bool show=false;
if(show){
printf("precession matrix: -------------------\n");
printf("%16.8lf%16.8lf%16.8lf\n",p11,p12,p13);
printf("%16.8lf%16.8lf%16.8lf\n",p21,p22,p23);
printf("%16.8lf%16.8lf%16.8lf\n",p31,p32,p33);
printf("--------------------------------------\n");
}
// Finally, calculate the return values X, Y and Z:
// (= the components of the position vector; wrt J2000.0 mean ecliptic and equinox)
X = p11*x0 + p12*y0 + p13*z0;
Y = p21*x0 + p22*y0 + p23*z0;
Z = p31*x0 + p32*y0 + p33*z0;
}
//---------------------------------------------------------------------------
// Demo:
//---------------------------------------------------------------------------
// Demonstration output for our ELP/MPP02 series approximation
// Compare with the high precision version on page 6 in
// https://github.com/ytliu0/ElpMpp02/blob/master/docs/ElpMpp02.pdf
void showTableHeader(void)
{
printf("\n");
printf("Lunar geocentric rectangular coordinates (in km), fitted to DE405/DE406\n");
printf("-------------------------------------------------------------------------------------------\n");
printf("Jul. Date Date Time X Y Z R\n");
printf("-------------------------------------------------------------------------------------------\n");
}
void showMoonPosition(int D, int M, int Y, int h, int m, double s)
{
double JD=JulianDate(D,M,Y,h,m,s); // calculate Julian Date
double x,y,z,R; // lunar coordinates and distance
getX2000(JD,x,y,z); // calculate lunar coordinates on JD
R=sqrt(x*x+y*y+z*z);
printf("%.2lf %5d-%02d-%02d %02d:%02d:%02.0lf %10.0lf %10.0lf %10.0lf %10.0lf\n",JD,Y,M,D,h,m,s,x,y,z,R);
}
// Table at bottom of pag. 6:
void showDemoTable(void)
{
showTableHeader();
showMoonPosition(13,6,2192,4,4,48);
showMoonPosition(7,12,1490,19,55,12);
showMoonPosition(16,6,789,11,45,36);
showMoonPosition(25,12,87,3,36,0);
showMoonPosition(3,7,-614,19,26,24);
printf("\n");
}
//---------------------------------------------------------------------------
// <<< end of Yuk Tung Liu's Lunar Ephemeris
//---------------------------------------------------------------------------
Code: Selecteer alles
//---------------------------------------------------------------------------
// standard C libs:
#include <stdio.h>
#include <math.h>
// local libs:
#include "DateTime.h"
#include "MoonPosition.h"
#pragma hdrstop
//---------------------------------------------------------------------------
// global constants:
const double PI = 3.14159265358979323846264338327950;
const double p = PI/180.0;
// global variables:
double Day, Month, Year, hours, minutes, seconds, JD; // current date and time (UT)
double latitude, longitude; // current location on Earth
//---------------------------------------------------------------------------
// reduce x modulo mod: x goes to interval [0, mod>
double reduce(double x, double mod)
{
x=fmod(x,mod);
if(x<0.0) x+=mod;
return(x);
}
//---------------------------------------------------------------------------
// set latitude and longitude of the current observer location:
// lat and lon in degrees
void setLocation(double lat, double lon)
{
latitude = lat;
longitude = lon;
}
// set date and time:
double setDateTime(double D, double M, double Y, double h, double m, double s)
{
Day=D; Month=M; Year=Y;
hours=h; minutes=m; seconds=s;
JD = JulianDate(D,M,Y,h,m,s);
return(JD);
}
//---------------------------------------------------------------------------
// calculate geocentric Moon position
// pre: JD = Julian Date has to be known
// out: vec = geocentric Moon position as vector [x, y, z] in km
// R = distance Earth to Moon in km
void getMoonPosition(double *vec, double &R)
{
double x,y,z;
getX2000(JD, x,y,z);
vec[0]=x;
vec[1]=y;
vec[2]=z;
R=sqrt(x*x+y*y+z*z);
}
//---------------------------------------------------------------------------
// convert geocentric to topocentric location:
// in: lat = latitude [deg]
// st = sidereal time [hr]
// u = geocentric location vector [x, y, z] in km
// out: v = topocentric location vector [x, y, z] in km
void geocentric2topocentric(double lat, double st, double *u, double *v)
{
double a=6378.1366; // major earth axis in km
double f=1.0/298.25642;
double phi=lat*PI/180.0;
double cp = cos(phi);
double fsp = (1.0-f)*sin(phi);
double ts = 15.0*st;
double C=sqrt(cp*cp + fsp*fsp);
double S=(1-f)*(1-f)*C;
v[0] = u[0] - a*C*cp*cos(ts);
v[1] = u[1] - a*C*cp*sin(ts);
v[2] = u[2] - a*S*sin(phi);
}
//---------------------------------------------------------------------------
// get obliquity of the ecliptic at Julian Date JD:
// out: epsilon in degrees
double calcEpsilon(double JD)
{
double t=JD-2451545.0; // start of J2000
double T=t/36525.0; // Julian century
double T2=T*T;
double T3=T2*T;
double Omega = 125.04452 - 1934.136261*T + 0.0020708*T2 + T3/450000;
Omega=reduce(Omega,360.0);
double L = reduce(280.4665 + 36000.7698*T, 360.0);
double Lacc = reduce(218.3165 + 481267.8813*T, 360.0);
// --- phi correction in seconds --- skipped: far below current precision
//double Dphi = -17.20*sin(Omega*p) - 1.32*sin(2.0*L*p) - 0.23*sin(2.0*Lacc*p) + 0.21*sin(2.0*Omega*p);
//printf("Dphi = %.8lf\n",Dphi);
//Dphi/=3600.0;
//printf("Dphi = %.8lf\n",Dphi);
double Deps = 9.20*cos(Omega*p) + 0.57*cos(2.0*L*p); + 0.10*cos(2.0*Lacc*p) - 0.09*cos(2.0*Omega*p);
Deps/=3600.0;
double U=T/100.0;
double Ufac[11]={2.45, 5.79, 27.87, 7.12, -39.05, -249.67, -51.38, 1999.25, -1.55, -4680.93, 3600.0*23.43929111111111111};
double eps0=0.0;
for(int i=0; i<11; i++)
eps0 = eps0*U + Ufac[i];
eps0/=3600.0;
double eps = eps0 + Deps;
return(eps);
}
//---------------------------------------------------------------------------
// convert ecliptic to equatorial reference system:
// in: lambda = ecliptic longitude [deg]
// beta = ecliptic latitude [deg]
// epsilon = sidereal time [deg]
// out: alpha = Right Ascension [hours]
// delta = declination [deg]
void ecliptic2equatorial(double lambda, double beta, double epsilon, double &alpha, double &delta)
{
double l=lambda*p;
double b=beta*p;
double e=epsilon*p;
delta = asin(cos(e)*sin(b)+sin(e)*cos(b)*sin(l))/p;
alpha = atan2(cos(b)*sin(l)*cos(e)-sin(b)*sin(e),cos(b)*cos(l))/p;
alpha = reduce(alpha, 360.0)/15.0; // from degrees to [0,24> hours
}
//---------------------------------------------------------------------------
// in: alpha = right ascension [hours]
// delta = declination [deg]
// ts as time [hours]
// lat = latitude [deg]
// out: Az = Azimuth [deg]
// alt = altitude [deg]
void topoEquatorial2AzAlt(double alpha, double delta, double ts, double lat, double &Az, double &alt)
{
double p = PI/180.0;
double H = ts - alpha;
H*=15;
H = reduce(H, 360.0);
alt = asin( sin(delta*p)*sin(lat*p) + cos(delta*p)*cos(H*p)*cos(lat*p) );
double sA = -cos(delta*p)*sin(H*p)/cos(alt);
double cA = ( sin(delta*p)*cos(lat*p) - cos(delta*p)*cos(H*p)*sin(lat*p) );
Az = atan2(sA, cA)/p;
Az = reduce(Az, 360.0);
alt /= p;
}
//---------------------------------------------------------------------------
// OUTPUT TO SCREEN:
//---------------------------------------------------------------------------
// print hour H as hh:mm
void printHM(double H)
{
int min=(int)(60.0*H+0.5);
printf("%4d:%02d",min/60, min%60);
}
// print hour H as hh:mm:ss
void printHMS(double H)
{
int t=(int)(3600.0*H+0.5);
int sec=t%60;
t/=60;
printf("%2d:%02d:%02d",t/60, t%60, sec);
}
// show table header:
void showLocationHeader(void)
{
printf("---------------------------------------------------------------------------------------------\n");
printf(" Date UnivTime ST R RA Dec Az Alt \n");
//printf("dd-mm-yyyy hh:mm:ss hh:mm:ss km hh:mm deg deg deg \n");
printf("---------------------------------------------------------------------------------------------\n");
}
// show position of Moon as line in table:
void showMoon(double D, double M, double Y, double h, double m, double s)
{
setDateTime(D,M,Y,h,m,s);
// date and time:
printf("%2.0lf-%2.0lf-%.0lf", Day, Month, Year);
printf("%5.0lf:%02.0lf:%02.0lf", hours, minutes, seconds);
// calculate sidereal time:
double st=siderealTime(longitude, Day,Month,Year, hours,minutes,seconds);
printf(" "); printHMS(st);
// get geocentric position:
double u[3], v[3], R;
getMoonPosition(u, R);
printf("%11.0lf",R); // R = distance Earth to Moon
// convert geocentric to topocentric position:
geocentric2topocentric(latitude, st, u, v);
// topocentric polar representation:
double r=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
double lambda=atan2(v[1], v[0])/p;
double beta = asin(v[2]/r)/p;
// obliquity of the ecliptic:
double epsilon=calcEpsilon(JD);
// Right ascension and Declination:
double alpha,delta;
ecliptic2equatorial(lambda, beta, epsilon, alpha, delta);
printf(" "); printHM(alpha);
printf("%11.2lf", delta);
// Azimuth and Altitude:
double Azimuth, Altitude;
topoEquatorial2AzAlt(alpha, delta, st, 52.085555555, Azimuth, Altitude);
printf("%12.2lf", Azimuth);
printf("%12.2lf\n", Altitude);
}
//---------------------------------------------------------------------------
// MAIN PROGRAM:
//---------------------------------------------------------------------------
#pragma argsused
int main(int argc, char* argv[])
{
// set location: Utrecht = 52.086 N, 5.129 E (= standaardlocatie Sterrengids 2023)
setLocation(52.086, 5.129);
showLocationHeader();
showMoon(9,12,2023, 0,0,0);
showMoon(9,12,2023, 3,23,0); // opkomst
showMoon(9,12,2023, 8,29,0); // doorgang
showMoon(9,12,2023, 9,20,0);
showMoon(9,12,2023, 13,21,0); // onder
return 0;
}
//---------------------------------------------------------------------------OUTPUT:
Code: Selecteer alles
---------------------------------------------------------------------------------------------
Date UnivTime ST R RA Dec Az Alt
---------------------------------------------------------------------------------------------
9-12-2023 0:00:00 5:30:27 393682 13:43 -11.82 72.59 -29.30
9-12-2023 3:23:00 8:54:00 393041 13:50 -12.68 110.31 -0.46
9-12-2023 8:29:00 14:00:51 392054 14:00 -13.83 180.34 24.08
9-12-2023 9:20:00 14:51:59 391887 14:01 -14.00 194.56 23.00
9-12-2023 13:21:00 18:53:39 391093 14:08 -14.76 246.47 -0.69
Ter vergelijking: (tussen rechte haken steeds de waarden uit bovenstaande tabel):
https://hemel.waarnemen.com/applets/maan.cgi
Utrecht, 9-12-2023, 10:20 uur wintertijd = 9:20 UT:
alt = 22.9° [23.00]
Az = 192.8° [194.56]
R = 391818 km [391887]
opkomst: 4:23 uur, Az = 110° [110.31]
doorgang: 9:29 uur, alt = 24.8° [24.08]
(noot: de doorgang op deze site wijkt af van de Sterrengids, zie hieronder)
ondergang: 14:21 uur, Az = 246° [246.47]
Sterrengids 2023:
opkomst: Az = 109.8° [110.31]
doorgang: alt = 24.0° [24.08]
ondergang: Az = 246.3° [246.47]
0h UT:
R = 393614 km [393682]
RA = 13:43.9 h [13:43]
Dec = -11° 11' [-11 49']
Re: Maan Positie
Alvast bedankt voor de moeite !!!
Re: Maan Positie
Hieronder 4 bestanden in C:
Moon.cpp: de diverse rekenvoorschriften uit Yuk Tung Li - Calculations in Star Charts
Moon.h: bevat alle globale constanten en variabelen voor de (tussen-)resultaten
AstroLibrary.cpp: enkele algemene astronomie functies
MoonPosition.cpp: de Lunar Ephemeris: geeft de J2000 ecliptische coordinaten [X, Y, Z] van de maan (zoals hierboven al beschreven).
Moon.cpp roept in de startfunctie main() een demonstratie aan van het gebruik, waarbij:
setObserver() en getMoonPosition() alle invoergegevens krijgen
en
showMoonData() alle tussenresultaten en showShort() de belangrijkste resultaten op het scherm tonen.
Moon.cpp
Moon.h
AstroLibrary.cpp
MoonPosition.cpp
Output:
Hier de ouput van bovenstaande:
Ter vergelijk:
Sterrengids 2023 (22-12-2023):
0 uur UT: Geocentrische RA = 1 uur 44.7 min, decl = 11°23', R = 376064 km
Ondergang: 2:30, Az = 289.6 graden
Opkomst: 12:31, Az = 66.3 graden
Doorgang: 19:14, Altitude = 53.5 graden
En via deze site: https://ytliu0.github.io/starCharts/:
Heel wat informatie, stel gerust vragen over onduidelijkheden in bovenstaande.
Ik hoop dat het al voldoende nauwkeurig is, en dat de code eenvoudig is om te zetten naar Delphi.
PS: hier nog de tussenresultaten van 22-12-2023 0 uur UT, verkregen via showMoonData();
Moon.cpp: de diverse rekenvoorschriften uit Yuk Tung Li - Calculations in Star Charts
Moon.h: bevat alle globale constanten en variabelen voor de (tussen-)resultaten
AstroLibrary.cpp: enkele algemene astronomie functies
MoonPosition.cpp: de Lunar Ephemeris: geeft de J2000 ecliptische coordinaten [X, Y, Z] van de maan (zoals hierboven al beschreven).
Moon.cpp roept in de startfunctie main() een demonstratie aan van het gebruik, waarbij:
setObserver() en getMoonPosition() alle invoergegevens krijgen
en
showMoonData() alle tussenresultaten en showShort() de belangrijkste resultaten op het scherm tonen.
Moon.cpp
Code: Selecteer alles
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
// CALCULATIONS IN STAR CHARTS
// Yuk Tung Liu, 24-3-2023
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
#pragma hdrstop
//-------------------------------------------------------------------------------
// standard C libs:
#include <stdio.h>
#include <math.h>
// local libs:
#include "Moon.h"
#include "AstroLibrary.cpp"
#include "MoonPosition.cpp"
//-------------------------------------------------------------------------------
// Yuk - 7.5 Nutation Matrix (pg. 13)
//-------------------------------------------------------------------------------
// calculate nutation parameters [all in degrees]:
// uses: JD = Julian Date
// out: EpsilonMean = mean obliquity of the ecliptic of date (= eA in the text)
// EpsilonTrue = true obliquity of the ecliptic of date
// DeltaEpsilon = EpsilonTrue - EpsilonMean
// DeltaPsi = nutation in longitude
// EqEq = Equation of the Equinoxes
// notes:
// when used in time span between 3000 BCE and 3000 CE (-50 < T < 10):
// maximum error in (30) is 0.17" and the rms error is 0.044"
// maximum error in (31) is about 0.11" and the rms error is 0.032"
void setNutationParameters()
{
double L, Lacc, F, D, Om;
double T=(JD-2451545.0)/36525.0;
double T2=T*T;
double T3=T*T2;
double T4=T*T3;
double T5=T*T4;
// mean obliquity of the ecliptic of date:
EpsilonMean = (84381.406 - 46.836769*T - 0.0001831*T2 + 0.00200340*T3 - 0.000000576*T4 - 0.0000000434*T5)/3600.0; // (28)
// L = the mean anomaly of the Moon
// Lacc = is the mean anomaly of the Sun
// F = ML - Om, where ML = the mean longitude of the Moon
// D = the mean elongation of the Moon from the Sun
// Om = is the longitude of the ascending node of the Moon's mean orbit on the
// ecliptic measured from the mean equinox of date
L = (485868.249036 + 1717915923.2178*T + 31.8792*T2 + 0.051635*T3 - 0.00024470*T4)/3600.0; // (32)
Lacc = (1287104.79305 + 129596581.0481*T - 0.5532*T2 + 0.000136*T3 - 0.00001149*T4)/3600.0; // (33)
F = (335779.526232 + 1739527262.8478*T - 12.7512*T2 - 0.001037*T3 + 0.00000417*T4)/3600.0; // (34)
D = (1072260.70369 + 1602961601.2090*T - 6.3706*T2 + 0.006593*T3 - 0.00003169*T4)/3600.0; // (35)
Om = (450160.398036 - 6962890.5431*T + 7.4722*T2 + 0.007702*T3 - 0.00005939*T4)/3600.0; // (36)
DeltaPsi= -(17.2064161 + 0.0174666*T)*sin(Om*D2R) - 1.3170906*sin((2.0*F-2.0*D+2.0*Om)*D2R)
-0.2276413*sin((2.0*F+2.0*Om)*D2R) + 0.2074554*sin(2.0*Om*D2R) + 0.1475877*sin(Lacc*D2R)
-0.0516821*sin((Lacc + 2.0*F - 2.0*D + 2*Om)*D2R) + 0.0711159*sin(L*D2R); // (30)
DeltaEpsilon = 9.2052331*cos(Om*D2R) + 0.5730336*cos((F-D+Om)*2.0*D2R)
+ 0.0978459*cos((F+Om)*2.0*D2R) - 0.0897492*cos(2.0*Om*D2R); // (31)
DeltaPsi/=3600.0; // seconds to degrees
DeltaEpsilon/=3600.0; // seconds to degrees
EpsilonTrue = EpsilonMean + DeltaEpsilon; // {29)
EqEq = DeltaPsi * cos(EpsilonMean*D2R);
}
//-------------------------------------------------------------------------------
// SIDEREAL TIMES:
//-------------------------------------------------------------------------------
// Yuk - 2 Sidereal Time, UT1, UTC, and Civil Time (pg. 3)
// get GMST [hours]:
// DU = JD - 2451545
// T = DU / 36525
// Theta(DU) = (24*(0.7790572732640 + 1.00273781191135448*DU)) mod 24 // (2)
// Eprec(T) = -0.014506 - 4612.16534*T - 1.3915817*T^2 + 4.4E-7*T^3 + 2.9956E-5*T^4 // (4)
// GMST(JD) = Theta(DU) - Eprec(T) // (3)
double getGMST(void){
double DU, DU0, T, h, fp, gmst;
DU = JD-2451545.0;
T=DU/36525.0;
DU0=trunc(DU-0.5)+0.5;
fp=frac(DU);
if(fp<0.5)
h = 24.0*fp + 12.0; // h = civil hours past on current day
else
h = 24.0*fp - 12.0;
gmst = 6.697374558336001 + 0.06570748587250752*DU0 + 1.00273781191135448*h
+ 2.686296296296296E-7 + 0.08541030618518518*T + 2.577003148148148E-5*T*T; // (6)
gmst = upmod(gmst, 24.0);
return(gmst);
}
// set GMST, GAST, LMST, LAST [in hours] (GAST and LAST including nutation correction)
// uses: JD = Julian Date
// EqEq = Equation of the Equinoxes
// Mlongitude = longitude of observer on Earth [deg]
void setSiderealTimes(void)
{
GMST = getGMST(); // (6)
GAST = upmod(GMST + EqEq, 24.0); // (3)
LMST = upmod(GMST + Mlongitude/15.0, 24.0); // (1)
LAST = upmod(GAST + Mlongitude/15.0, 24.0); // (1)
}
//-------------------------------------------------------------------------------
// PRECESSION:
//-------------------------------------------------------------------------------
void setPrecessionMatrix(void)
{
double T = (JD-2451545.0)/36525.0;
double T2=T*T;
double T3=T*T2;
double a = 174.8763838890*D2R + (0.035360*T2 - 869.80890*T)*SEC2RAD;
double b = (6.00E-6*T3 - 1.111130*T2 - 5029.09660*T)*SEC2RAD - a;
double e = (6.00E-5*T3 - 0.033020*T2 + 47.00290*T)*SEC2RAD;
Mprec[0][0] = cos(b)*cos(a) - cos(e)*sin(b)*sin(a);
Mprec[0][1] = cos(e)*sin(b)*cos(a) + cos(b)*sin(a);
Mprec[0][2] = sin(e)*sin(b);
Mprec[1][0] = -sin(b)*cos(a) - cos(e)*cos(b)*sin(a);
Mprec[1][1] = cos(e)*cos(b)*cos(a) - sin(b)*sin(a);
Mprec[1][2] = sin(e)*cos(b);
Mprec[2][0] = sin(e)*sin(a);
Mprec[2][1] = -sin(e)*cos(a);
Mprec[2][2] = cos(e);
}
//-------------------------------------------------------------------------------
// ECLIPTIC TO EQUATORIAL REFERENCE SYSTEM:
//-------------------------------------------------------------------------------
// set rotation matrix for conversion ecliptic to equatorial system:
void setMecl2eqt(void)
{
double T = (JD-2451545.0)/36525.0;
double eps = ( 23.43929111-(46.8150+(0.00059-0.001813*T)*T)*T/3600.0 ) * D2R;
Mecl2eqt[0][0]=1;
Mecl2eqt[0][1]=0;
Mecl2eqt[0][2]=0;
Mecl2eqt[1][0]=0;
Mecl2eqt[1][1]=cos(eps);
Mecl2eqt[1][2]=-sin(eps);
Mecl2eqt[2][0]=0;
Mecl2eqt[2][1]=sin(eps);
Mecl2eqt[2][2]=cos(eps);
}
//-------------------------------------------------------------------------------
// COORDINATES FOR OBSERVER ON EARTH:
//-------------------------------------------------------------------------------
// Yuk - 8 Topocentric Coordinates (pg. 13)
// assume: height = 0m above sea level
void setXgeotopo(void)
{
double phi = Mlatitude*D2R;
double ts = (LMST*15)*D2R;
double a = 6378.1366; // [km]
double f = 1.0/298.25642;
double cx = cos(phi);
double cy = (1.0-f)*sin(phi);
double C = 1.0/sqrt(cx*cx + cy*cy);
double S = (1.0-f)*(1.0-f) * C;
// position of observer on earth:
Xgeo[0] = a*C*cos(phi)*cos(ts);
Xgeo[1] = a*C*cos(phi)*sin(ts);
Xgeo[2] = a*S*sin(phi);
// topocentric position of the moon:
Xtopo[0] = Xeqt[0] - Xgeo[0];
Xtopo[1] = Xeqt[1] - Xgeo[1];
Xtopo[2] = Xeqt[2] - Xgeo[2];
}
//-------------------------------------------------------------------------------
// AZIMUTH AND ALTITUDE:
//-------------------------------------------------------------------------------
// Yuk - 10 Horizontal Coordinates (pg. 16)
// in: alpha = right ascension [hours]
// delta = declination [deg]
// ts as time [hours]
// lat = latitude [deg]
// out: Azimuth [deg]
// Altitude [deg]
void topo2AzAlt(double alpha, double delta, double ts, double lat)
{
double H = ts - alpha; // H = hour angle
H *= 15.0;
H = upmod(H, 360.0);
AltitudeRaw = asin( sin(delta*D2R)*sin(lat*D2R) + cos(delta*D2R)*cos(H*D2R)*cos(lat*D2R) ); // rad
double sA = -cos(delta*D2R)*sin(H*D2R)/cos(AltitudeRaw);
double cA = ( sin(delta*D2R)*cos(lat*D2R) - cos(delta*D2R)*cos(H*D2R)*sin(lat*D2R) )/cos(AltitudeRaw);
Azimuth = atan2(sA, cA)*R2D; // deg
Azimuth = upmod(Azimuth, 360.0);
AltitudeRaw *= R2D; // deg
}
//-------------------------------------------------------------------------------
// ATMOSPHERIC REFRACTION
//-------------------------------------------------------------------------------
// Yuk - 11 Atmospheric Refraction
void setRefraction(void)
{
DeltaAltitude = 1.02*(Mpressure/101.0)*(283.0/Mtemperature) /
(tan( (AltitudeRaw + 10.3/(AltitudeRaw+5.11) )*D2R)*60.0); // (45)
Altitude = AltitudeRaw + DeltaAltitude;
}
void showRefraction(void)
{
printf("Atmospheric Refraction:\n");
printf("Delta Altitude = %13.10lf = ", DeltaAltitude); printdeg(DeltaAltitude);
printf("\nObserved Altitude = %13.10lf = ", Altitude);
printdeg(Altitude); printf("\n\n");
}
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
// SET OBSERVER LOCATION AND WEATHER:
// longitude: East=positive
// latitude: North=positive
// pressure in kPa
// temperature in degrees Celsius
void setObserver(double longitude, double latitude, double pressure, double temperature)
{
Mlongitude = longitude;
Mlatitude = latitude;
Mpressure = pressure;
Mtemperature = temperature+273.15;;
}
// GIVEN OBSERVER DATE, TIME AND PLACE, CALCULATE THE MOON'S POSITION:
// time in UT
void getMoonPosition(int day, int month, int year, int hour, int minutes, double seconds)
{
// save global variables:
Mday=day;
Mmonth=month;
Myear=year;
Mhour=hour;
Mminutes=minutes;
Mseconds=seconds;
// terrestial times:
MHour = hms2H(hour,minutes,seconds); // UT
JD = JulianDate(day, month, year, hour, minutes, seconds);
// sidereal times:
setNutationParameters();
setSiderealTimes();
// POSITION MOON IN GCRS (=ECLIPTIC) COORDINATES:
getX2000(JD, Xecl);
// Precession:
setPrecessionMatrix();
MatVecProd(XeclOD, Mprec, Xecl);
// Ecliptic to Equatorial:
setMecl2eqt();
MatVecProd(Xeqt, Mecl2eqt, XeclOD);
Cart2Polar3D(Xeqt);
RA_hr = Vector_Hour;
RA_deg = Vector_Alpha;
Dec_deg = Vector_Delta;
Reqt = Vector_Length;
// Topocentric:
setXgeotopo();
Cart2Polar3D(Xtopo);
topoRA_hr = Vector_Hour;
topoRA_deg = Vector_Alpha;
topoDec_deg = Vector_Delta;
Rtopo = Vector_Length;
// Azimuth and Altitude:
topo2AzAlt(topoRA_hr, topoDec_deg, LMST, Mlatitude);
setRefraction();
}
//-------------------------------------------------------------------------------
// show all results:
void showMoonData(void)
{
printf("-------------------------------------------------------------------------------\n");
printf("\nLocation: longitude = %.6lf, latitude = %.6lf\n", Mlongitude, Mlatitude);
printf("date = "); printdate(Mday, Mmonth, Myear); printf("\n");
printf("time (UT) = "); printtime(MHour); printf("\n");
printf("Julian Date = %.4lf\n",JD );
printf("\nObliquity of the Ecliptic:\n");
printf("Epsilon True = %015.12lf deg = ", EpsilonTrue); printdeg(EpsilonTrue);
printf("\nEpsilon Mean = %015.12lf deg = ", EpsilonMean); printdeg(EpsilonMean);
printf("\n\nSidereal Times:");
printf("\nGMST = %15.12lf = ", GMST); printtime(GMST);
printf("\nGAST = %15.12lf = ", GAST); printtime(GAST);
printf("\nLMST = %15.12lf = ", LMST); printtime(LMST);
printf("\nLAST = %15.12lf = ", LAST); printtime(LAST);
printf("\n\nEcliptic Coordinate System:\n");
printf("Xecl = "); printvector3D(Xecl);
Cart2Polar3D(Xecl);
printPolar3D();
printf("\n\nPrecession Matrix:\n");
printMatrix33(Mprec);
printf("\nEcliptic Coordinate System of date:\n");
printf("XeclOD = "); printvector3D(XeclOD);
Cart2Polar3D(XeclOD);
printPolar3D();
printf("\n\nRotation Ecliptic to Equatorial system:\n");
printMatrix33(Mecl2eqt);
printf("\nEquatorial Coordinate System:\n");
printf("Xeqt = "); printvector3D(Xeqt);
Cart2Polar3D(Xeqt);
printPolar3D();
printf("\n\nObserver position:\n");
printf("Xgeo = "); printvector3D(Xgeo);
Cart2Polar3D(Xgeo);
printPolar3D();
printf("\n\nApparent topocentric position of the Moon:\n");
printf("Xtopo = "); printvector3D(Xtopo);
Cart2Polar3D(Xtopo);
printPolar3D();
printf("\n\nAzimuth = %12.8lf = ",Azimuth); printdeg(Azimuth);
printf("\nAltitude = %12.8lf = ",AltitudeRaw); printdeg(AltitudeRaw);
printf("\n\n");
showRefraction();
}
// show main results:
void showShort(void)
{
printf("-------------------------------------------------------------------------------\n");
printf("Date = "); printdate(Mday, Mmonth, Myear);
printf(" Time = %2d:%02d:%05.2lf (UT)\n",Mhour,Mminutes, Mseconds);
printf(" Julian Date = %.6lf\n",JD);
printf(" Apparent Sidereal Time: "); printtime(LMST); printf("\n");
printf("Position Moon (Equatorial, of Date):\n");
printf(" Geocentric: RA = "); printtime(RA_hr); printf(" decl = "); printdeg(Dec_deg);
printf(" R_eqt = %.0lf km\n",Reqt);
printf(" Topocentric: RA = "); printtime(topoRA_hr); printf(" decl = "); printdeg(topoDec_deg);
printf(" R_topo = %.0lf km\n",Rtopo);
printf(" Altitude = %.2lf%c Azimuth = %.2lf%c\n",Altitude, 248, Azimuth, 248); // 248 = graden-symbool
printf("-------------------------------------------------------------------------------\n\n");
}
//-------------------------------------------------------------------------------
void demoMoon(void)
{
// Yuk table ELP/MPP02 series approximation -> MoonPosition.cpp:
//showDemoTable();
// set Moon input:
setObserver(5.129, 52.086, 101.325, 10.0); // Utrecht, mean pressure at sea level, 10 degrees Celsius
//getMoonPosition(9,12,2023, 9,20,0.0); // 9 dec 2023, 9:20:00 UT
//showMoonData();
//showShort();
getMoonPosition(22,12,2023, 0,0,0); // 22-12-2023, 00:00:00 uur UT (= 1:00:00 uur CET)
showShort();
printf("Ondergang:\n");
getMoonPosition(22,12,2023, 2,30,0); // ondergang
showShort();
printf("Opkomst:\n");
getMoonPosition(22,12,2023, 12,31,0); // opkomst
showShort();
printf("Doorgang:\n");
getMoonPosition(22,12,2023, 20,02,0); // doorgang
showShort();
}
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
#pragma argsused
int main(int argc, char* argv[])
{
demoMoon();
return 0;
}
//-------------------------------------------------------------------------------
Moon.h
Code: Selecteer alles
#ifndef MOONHDEF
#define MOONHDEF
//-------------------------------------------------------------------------------
// CONSTANTS:
//-------------------------------------------------------------------------------
const double PI = 3.1415926535897932384626433832795;
double TPI = 2.0*PI;
const double D2R = PI/180.0; // degrees => radians
const double R2D = 180.0/PI; // radians => degrees
const double SEC2RAD = PI/(3600.0*180.0);
//-------------------------------------------------------------------------------
// VARIABLES:
//-------------------------------------------------------------------------------
// input variables:
int Mday, Mmonth, Myear, Mhour, Mminutes;
double Mseconds, MHour, JD;
double Mlongitude, Mlatitude, Mpressure, Mtemperature;
// Nutation parameters: [all in degrees]
// Mean and true obliquity of the ecliptic:
// Epsilon_true = Epsilon_mean + Delta_Epsilon
// DeltaPsi = nutation in longitude
// EqEq = Equation of the Equinoxes
double EpsilonMean, EpsilonTrue, DeltaEpsilon, DeltaPsi, EqEq;
// Sidereal Times:
double GMST, GAST, LMST, LAST;
// GCRS ecliptic position of the Moon:
double Xecl[3];
// Precession matrix:
double Mprec[3][3];
// GCRS ecliptic position of the Moon of Date:
double XeclOD[3];
// Rotation from ecliptic to equatorial coordinate system:
double Mecl2eqt[3][3];
// Equatorial position of the Moon:
double Xeqt[3];
double RA_hr, RA_deg, Dec_deg, Reqt;
// Observer position on Earth (in Equatorial coordinates):
double Xgeo[3];
// Topocentric position of the Moon:
double Xtopo[3];
// apparent topocentric RA and decl:
double topoRA_hr, topoRA_deg, topoDec_deg, Rtopo;
// Azimuth and Altitude without Earth atmosphere:
double Azimuth, AltitudeRaw;
// Including atmospheric refraction:
double DeltaAltitude, Altitude;
// polar values for a 3D vector:
double Vector_Length, Vector_Delta, Vector_Alpha, Vector_Hour;
//-------------------------------------------------------------------------------
#endif
AstroLibrary.cpp
Code: Selecteer alles
#include <math.h>
//---------------------------------------------------------------------------
// MATH FUNCTIONS:
//---------------------------------------------------------------------------
// truncate x = drop fractional part of x (trunc(-1.2) = -1):
double trunc(double x)
{
int sign=1;
if(x<0.0){
x=-x;
sign=-1;
}
int n = floor(x);
return((double)(sign*n));
}
// get fractional part of x (frac(-1.2) = -0.2):
double frac(double x)
{
double sign=1.0;
if(x<0){
sign=-1.0;
x = -x;
}
return(sign*(x-floor(x)));
}
// for all m>0: for all x: upmod(x, m) is in interval [0, m>
double upmod(double x, double m)
{
double r = fmod(x, m);
if(r < 0.0)
return(r + m);
else
return(r);
}
//---------------------------------------------------------------------------
// DATA OUTPUT:
//---------------------------------------------------------------------------
// print time as h:mm:ss.sss
void printtime(double T){
int h = (int)T;
T=(T-h)*60.0;
int m = (int)T;
T=(T-m)*60;
printf("%2d:%02d:%05.2lf",h,m,T);
}
// print date as dd-mm-year
void printdate(int day, int month, int year)
{
if(year < 0)
printf("%2d-%02d-(%d)", day, month, year);
else
printf("%2d-%02d-%d", day, month, year);
}
// print D degrees as d* mm' ss.sss"
void printdeg(double D){
int sign=0;
if(D<0.0){
sign=1;
D=-D;
}
int d = (int)D;
D=(D-d)*60.0;
int m = (int)D;
D=(D-m)*60;
if(sign==1)
printf("-");
printf("%d%c%02d\'%05.2lf\"",d,248,m,D); // 248 = degree symbol
}
// print vector double vec[3]:
void printvector3D(double *vec)
{
printf("[ %.4lf, %.4lf, %.4lf ]\n",vec[0], vec[1], vec[2]);
}
// print results of Cart2Polar3D(double *vec):
void printPolar3D(void)
{
printf("R = %.8lf\n", Vector_Length);
printf("D = %15.10lf = ", Vector_Delta); printdeg(Vector_Delta);
printf("\nA = %15.10lf = ", Vector_Alpha); printdeg(Vector_Alpha);
printf("\nH = %15.10lf = ", Vector_Hour); printtime(Vector_Hour);
}
// print matrix double M[3][3]:
void printMatrix33(double M[3][3])
{
printf("[%18.12lf %18.12lf %18.12lf ]\n",M[0][0],M[0][1],M[0][2]);
printf("[%18.12lf %18.12lf %18.12lf ]\n",M[1][0],M[1][1],M[1][2]);
printf("[%18.12lf %18.12lf %18.12lf ]\n",M[2][0],M[2][1],M[2][2]);
}
//---------------------------------------------------------------------------
// CONVERSIONS:
//---------------------------------------------------------------------------
// time in h:m:s to hour H (as double value)
double hms2H(int h, int m, double s)
{
return((double)h + (double)m/60.0 + s/3600.0);
}
// Cartesian vector vec[3] to polar, equatorial map style:
void Cart2Polar3D(double *vec)
{
Vector_Length = sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]); // = R
Vector_Delta = asin(vec[2]/Vector_Length)*R2D; // latitude [-90, 90]
Vector_Alpha = atan2(vec[1],vec[0])*R2D; // longitude in degrees <-180, 180]
Vector_Hour = upmod(Vector_Alpha, 360.0)/15.0; // longitude in hours [0, 24>
}
// Matrix vector product: Vout[3] = M[3][3] * Vin[3]
// in: matrix M[3][3]
// vector Vin[3]
// out: vector Vout[3]
void MatVecProd(double *Vout, double M[3][3], double *Vin)
{
Vout[0] = M[0][0]*Vin[0] + M[0][1]*Vin[1] + M[0][2]*Vin[2];
Vout[1] = M[1][0]*Vin[0] + M[1][1]*Vin[1] + M[1][2]*Vin[2];
Vout[2] = M[2][0]*Vin[0] + M[2][1]*Vin[1] + M[2][2]*Vin[2];
}
//---------------------------------------------------------------------------
// JULIAN DATE:
//---------------------------------------------------------------------------
// calculate Julian Date for date D-M-Y time h:m:s
double JulianDate(int D, int M, int Y, double h, double m, double s)
{
double day=(double)D + (h+m/60.0+s/3600.0)/24.0;
double year, month, tmp, a, b, c;
if(M<=2){
year = Y-1;
month = M+12;
}
else{
year = Y;
month = M;
}
// 15-10-1582 = start Gregorian calendar:
if((Y<1582) || (Y==1582 && M<10) || (Y==1582 && M==10 && D<15)) // before Gregorian cal:
a = 0;
else{ // Gregorian cal:
tmp = trunc(year/100.0);
a = 2 - tmp + floor(tmp/4.0);
}
if(year<0)
b = trunc(365.25*year-0.75);
else
b = trunc(365.25*year);
c = trunc(30.6001 * (month+1));
return(a + b + c + day + 1720994.5);
}
MoonPosition.cpp
Code: Selecteer alles
//---------------------------------------------------------------------------
// Lunar Ephemeris
//---------------------------------------------------------------------------
// Source:
// https://ytliu0.github.io/starCharts/docs/star_charts.pdf
// chapter 17
// Based on truncated ELP/MPP02 series using the following parameters:
// corr = 1 = Parameters fitted to JPL ephemeris DE405/DE406
// AthU = 1.454441e-04 radians = 30.000000 arcsecs
// AthV = 1.454441e-04 radians = 30.000000 arcsecs
// AthR = 100.000000 km
// tau = 50.000000
// Total number of terms in the original series: 35901
// Total number of terms in the truncated series: 42
// Statistical evaluation:
// Estimated error of the truncated series by comparing it to the original series
// at 10000 randomly chosen T between -50.000000 and 10.000000:
// Maximum deviation found in ecliptic longitude: 213.496798 arcsecs
// Root mean square deviation in ecliptic longitude: 47.665830 arcsecs
// Maximum deviation found in ecliptic latitude: 158.214315 arcsecs
// Root mean square deviation in ecliptic latitude: 34.778046 arcsecs
// Maximum deviation found in geocentric distance: 336.995129 km
// Root mean square deviation in geocentric distance: 85.856277 km
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
// standard C libs:
#include <stdio.h>
#include <math.h>
// local libs:
#include "Moon.h"
//---------------------------------------------------------------------------
// constants:
// extern const double PI;
// extern const double TPI;
//---------------------------------------------------------------------------
// Arguments for the ELP/MPP02 series:
double W1, D, F, L, Lp;
// coefficients for the ELP/MPP02 main series:
int n_main_long=19;
int i_main_long[19][4]={
{0,2,0,0},{0,-2,1,0},{0,0,1,0},{0,2,1,0},{0,0,2,0},{0,0,3,0},{0,0,-1,1},{0,0,0,1},
{0,0,1,1},{1,0,0,0},{2,0,-1,-1},{2,0,0,-1},{2,0,-2,0},{2,0,-1,0},{2,-2,0,0},{2,0,0,0},
{2,0,1,0},{4,0,-2,0},{4,0,-1,0}};
double A_main_long[19]={
-1.99547249801829270800e-03, 1.91662856547817478400e-04, 1.09759808626191612800e-01, -2.18649080896658580300e-04,
3.72834182278018950000e-03, 1.75133187917569992000e-04, -7.14234214377495904400e-04, -3.23088338567132336800e-03,
-5.30290959258585535000e-04, -6.05959497615181906000e-04, 9.95981590587821435600e-04, 7.98626872936074750500e-04,
1.02613505456892778600e-03, 2.22356802322449947900e-02, 2.67505934260747291300e-04, 1.14896669562630695400e-02,
9.30629894589579480800e-04, 1.49189651057941429400e-04, 1.86313093175213899100e-04};
int n_main_lat=10;
int i_main_lat[10][4]={
{0,1,0,0},{0,-1,1,0},{0,1,1,0},{0,-1,2,0},{0,1,2,0},{2,-1,-1,0},{2,1,-1,0},{2,-1,0,0},{2,1,0,0},{2,-1,1,0}};
double A_main_lat[10]={
8.95026190647627839500e-02, 4.84665164815963230000e-03, 4.89742857492255542000e-03, 1.53975229251259538900e-04,
3.00157605100467080500e-04, 8.07574106044906285200e-04, 9.67124565462755482600e-04, 3.02355257133907857600e-03,
5.68495899772767385500e-04, 1.61720283648969048700e-04};
int n_main_dist=12;
int i_main_dist[12][4]={
{0,0,0,0},{0,0,1,0},{0,0,2,0},{0,0,-1,1},{0,0,1,1},{1,0,0,0},
{2,0,-1,-1},{2,0,0,-1},{2,0,-2,0},{2,0,-1,0},{2,0,0,0},{2,0,1,0}};
double A_main_dist[12]={
3.85000529032284044600e+05, -2.09053549352050904400e+04, -5.69925115335094687900e+02, -1.29620222172050574700e+02,
1.04755292674270279460e+02, 1.08742701272463051500e+02, -1.52137809425890736800e+02, -2.04586116460873086000e+02,
2.46158474853557947900e+02, -3.69911089639807642000e+03, -2.95596755797241485200e+03, -1.70733077124770574100e+02};
double main_long, main_lat, main_dist;
double pert_longT1;
//---------------------------------------------------------------------------
// restrict x to [-pi,pi>
double mod2pi(double x)
{
return(x - TPI*floor((x + PI)/TPI));
}
// Compute the lunar and planetary arguments used in the ELP/MPP02 series:
// set Elp arguments: W1, D, F, L, Lp
void compute_Elp_arguments(double T)
{
double T2 = T*T;
double T3 = T*T2;
double T4 = T2*T2;
W1 = -2.4728412163487064 + mod2pi(8399.6847300719273*T) + mod2pi(-3.3189520425500949e-05*T2)
+ mod2pi(3.1102494491060616e-08*T3) + mod2pi(-2.0328237648922845e-10*T4);
double W2 = 1.4547895404440774 + mod2pi(70.993305448479248*T) + mod2pi(-0.00018548192818782712*T2)
+ mod2pi(-2.1961637966359414e-07*T3) + mod2pi(1.0327016221314225e-09*T4);
double W3 = 2.1824388474237688 + mod2pi(-33.781427419672326*T) + mod2pi(3.0816644950982666e-05*T2)
+ mod2pi(3.6447710769397583e-08*T3) + mod2pi(-1.738541860458796e-10*T4);
double Ea = 1.7534699452640699 + mod2pi(628.30758508103168*T) + mod2pi(-9.7932363584126268e-08*T2)
+ mod2pi(4.3633231299858238e-11*T3) + mod2pi(7.2722052166430391e-13*T4);
double pomp = 1.7965956331206001 + mod2pi(0.0056298669711442699*T) + mod2pi(2.5659491293243858e-06*T2)
+ mod2pi(-5.7275888286280579e-10*T3) + mod2pi(5.5166948773454105e-11*T4);
// Mean longitude of the Moon:
W1 = mod2pi(W1);
// Arguments of Delaunay:
D = mod2pi(W1 - Ea + PI);
F = mod2pi(W1 - W3);
L = mod2pi(W1 - W2);
Lp = mod2pi(Ea - pomp);
}
//---------------------------------------------------------------------------
// Sum the ELP/MPP02 main and perturbation series:
void compute_Elp_sums(void)
{
int i;
double phase;
// main sums (long and lat: sin-sums; dist: cos-sum):
main_long=0.0;
for(i=0;i<n_main_long;i++){
phase=i_main_long[i][0]*D + i_main_long[i][1]*F + i_main_long[i][2]*L + i_main_long[i][3]*Lp;
main_long += A_main_long[i]*sin(phase);
}
main_lat=0.0;
for(i=0;i<n_main_lat;i++){
phase=i_main_lat[i][0]*D + i_main_lat[i][1]*F + i_main_lat[i][2]*L + i_main_lat[i][3]*Lp;
main_lat += A_main_lat[i]*sin(phase);
}
main_dist=0.0;
for(i=0;i<n_main_dist;i++){
phase=i_main_dist[i][0]*D + i_main_dist[i][1]*F + i_main_dist[i][2]*L + i_main_dist[i][3]*Lp;
main_dist += A_main_dist[i]*cos(phase);
}
// perturbations sums: pert_longT1 is the only non-zero for our approximation:
pert_longT1 = 8.1293558048447e-06 * sin(Lp);
}
//---------------------------------------------------------------------------
// Calculate the Moon's geocentric X,Y,Z coordinates with respect to
// J2000.0 mean ecliptic and equinox.
void getX2000(double JD, double &X, double &Y, double &Z)
{
// T is the TDB Julian century from J2000.0 = (JD - 2451545)/36525
double T=(JD - 2451545.0)/36525.0;
double T2 = T*T;
double T3 = T*T2;
double T4 = T2*T2;
double T5 = T2*T3;
compute_Elp_arguments(T);
compute_Elp_sums();
// Moon's longitude, latitude and distance
double longM = W1 + main_long + mod2pi(pert_longT1*T);
double latM = main_lat;
// ra0 = a0(DE405)/a0(ELP) = 384747.961370173/384747.980674318 = 0.99999994982652041.
double r = 0.99999994982652041*main_dist;
double x0 = r*cos(longM)*cos(latM);
double y0 = r*sin(longM)*cos(latM);
double z0 = r*sin(latM);
// Precession matrix
double P = 0.10180391e-4*T + 0.47020439e-6*T2 - 0.5417367e-9*T3 - 0.2507948e-11*T4 + 0.463486e-14*T5;
double Q = -0.113469002e-3*T + 0.12372674e-6*T2 + 0.12654170e-8*T3 - 0.1371808e-11*T4 - 0.320334e-14*T5;
double sq = sqrt(1 - P*P - Q*Q);
double p11 = 1 - 2*P*P;
double p12 = 2*P*Q;
double p13 = 2*P*sq;
double p21 = 2*P*Q;
double p22 = 1-2*Q*Q;
double p23 = -2*Q*sq;
double p31 = -2*P*sq;
double p32 = 2*Q*sq;
double p33 = 1 - 2*P*P - 2*Q*Q;
bool show=false;
if(show){
printf("precession matrix: -------------------\n");
printf("%16.8lf%16.8lf%16.8lf\n",p11,p12,p13);
printf("%16.8lf%16.8lf%16.8lf\n",p21,p22,p23);
printf("%16.8lf%16.8lf%16.8lf\n",p31,p32,p33);
printf("--------------------------------------\n");
}
// Finally, calculate the return values X, Y and Z:
// (= the components of the position vector; wrt J2000.0 mean ecliptic and equinox)
X = p11*x0 + p12*y0 + p13*z0;
Y = p21*x0 + p22*y0 + p23*z0;
Z = p31*x0 + p32*y0 + p33*z0;
}
// return the getX2000 results as double vector v[3] = [X, Y, Z] :
void getX2000(double JD, double *v)
{
double x, y, z;
getX2000(JD, x, y, z);
v[0]=x;
v[1]=y;
v[2]=z;
}
//---------------------------------------------------------------------------
// Demo:
//---------------------------------------------------------------------------
// Demonstration output for our ELP/MPP02 series approximation
// Compare with the high precision version on page 6 in
// https://github.com/ytliu0/ElpMpp02/blob/master/docs/ElpMpp02.pdf
void showTableHeader(void)
{
printf("\n");
printf("Lunar geocentric rectangular coordinates (in km), fitted to DE405/DE406\n");
printf("-------------------------------------------------------------------------------------------\n");
printf("Jul. Date Date Time X Y Z R\n");
printf("-------------------------------------------------------------------------------------------\n");
}
void showMoonPosition(int D, int M, int Y, int h, int m, double s)
{
double JD=JulianDate(D,M,Y,h,m,s); // calculate Julian Date
double x,y,z,R; // lunar coordinates and distance
getX2000(JD,x,y,z); // calculate lunar coordinates on JD
R=sqrt(x*x+y*y+z*z);
printf("%.2lf %5d-%02d-%02d %02d:%02d:%02.0lf %10.0lf %10.0lf %10.0lf %10.0lf\n",JD,Y,M,D,h,m,s,x,y,z,R);
}
// Table at bottom of pag. 6:
void showDemoTable(void)
{
showTableHeader();
showMoonPosition(13,6,2192,4,4,48);
showMoonPosition(7,12,1490,19,55,12);
showMoonPosition(16,6,789,11,45,36);
showMoonPosition(25,12,87,3,36,0);
showMoonPosition(3,7,-614,19,26,24);
printf("\n");
}
//---------------------------------------------------------------------------
// <<< end of Yuk Tung Liu's Lunar Ephemeris
//---------------------------------------------------------------------------
Output:
Hier de ouput van bovenstaande:
Code: Selecteer alles
-------------------------------------------------------------------------------
Date = 22-12-2023 Time = 0:00:00.00 (UT)
Julian Date = 2460300.500000
Apparent Sidereal Time: 6:21:42.03
Position Moon (Equatorial, of Date):
Geocentric: RA = 1:44:39.63 decl = 11°22'29.58" R_eqt = 376070 km
Topocentric: RA = 1:42:22.13 decl = 10°39'49.11" R_topo = 373766 km
Altitude = 20.79° Azimuth = 260.55°
-------------------------------------------------------------------------------
Ondergang:
-------------------------------------------------------------------------------
Date = 22-12-2023 Time = 2:30:00.00 (UT)
Julian Date = 2460300.604167
Apparent Sidereal Time: 8:52:06.68
Position Moon (Equatorial, of Date):
Geocentric: RA = 1:49:54.67 decl = 11°59'49.83" R_eqt = 376332 km
Topocentric: RA = 1:47:33.74 decl = 11°13'02.63" R_topo = 376375 km
Altitude = -0.19° Azimuth = 289.55°
-------------------------------------------------------------------------------
Opkomst:
-------------------------------------------------------------------------------
Date = 22-12-2023 Time = 12:31:00.00 (UT)
Julian Date = 2460301.021528
Apparent Sidereal Time: 18:54:45.40
Position Moon (Equatorial, of Date):
Geocentric: RA = 2:11:08.19 decl = 14°25'02.75" R_eqt = 377406 km
Topocentric: RA = 2:13:27.32 decl = 13°37'54.46" R_topo = 377457 km
Altitude = -0.24° Azimuth = 66.23°
-------------------------------------------------------------------------------
Doorgang:
-------------------------------------------------------------------------------
Date = 22-12-2023 Time = 20:02:00.00 (UT)
Julian Date = 2460301.334722
Apparent Sidereal Time: 2:26:59.49
Position Moon (Equatorial, of Date):
Geocentric: RA = 2:27:17.26 decl = 16°08'50.66" R_eqt = 378240 km
Topocentric: RA = 2:27:17.45 decl = 15°34'34.66" R_topo = 373093 km
Altitude = 53.50° Azimuth = 179.88°
-------------------------------------------------------------------------------
Ter vergelijk:
Sterrengids 2023 (22-12-2023):
0 uur UT: Geocentrische RA = 1 uur 44.7 min, decl = 11°23', R = 376064 km
Ondergang: 2:30, Az = 289.6 graden
Opkomst: 12:31, Az = 66.3 graden
Doorgang: 19:14, Altitude = 53.5 graden
En via deze site: https://ytliu0.github.io/starCharts/:
Heel wat informatie, stel gerust vragen over onduidelijkheden in bovenstaande.
Ik hoop dat het al voldoende nauwkeurig is, en dat de code eenvoudig is om te zetten naar Delphi.
PS: hier nog de tussenresultaten van 22-12-2023 0 uur UT, verkregen via showMoonData();
Code: Selecteer alles
Location: longitude = 5.129000, latitude = 52.086000
date = 22-12-2023
time (UT) = 0:00:00.00
Julian Date = 2460300.5000
Obliquity of the Ecliptic:
Epsilon True = 23.438381877115 deg = 23°26'18.17"
Epsilon Mean = 23.436160738023 deg = 23°26'10.18"
Sidereal Times:
GMST = 6.019742784719 = 6:01:11.07
GAST = 6.018139356346 = 6:01:05.30
LMST = 6.361676118052 = 6:21:42.03
LAST = 6.360072689679 = 6:21:36.26
Ecliptic Coordinate System:
Xecl = [ 331941.9979, 176725.9009, 3380.9746 ]
R = 376070.42562091
D = 0.5151114507 = 0°30'54.40"
A = 28.0308346588 = 28°01'51.00"
H = 1.8687223106 = 1:52:07.40
Precession Matrix:
[ 0.999982918492 -0.005844886866 -0.000004614425 ]
[ 0.005844886606 0.999982917024 -0.000054420371 ]
[ 0.000004932427 0.000054392470 0.999999998509 ]
Ecliptic Coordinate System of date:
XeclOD = [ 330903.3694, 178662.8612, 3392.2244 ]
R = 376070.42562091
D = 0.5168254756 = 0°31'00.57"
A = 28.3657003303 = 28°21'56.52"
H = 1.8910466887 = 1:53:27.77
Rotation Ecliptic to Equatorial system:
[ 1.000000000000 0.000000000000 0.000000000000 ]
[ 0.000000000000 0.917503702391 -0.397727238317 ]
[ 0.000000000000 0.397727238317 0.917503702391 ]
Equatorial Coordinate System:
Xeqt = [ 330903.3694, 162574.6566, 74171.4648 ]
R = 376070.42562091
D = 11.3748828408 = 11°22'29.58"
A = 26.1651418862 = 26°09'54.51"
H = 1.7443427924 = 1:44:39.63
Observer position:
Xgeo = [ -371.3181, 3909.8227, 5008.6886 ]
R = 6364.86861879
D = 51.8992790853 = 51°53'57.40"
A = 95.4251417708 = 95°25'30.51"
H = 6.3616761181 = 6:21:42.03
Apparent topocentric position of the Moon:
Xtopo = [ 331274.6875, 158664.8339, 69162.7762 ]
R = 373765.88621265
D = 10.6636417919 = 10°39'49.11"
A = 25.5922063531 = 25°35'31.94"
H = 1.7061470902 = 1:42:22.13
Azimuth = 260.54716575 = 260°32'49.80"
Altitude = 20.74288133 = 20°44'34.37"
Atmospheric Refraction:
Delta Altitude = 0.0440802745 = 0°02'38.69"
Observed Altitude = 20.7869616015 = 20°47'13.06"
Re: Maan Positie
Heel erg bedankt !!!
Ik durf het bijna niet te vragen, maar is het ook mogelijk om het verlichte percentage van de maan te berekenen, en in welke fase de maan is (Nieuwe maan ..... Volle maan) ?
Ik durf het bijna niet te vragen, maar is het ook mogelijk om het verlichte percentage van de maan te berekenen, en in welke fase de maan is (Nieuwe maan ..... Volle maan) ?
Re: Maan Positie
Zie Yuk Tung Liu - Calculations in Star Charts:
- paragraaf 17.3: elongatie, fasehoek etc, zie ook figuur 8 op pagina 41
Hierbij hebben we ook nodig:
- paragraaf 16.1: de heliocentrische positie van de planeten
Voor dit laatste verwijst Yuk naar https://ssd.jpl.nasa.gov/planets/approx_pos.html.
Bovenstaande in twee functies gevat:
getPlanetPosition(): het nasa algoritme voor 3000 BC .. 3000 AD; wij hebben alleen planet=2=Aarde nodig,
getAppearanceMoon(): let wel op: hier is:
- de invoer van de heliocentrische positie van de Aarde in AU's (= de nasa output),
- de invoer van de geocentrische positie van de Maan in km's (= de Yuk output)
Aanroep (opnieuw met de gegevens uit mijn eerdere berichten):
Uitvoer:
En voor 22-12-2023 0:00:00 UT (zoals het plaatje in mijn vorige post):
Elongatie > 0 = maan ten westen, elongatie < 0 = maan ten oosten van de zon, dus
180° = volle maan
90° = laatste kwartier
0° = nieuwe maan
-90° = eerste kwartier
Noot: je kan voor hogere nauwkeurigheid ook tabel 1 in de nasa-link kiezen, maar die is slechts bruikbaar
tussen 1800 AD en 2050 AD.
- paragraaf 17.3: elongatie, fasehoek etc, zie ook figuur 8 op pagina 41
Hierbij hebben we ook nodig:
- paragraaf 16.1: de heliocentrische positie van de planeten
Voor dit laatste verwijst Yuk naar https://ssd.jpl.nasa.gov/planets/approx_pos.html.
Bovenstaande in twee functies gevat:
getPlanetPosition(): het nasa algoritme voor 3000 BC .. 3000 AD; wij hebben alleen planet=2=Aarde nodig,
getAppearanceMoon(): let wel op: hier is:
- de invoer van de heliocentrische positie van de Aarde in AU's (= de nasa output),
- de invoer van de geocentrische positie van de Maan in km's (= de Yuk output)
Code: Selecteer alles
//------------------------------------------------------------------------------------------------------
// used libraries:
//------------------------------------------------------------------------------------------------------
#include <math.h>
#include "moon.h"
//------------------------------------------------------------------------------------------------------
// used external functions:
//------------------------------------------------------------------------------------------------------
extern void setPrecessionMatrix(void);
extern void MatVecProd(double *Vout, double M[3][3], double *Vin);
//------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------
// APPROXIMATE POSITIONS OF THE PLANETS
// https://ssd.jpl.nasa.gov/planets/approx_pos.html
// in: planet: 0=Mercury .. 7=Neptune
// JD: Julian Date
// out: PeclOD[3]: planet's ecliptic of date position [AU]
void getPlanetPosition(int planet, double JD, double PeclOD[3])
{
// Table 2a: Keplerian elements and their rates, with respect to the mean ecliptic and equinox of J2000,
// valid for the time-interval 3000 BC .. 3000 AD.
//------------------------------------------------------------------------------------------------------
// a e I L long.peri. long.node.
// au, au/Cy rad, rad/Cy deg, deg/Cy deg, deg/Cy deg, deg/Cy deg, deg/Cy
//------------------------------------------------------------------------------------------------------
double KepEl[8][2][6]={
// TABLE 2A:
{{ 0.38709843, 0.20563661, 7.00559432, 252.25166724, 77.45771895, 48.33961819}, // Mercury
{ 0.00000000, 0.00002123, -0.00590158, 149472.67486623, 0.15940013, -0.12214182}},
{{ 0.72332102, 0.00676399, 3.39777545, 181.97970850, 131.76755713, 76.67261496}, // Venus
{ -0.00000026, -0.00005107, 0.00043494, 58517.81560260, 0.05679648, -0.27274174}},
{{ 1.00000018, 0.01673163, -0.00054346, 100.46691572, 102.93005885, -5.11260389}, // EM Bary
{ -0.00000003, -0.00003661, -0.01337178, 35999.37306329, 0.31795260, -0.24123856}},
{{ 1.52371243, 0.09336511, 1.85181869, -4.56813164, -23.91744784, 49.71320984}, // Mars
{ 0.00000097, 0.00009149, -0.00724757, 19140.29934243, 0.45223625, -0.26852431}},
{{ 5.20248019, 0.04853590, 1.29861416, 34.33479152, 14.27495244, 100.29282654}, // Jupiter
{ -0.00002864, 0.00018026, -0.00322699, 3034.90371757, 0.18199196, 0.13024619}},
{{ 9.54149883, 0.05550825, 2.49424102, 50.07571329, 92.86136063, 113.63998702}, // Saturn
{ -0.00003065, -0.00032044, 0.00451969, 1222.11494724, 0.54179478, -0.25015002}},
{{ 19.18797948, 0.04685740, 0.77298127, 314.20276625, 172.43404441, 73.96250215}, // Uranus
{ -0.00020455, -0.00001550, -0.00180155, 428.49512595, 0.09266985, 0.05739699}},
{{ 30.06952752, 0.00895439, 1.77005520, 304.22289287, 46.68158724, 131.78635853}, // Neptune
{ 0.00006447, 0.00000818, 0.00022400, 218.46515314, 0.01009938, -0.00606302}}
};
//------------------------------------------------------------------------------------------------------
// Additional terms which must be added to the computation of M for Jupiter through Neptune,
// 3000 BC to 3000 AD, as described in the related document:
// M += b*T^2 + c*cos(f*T) + s*sin(f*T)
//-----------------------------------------------------------------------
// b c s f
// deg/cty^2 deg deg deg/cty
//-----------------------------------------------------------------------
double K2b[4][4]={
{ -0.00012452, 0.06064060, -0.35635438, 38.35125000}, // Jupiter
{ 0.00025899, -0.13434469, 0.87320147, 38.35125000}, // Saturn
{ 0.00058331, -0.97731848, 0.17689245, 7.67025000}, // Uranus
{ -0.00041348, 0.68346318, -0.10162547, 7.67025000} // Neptune
};
//-----------------------------------------------------------------------
// step 1: compute planet parameters:
double T = (JD - 2451545.0) / 36525.0;
double a = KepEl[planet][0][0] + T*KepEl[planet][1][0];
double e = KepEl[planet][0][1] + T*KepEl[planet][1][1];
double I = KepEl[planet][0][2] + T*KepEl[planet][1][2];
double L = KepEl[planet][0][3] + T*KepEl[planet][1][3];
double wt = KepEl[planet][0][4] + T*KepEl[planet][1][4];
double Om = KepEl[planet][0][5] + T*KepEl[planet][1][5];
// step 2: Compute w (= argument of perihelion) and M (= the mean anomaly):
double w = wt - Om;
double M = L - wt;
if(planet>=4){
int p=planet-4;
M += K2b[p][0]*T*T + K2b[p][1]*cos(K2b[p][3]*T*D2R) + K2b[p][2]*sin(K2b[p][3]*T*D2R);
}
// step 3a: modulus the mean anomaly M:
M=fmod(M,360.0);
if(M <= -180.0) M+=360.0;
if(M > 180.0) M-=360.0;
// now -180 < M <= 180
// step 3b: solve Kepler's equation:
double estar = e*R2D;
double E = M - estar*sin(M*D2R);
double tol=1e-10;
double DE=tol+0.01;
double DM;
int iter=0;
while(fabs(DE)>tol && iter<1000){ // just for safety: max 1000 iterations to prevent endless looping
DM=M-E+estar*sin(E*D2R);
DE=DM/(1.0-e*cos(E*D2R));
E+=DE;
iter++;
}
// step 4: Compute the planet's heliocentric coordinates, aligned from the focus to the perihelion:
double xacc = a*(cos(E*D2R)-e);
double yacc = a*sqrt(1.0-e*e)*sin(E*D2R);
// step 5: Compute Pecl[3] = the planet's coordinates in the J2000 ecliptic plane,
// with the x-axis aligned toward the equinox:
double Pecl[3];
Pecl[0] = (cos(w*D2R)*cos(Om*D2R)-sin(w*D2R)*sin(Om*D2R)*cos(I*D2R)) * xacc;
Pecl[0] += (-sin(w*D2R)*cos(Om*D2R)-cos(w*D2R)*sin(Om*D2R)*cos(I*D2R)) * yacc;
Pecl[1] = (cos(w*D2R)*sin(Om*D2R)+sin(w*D2R)*cos(Om*D2R)*cos(I*D2R)) * xacc;
Pecl[1] += (-sin(w*D2R)*sin(Om*D2R)+cos(w*D2R)*cos(Om*D2R)*cos(I*D2R)) * yacc;
Pecl[2] = (sin(w*D2R)*sin(I*D2R)) * xacc;
Pecl[2] += (cos(w*D2R)*sin(I*D2R)) * yacc;
// rotate J2000 grid to of date grid:
setPrecessionMatrix();
MatVecProd(PeclOD, Mprec, Pecl);
}
//------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------
// Get the elongation, phase angle, magnitude and illumination of the Moon:
// in: vse: vector from sun to earth [AU]
// vem: vector from earth to moon [km]
// out: E = Moon's solar elongation [deg]
// PA = phase angle [deg]
// Mag = Apparent magnitude
// Illum = illumination
void getAppearanceMoon(double vse[3], double vem[3], double &E, double &PA, double &Mag, double &Illum)
{
double re=vectorlength(vse); // [AU]
double rm=vectorlength(vem); // [km]
double dm=asin(Xecl[2]/rm); // [rad]
double am=atan2(Xecl[1],Xecl[0]); // [rad]
double ds=asin(-vse[2]/re); // decl earth to sun [rad]
double as=atan2(-vse[1], -vse[0]); // ra earth to sun [rad]
double cosE=sin(ds)*sin(dm)+cos(ds)*cos(dm)*cos(as-am);
E=acos(cosE)*R2D;
if(vse[0]*Xecl[1] < vse[1]*Xecl[0]) // => East => E negative:
E*=(-1.0);
double rsun=re*AU;
double sigma=rm/rsun;
double cosP=(sigma-cosE)/sqrt(1.0-2.0*sigma*cosE+sigma*sigma);
PA=acos(cosP)*R2D;
Mag = -12.73 + 0.026*PA + 4E-9*(PA*PA*PA*PA) + 5*log(re*rm/384400.0)/log(10.0);
Illum = (1.0+cos(PA*D2R))/2.0;
}
Aanroep (opnieuw met de gegevens uit mijn eerdere berichten):
Code: Selecteer alles
setObserver(5.129, 52.086, 101.325, 10.0); // Utrecht, mean pressure at sea level, 10 degrees Celsius
getMoonPosition(9,12,2023, 9,20,0.0); // 9 dec 2023, 9:20:00 UT
double EarthPos[3];
getPlanetPosition(2, JD, EarthPos); // planet 2 = Earth
printf("Earth Position:\n");
for(int i=0;i<3;i++)
printf("EarthPos[%d] = %13.10lf AU\n",i, EarthPos[i]);
double elon, pa, mag, illum;
getAppearanceMoon(EarthPos, Xecl, elon, pa, mag, illum);
printf("\nMoon Appearance:\n");
printf("Elongation = %.2lf%c\n",elon,CHR_DEGREE);
printf("Phase angle = %.2lf%c\n",pa,CHR_DEGREE);
printf("Magnitude = %.2lf\n",mag);
printf("Illumination = %.4lf\n",illum);
Uitvoer:
Code: Selecteer alles
Earth Position:
EarthPos[0] = 0.2210223296 AU
EarthPos[1] = 0.9598364534 AU
EarthPos[2] = -0.0000105409 AU
Moon Appearance:
Elongation = 44.36°
Phase angle = 135.53°
Magnitude = -7.85
Illumination = 0.1432
En voor 22-12-2023 0:00:00 UT (zoals het plaatje in mijn vorige post):
Code: Selecteer alles
Moon Appearance:
Elongation = -118.17°
Phase angle = 61.70°
Magnitude = -11.15
Illumination = 0.7370
Elongatie > 0 = maan ten westen, elongatie < 0 = maan ten oosten van de zon, dus
180° = volle maan
90° = laatste kwartier
0° = nieuwe maan
-90° = eerste kwartier
Noot: je kan voor hogere nauwkeurigheid ook tabel 1 in de nasa-link kiezen, maar die is slechts bruikbaar
tussen 1800 AD en 2050 AD.