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