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
Heel erg bedankt voor deze info, maar dit is echt te hoog gegrepen voor mij.
Dit krijg ik helaas nooit omgezet in de manier van programmeren (Pascal) die ik gewend ben.
Dit krijg ik helaas nooit omgezet in de manier van programmeren (Pascal) die ik gewend ben.
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
Mocht er iemand de kennis, tijd en zin hebben om ook deze formule(s) om te zetten in een pascal achtige code, dan zou ik daar erg blij mee zijn.
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
Lukt het een beetje ?
Re: Maan Positie
Ja, er is nog hoop.
Ik begin nu overzicht te krijgen over alle tijd- en coordinaat-systemen.
Maar ik heb nog even nodig om alles werkende te krijgen en te controleren aan de hand van bekende waarden.
Ik verwacht over 1 a 2 weken meer te kunnen melden.
Ik begin nu overzicht te krijgen over alle tijd- en coordinaat-systemen.
Maar ik heb nog even nodig om alles werkende te krijgen en te controleren aan de hand van bekende waarden.
Ik verwacht over 1 a 2 weken meer te kunnen melden.