Maan Positie

[MMSoft]

Bedankt voor al deze info !!!

[arie]

Ik kwam nog een behoorlijk goede benadering tegen voor de positie van de Zon en de Maan, zie de code hieronder.
Je hebt alleen de Julian date (JD) nodig, het algoritme bepaalt daarmee direct de Equatoriale Coordinaten (RA = Right ascension en decl = declinatie).
De nauwkeurigheid lijkt binnen een boogminuut te liggen, ik moet dit nog wel even goed bekijken, evenals over welke eeuwen deze benadering nog bruikbare resultaten levert.

Dit is een sneller en eenvoudiger alternatief voor de berekeningen van Yuk hierboven.
We hebben dan alleen nog het gedeelte vanaf de omzetting van Equatoriale naar Topocentrische coordinaten nodig.

Code: Selecteer alles

#include <math.h>

double fraction(double x)
{
return(x-floor(x));
}

void QuickSun(double JD, double& RA, double& decl)
{
double TPI = 6.283185307179586476925;
double T   = (JD - 2451545.0) / 36525.0;
// Mean anomaly:
double MFrac = fraction(0.993133 + 99.997361*T);
double M     = TPI * MFrac;
// Ecliptic longitude:
double lambda = TPI * fraction(0.7859453 + MFrac + (6.893*sin(M)+0.072*sin(2.0*M)+6.1912*T)/1296.0);
// obliquity ecliptic:
double eps = 0.40909280420293639;	// [rad]
decl = asin(sin(eps)*sin(lambda));
RA   = atan2(cos(eps)*sin(lambda)/cos(decl), cos(lambda)/cos(decl));
}


void QuickMoon(double JD, double& RA, double& decl)
{
double TPI = 6.283185307179586476925;
double SR  = 4.8481368110953599358991E-4;
double T   = (JD - 2451545.0) / 36525.0;

// Mean elements of lunar orbit
double Lm = fraction(0.606433 + 1336.855225*T);      // mean longitude
double l  = TPI*fraction(0.374897 + 1325.552410*T);  // Moon's mean anomaly
double ls = TPI*fraction(0.993133 +   99.997361*T);  // Sun's mean anomaly
double D  = TPI*fraction(0.827361 + 1236.853086*T);  // Diff. long. Moon-Sun
double F  = TPI*fraction(0.259086 + 1342.227825*T);  // Dist. from ascending node

// Perturbations in longitude and latitude
double dL = 226.40*sin(l)      + 7.69*sin(2*l)      -  1.10*sin(l+ls)  + 1.48*sin(l-ls)
			- 6.68*sin(ls)     + 1.92*sin(l+2*D)    - 45.86*sin(l-2*D) - 2.12*sin(2*l-2*D)
			- 1.65*sin(ls-2*D) - 1.25*sin(D)	    + 23.70*sin(2*D)   - 0.55*sin(2*F-2*D)
			- 4.12*sin(2*F)    - 2.06*sin(l+ls-2*D);

double S = F + (dL+4.12*sin(2*F)+5.41*sin(ls))*SR;
double h = F - 2*D;
double N = - 5.26*sin(h)     + 0.44*sin(l+h)    - 0.31*sin(-l+h) - 0.23*sin(ls+h)
		   + 0.11*sin(-ls+h) - 0.25*sin(-2*l+F) + 0.21*sin(-l+F);

// Ecliptic longitude and latitude
double lambda = TPI * fraction(Lm + dL/12960.0); 	// [rad]
double beta   = (185.2*sin(S) + N) * SR;   			// [rad]

double eps    = 0.40909280420293639;				// [rad]

decl = asin(sin(beta)*cos(eps)+sin(lambda)*cos(beta)*sin(eps));
double x = cos(lambda)*cos(beta)/cos(decl);
double y = (-sin(beta)*sin(eps)+sin(lambda)*cos(beta)*cos(eps))/cos(decl);
RA = atan2(y, x);
}