Maan Positie

Wiskunde is niet alleen een vak op school. Kom je ergens in de praktijk (bijvoorbeeld tijdens je werk) een wiskundig probleem tegen dan kun je hier om hulp vragen.
MMSoft
Nieuw lid
Nieuw lid
Berichten: 15
Lid geworden op: 06 jun 2014, 09:20

Maan Positie

Bericht door MMSoft » 31 okt 2023, 10:13

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 ?

arie
Moderator
Moderator
Berichten: 3928
Lid geworden op: 09 mei 2008, 09:19

Re: Maan Positie

Bericht door arie » 31 okt 2023, 12:17


MMSoft
Nieuw lid
Nieuw lid
Berichten: 15
Lid geworden op: 06 jun 2014, 09:20

Re: Maan Positie

Bericht door MMSoft » 31 okt 2023, 22:32

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));

arie
Moderator
Moderator
Berichten: 3928
Lid geworden op: 09 mei 2008, 09:19

Re: Maan Positie

Bericht door arie » 01 nov 2023, 09:01

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)

MMSoft
Nieuw lid
Nieuw lid
Berichten: 15
Lid geworden op: 06 jun 2014, 09:20

Re: Maan Positie

Bericht door MMSoft » 02 nov 2023, 08:10

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.

arie
Moderator
Moderator
Berichten: 3928
Lid geworden op: 09 mei 2008, 09:19

Re: Maan Positie

Bericht door arie » 04 nov 2023, 10:34

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:

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
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:

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
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?

MMSoft
Nieuw lid
Nieuw lid
Berichten: 15
Lid geworden op: 06 jun 2014, 09:20

Re: Maan Positie

Bericht door MMSoft » 05 nov 2023, 14:16

Ik heb de code omgezet naar Delphi (Pascal):

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;

//------------------------------------------------------------------------------
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 ?

arie
Moderator
Moderator
Berichten: 3928
Lid geworden op: 09 mei 2008, 09:19

Re: Maan Positie

Bericht door arie » 05 nov 2023, 17:40

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

MMSoft
Nieuw lid
Nieuw lid
Berichten: 15
Lid geworden op: 06 jun 2014, 09:20

Re: Maan Positie

Bericht door MMSoft » 05 nov 2023, 19:43

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.

arie
Moderator
Moderator
Berichten: 3928
Lid geworden op: 09 mei 2008, 09:19

Re: Maan Positie

Bericht door arie » 07 nov 2023, 15:42

We weten nu de positie van de maan ten opzichte van de aarde als functie van tijd.
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
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:

Afbeelding

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).

arie
Moderator
Moderator
Berichten: 3928
Lid geworden op: 09 mei 2008, 09:19

Re: Maan Positie

Bericht door arie » 09 dec 2023, 12:06

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

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);
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
MoonPosition.h = Yuk Tung Liu's Lunar Ephemeris

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
//---------------------------------------------------------------------------

Moon.c = hoofdprogramma

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']

MMSoft
Nieuw lid
Nieuw lid
Berichten: 15
Lid geworden op: 06 jun 2014, 09:20

Re: Maan Positie

Bericht door MMSoft » 11 dec 2023, 11:30

Alvast bedankt voor de moeite !!!

arie
Moderator
Moderator
Berichten: 3928
Lid geworden op: 09 mei 2008, 09:19

Re: Maan Positie

Bericht door arie » 22 dec 2023, 23:09

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

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/:

Afbeelding


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"

MMSoft
Nieuw lid
Nieuw lid
Berichten: 15
Lid geworden op: 06 jun 2014, 09:20

Re: Maan Positie

Bericht door MMSoft » 02 jan 2024, 19:27

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) ?

arie
Moderator
Moderator
Berichten: 3928
Lid geworden op: 09 mei 2008, 09:19

Re: Maan Positie

Bericht door arie » 04 jan 2024, 16:08

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)

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.

Plaats reactie