astronomicalmath.cpp

一个基于BREW上的电子时钟程序,通过太阳历计算当前时间,并有闹铃和日志功能.

				+ 0.0018 * sin(DegreeToRadian(4.1 + 4013.29 * T))
				+ 0.0016 * sin(DegreeToRadian(242.2 + 18569.38 * T))
				+ 0.0012 * sin(DegreeToRadian(339.0 + 12678.71 * T))
				+ 0.0011 * sin(DegreeToRadian(276.5 + 19208.02 * T))
				+ 0.0009 * sin(DegreeToRadian(218.0 + 8586.0 * T))
				+ 0.0008 * sin(DegreeToRadian(188.0 + 14037.3 * T))
				+ 0.0008 * sin(DegreeToRadian(204.0 + 7906.7 * T))
				+ 0.0007 * sin(DegreeToRadian(140.0 + 4052.0 * T))
				+ 0.0007 * sin(DegreeToRadian(275.0 + 4853.3 * T))
				+ 0.0007 * sin(DegreeToRadian(216.0 + 278.6 * T))
				+ 0.0006 * sin(DegreeToRadian(128.0 + 1118.7 * T))
				+ 0.0005 * sin(DegreeToRadian(247.0 + 22582.7 * T))
				+ 0.0005 * sin(DegreeToRadian(181.0 + 19088.0 * T))
				+ 0.0005 * sin(DegreeToRadian(114.0 + 17450.7 * T))
				+ 0.0005 * sin(DegreeToRadian(332.0 + 5091.3 * T))
				+ 0.0004 * sin(DegreeToRadian(313.0 + 398.7 * T))
				+ 0.0004 * sin(DegreeToRadian(278.0 + 120.1 * T))
				+ 0.0004 * sin(DegreeToRadian(71.0 + 9584.7 * T))
				+ 0.0004 * sin(DegreeToRadian(20.0 + 720.0 * T))
				+ 0.0003 * sin(DegreeToRadian(83.0 + 3814.0 * T))
				+ 0.0003 * sin(DegreeToRadian(66.0 + 3494.7 * T))
				+ 0.0003 * sin(DegreeToRadian(147.0 + 18089.3 * T))
				+ 0.0003 * sin(DegreeToRadian(311.0 + 5492.0 * T))
				+ 0.0003 * sin(DegreeToRadian(161.0 + 40.7 * T))
				+ 0.0003 * sin(DegreeToRadian(280.0 + 23221.3 * T));

	lambda_m = fmod(lambda_m, 360.0);
	return (lambda_m >= 0.0) ? lambda_m : (lambda_m + 360.0);
}

// 寧偺墿堒傪寁嶼偡傞
Float64 AstronomicalMath::Moon::CalcEclipticLatitude(Float64 T) const
{
	Float64 Bm, beta_m;

	Bm = 0.0267 * sin(DegreeToRadian(234.95 + 19.341 * T))
			+ 0.0043 * sin(DegreeToRadian(322.1 + 19.36 * T))
			+ 0.0040 * sin(DegreeToRadian(119.5 + 1.33 * T))
			+ 0.0026 * sin(DegreeToRadian(55.0 + 19.34 * T))
			+ 0.0005 * sin(DegreeToRadian(307.0 + 19.4 * T));

	beta_m = 5.1282 * sin(DegreeToRadian(93.273 + 4832.0202 * T + Bm))
				+ 0.2806 * sin(DegreeToRadian(228.235 + 9604.0088 * T))
				+ 0.2777 * sin(DegreeToRadian(138.311 + 60.0316 * T))
				+ 0.1732 * sin(DegreeToRadian(142.427 + 4073.3220 * T))
				+ 0.0554 * sin(DegreeToRadian(194.01 + 8965.374 * T))
				+ 0.0463 * sin(DegreeToRadian(172.55 + 698.667 * T))
				+ 0.0326 * sin(DegreeToRadian(328.96 + 13737.362 * T))
				+ 0.0172 * sin(DegreeToRadian(3.18 + 14375.997 * T))
				+ 0.0093 * sin(DegreeToRadian(277.4 + 8845.31 * T))
				+ 0.0088 * sin(DegreeToRadian(176.7 + 4711.96 * T))
				+ 0.0082 * sin(DegreeToRadian(144.9 + 3713.33 * T))
				+ 0.0043 * sin(DegreeToRadian(307.6 + 5470.66 * T))
				+ 0.0042 * sin(DegreeToRadian(103.9 + 18509.35 * T))
				+ 0.0034 * sin(DegreeToRadian(319.9 + 4433.31 * T))
				+ 0.0025 * sin(DegreeToRadian(196.5 + 8605.38 * T))
				+ 0.0022 * sin(DegreeToRadian(331.4 + 13377.37 * T))
				+ 0.0021 * sin(DegreeToRadian(170.1 + 1058.66 * T))
				+ 0.0019 * sin(DegreeToRadian(230.7 + 9244.02 * T))
				+ 0.0018 * sin(DegreeToRadian(243.3 + 8206.68 * T))
				+ 0.0018 * sin(DegreeToRadian(270.8 + 5192.01 * T))
				+ 0.0017 * sin(DegreeToRadian(99.8 + 14496.06 * T))
				+ 0.0016 * sin(DegreeToRadian(135.7 + 420.02 * T))
				+ 0.0015 * sin(DegreeToRadian(211.1 + 9284.69 * T))
				+ 0.0015 * sin(DegreeToRadian(45.8 + 9964.00 * T))
				+ 0.0014 * sin(DegreeToRadian(219.2 + 299.96 * T))
				+ 0.0013 * sin(DegreeToRadian(95.8 + 4472.03 * T))
				+ 0.0013 * sin(DegreeToRadian(155.4 + 379.35 * T))
				+ 0.0012 * sin(DegreeToRadian(38.4 + 4812.68 * T))
				+ 0.0012 * sin(DegreeToRadian(148.2 + 4851.36 * T))
				+ 0.0011 * sin(DegreeToRadian(138.3 + 19147.99 * T))
				+ 0.0010 * sin(DegreeToRadian(18.0 + 12978.66 * T))
				+ 0.0008 * sin(DegreeToRadian(70.0 + 17870.7 * T))
				+ 0.0008 * sin(DegreeToRadian(326.0 + 9724.1 * T))
				+ 0.0007 * sin(DegreeToRadian(294.0 + 13098.7 * T))
				+ 0.0006 * sin(DegreeToRadian(224.0 + 5590.7 * T))
				+ 0.0006 * sin(DegreeToRadian(52.0 + 13617.3 * T))
				+ 0.0005 * sin(DegreeToRadian(280.0 + 8485.3 * T))
				+ 0.0005 * sin(DegreeToRadian(239.0 + 4193.4 * T))
				+ 0.0004 * sin(DegreeToRadian(311.0 + 9483.9 * T))
				+ 0.0004 * sin(DegreeToRadian(238.0 + 23281.3 * T))
				+ 0.0004 * sin(DegreeToRadian(81.0 + 10242.6 * T))
				+ 0.0004 * sin(DegreeToRadian(13.0 + 9325.4 * T))
				+ 0.0004 * sin(DegreeToRadian(147.0 + 14097.4 * T))
				+ 0.0003 * sin(DegreeToRadian(205.0 + 22642.7 * T))
				+ 0.0003 * sin(DegreeToRadian(107.0 + 18149.4 * T))
				+ 0.0003 * sin(DegreeToRadian(146.0 + 3353.3 * T))
				+ 0.0003 * sin(DegreeToRadian(234.0 + 19268.0 * T));

	return beta_m;
}

// 寧偺帇嵎傪寁嶼偡傞
// T : 宱夁儐儕僂僗擭
Float64 AstronomicalMath::Moon::CalcParallax(Float64 T) const
{
	Float64 p;

	p = 0.9507 * sin(DegreeToRadian(90.0))
			+ 0.0518 * sin(DegreeToRadian(224.98 + 4771.989 * T))
			+ 0.0095 * sin(DegreeToRadian(190.7 + 4133.35 * T))
			+ 0.0078 * sin(DegreeToRadian(325.7 + 8905.34 * T))
			+ 0.0028 * sin(DegreeToRadian(0.0 + 9543.98 * T))
			+ 0.0009 * sin(DegreeToRadian(100.0 + 13677.3 * T))
			+ 0.0005 * sin(DegreeToRadian(329.0 + 8545.4 * T))
			+ 0.0004 * sin(DegreeToRadian(194.0 + 3773.4 * T))
			+ 0.0003 * sin(DegreeToRadian(227.0 + 4412.0 * T));

	return p;
}

// 惎偺弌杤帪崗傪寁嶼偡傞偨傔偵昁梫側忦審傪愝掕偡傞
// gmtDate : 帪崗偺弶婜抣丅偙偙偐傜 24 帪娫偺斖埻偱寁嶼偡傞
// position : 娤應抧揰偺宱搙丒堒搙
// H : 娤應抧揰偺崅搙 (m)
// R : 戝婥嵎
Void AstronomicalMath::StarRiseSet::Set(DateConstRef gmtDate, AstronomicalPointConstRef position, Float64 H, Float64 R)
{
	_gmtDate = gmtDate;
	_E = 0.0353333 * sqrt(H);
	_R = R;
	_position = position;
	_riseDate.Clear();
	_setDate.Clear();
}

// 惎偺弌杤帪崗傪寁嶼偡傞偨傔偵昁梫側忦審傪愝掕偡傞
// gmtDate : 帪崗偺弶婜抣丅偙偙偐傜 24 帪娫偺斖埻偱寁嶼偡傞
// lambda : 娤應抧揰偺宱搙
// phi : 娤應抧揰偺堒搙
// position : 娤應抧揰偺宱搙丒堒搙
// H : 娤應抧揰偺崅搙 (m)
// R : 戝婥嵎
Void AstronomicalMath::StarRiseSet::Set(DateConstRef gmtDate, Float64 lambda, Float64 phi, Float64 H, Float64 R)
{
	AstronomicalPoint position;

	position.Longitude = lambda;
	position.Latitude = phi;
	Set(gmtDate, position, H, R);
}

AstronomicalMath::StarRiseSet::StarRiseSet(Void) : _E(0.0), _R(0.0)
{
	_position.Longitude = 0.0;
	_position.Latitude = 0.0;
}

// gmtDate : 帪崗偺弶婜抣丅偙偙偐傜 24 帪娫偺斖埻偱寁嶼偡傞
// position : 娤應抧揰偺宱搙丒堒搙
// H : 娤應抧揰偺崅搙 (m)
// R : 戝婥嵎
AstronomicalMath::StarRiseSet::StarRiseSet(DateConstRef gmtDate, AstronomicalPointConstRef position, Float64 H, Float64 R)
{
	Set(gmtDate, position, H, R);
}

// gmtDate : 帪崗偺弶婜抣丅偙偙偐傜 24 帪娫偺斖埻偱寁嶼偡傞
// lambda : 娤應抧揰偺宱搙
// phi : 娤應抧揰偺堒搙
// position : 娤應抧揰偺宱搙丒堒搙
// H : 娤應抧揰偺崅搙 (m)
// R : 戝婥嵎
AstronomicalMath::StarRiseSet::StarRiseSet(DateConstRef gmtDate, Float64 lambda, Float64 phi, Float64 H, Float64 R)
{
	Set(gmtDate, lambda, phi, H, R);
}

// 巜掕偝傟偨擔偺弌杤帪崗傪寁嶼偡傞
Void AstronomicalMath::StarRiseSet::Compute(Void)
{
	Date x[5], e[7];
	Float64 y[7];
	SInt32 i, e_count(0);

	_riseDate.Clear();
	_setDate.Clear();

	x[0] = _gmtDate;
	for (i = 1; i < 5; ++i) {
		x[i] = x[i - 1] + DateDuration(0, 6, 0, 0);
	}
	for (i = 0; i < 5; ++i) {
		y[i] = DiffEquation(x[i]);
	}

	e[e_count++] = x[0];
	for (i = 0; i < 4; ++i) {
		if (y[i + 1] == 0.0) {
			e[e_count++] = x[i + 1];
			continue;
		}
		if (DifferentSign(y[i], y[i + 1])) {
#if defined TARGET_COMPILER_ARMCPP
			e[e_count++] = SolveEquation(&StarRiseSet::DiffEquation, x[i], x[i + 1]);
#else
			e[e_count++] = SolveEquation(DiffEquation, x[i], x[i + 1]);
#endif
		}
	}
	if (e[e_count - 1] != x[4]) {
		e[e_count++] = x[4];
	}

	for (i = 0; i < e_count; ++i) {
		y[i] = Equation(e[i]);
	}
	--e_count;
	for (i = 0; i < e_count; ++i) {
		if (y[i] < 0.0 && y[i + 1] > 0.0) {
#if defined TARGET_COMPILER_ARMCPP
			_riseDate.Append(SolveEquation(&StarRiseSet::Equation, e[i], e[i + 1]));
#else
			_riseDate.Append(SolveEquation(Equation, e[i], e[i + 1]));
#endif
		}
		else if (y[i] > 0.0 && y[i + 1] < 0.0) {
#if defined TARGET_COMPILER_ARMCPP
			_setDate.Append(SolveEquation(&StarRiseSet::Equation, e[i], e[i + 1]));
#else
			_setDate.Append(SolveEquation(Equation, e[i], e[i + 1]));
#endif
		}
	}
}

// 梌偊傜傟偨曽掱幃傪 2 暘朄偱夝偔
AstronomicalMath::Date AstronomicalMath::StarRiseSet::SolveEquation(PFUNC_EQUATION equation, DateConstRef initial, DateConstRef final) const
{
	Date x1(initial), x2, x3(final);
	DateDuration d1, d2;
	Float64 y1, y2, y3;

	y1 = (this->*equation)(initial);
	if (y1 == 0.0) {
		return initial;
	}
	y3 = (this->*equation)(final);
	if (y3 == 0.0) {
		return final;
	}
	if (!DifferentSign(y1, y3)) {
		return initial;
	}

	for ( ; ; ) {
		x2 = Date::GetMidpoint(x1, x3);
		d1 = x2 - x1;
		d2 = x3 - x2;
		if (d1.Get() < ASTRONOMICALMATH_EQ_EPS || d2.Get() < ASTRONOMICALMATH_EQ_EPS) {
			break;
		}
		y2 = (this->*equation)(x2);
		if (y2 == 0.0) {
			break;
		}
		if (DifferentSign(y1, y2)) {
			x3 = x2;
			y3 = y2;
		}
		else {
			x1 = x2;
			y1 = y2;
		}
	}
	
	return x2;
}

AstronomicalMath::SunRiseSunSet::SunRiseSunSet(Void)
{
}

AstronomicalMath::SunRiseSunSet::SunRiseSunSet(DateConstRef gmtDate, AstronomicalPointConstRef position, Float64 H, Float64 R) : StarRiseSet(gmtDate, position, H, R)
{
}

AstronomicalMath::SunRiseSunSet::SunRiseSunSet(DateConstRef gmtDate, Float64 lambda, Float64 phi, Float64 H, Float64 R) : StarRiseSet(gmtDate, lambda, phi, H, R)
{
}

// 擔弌杤曽掱幃
Float64 AstronomicalMath::SunRiseSunSet::Equation(DateConstRef gmtDate) const
{
	Sun sun(gmtDate);
	AstronomicalPoint position;
	Float64 Theta, t, delta, phi, k;

	position = sun.GetEquatorialPosition();		// 懢梲偺愒摴嵗昗
	Theta = GetSiderealTime(sun._gmtDate, _position.Longitude /* 娤應抧揰偺宱搙 */);
	t = DegreeToRadian(Theta - position.Longitude);
	delta = DegreeToRadian(position.Latitude);
	phi = DegreeToRadian(_position.Latitude /* 娤應抧揰偺堒搙*/);
	k = DegreeToRadian(- sun._S - _E - _R + sun._PI);

	return sin(delta) * sin(phi) + cos(delta) * cos(phi) * cos(t) - sin(k);
}

AstronomicalMath::MoonRiseMoonSet::MoonRiseMoonSet(Void)
{
}

AstronomicalMath::MoonRiseMoonSet::MoonRiseMoonSet(DateConstRef gmtDate, AstronomicalPointConstRef position, Float64 H, Float64 R) : StarRiseSet(gmtDate, position, H, R)
{
}

AstronomicalMath::MoonRiseMoonSet::MoonRiseMoonSet(DateConstRef gmtDate, Float64 lambda, Float64 phi, Float64 H, Float64 R) : StarRiseSet(gmtDate, lambda, phi, H, R)
{
}

// 寧弌杤曽掱幃
Float64 AstronomicalMath::MoonRiseMoonSet::Equation(DateConstRef gmtDate) const
{
	Moon moon(gmtDate);
	AstronomicalPoint position;
	Float64 Theta, t, delta, phi, k;

	position = moon.GetEquatorialPosition();	// 寧偺愒摴嵗昗
	Theta = GetSiderealTime(moon._gmtDate, _position.Longitude /* 娤應抧揰偺宱搙 */);
	t = DegreeToRadian(Theta - position.Longitude);
	delta = DegreeToRadian(position.Latitude);
	phi = DegreeToRadian(_position.Latitude /* 娤應抧揰偺堒搙 */);
	k = DegreeToRadian(- _E - _R + moon._PI);

	return sin(delta) * sin(phi) + cos(delta) * cos(phi) * cos(t) - sin(k);
}

// 峆惎帪傪庢摼偡傞
// lambda : 宱搙
Float64 AstronomicalMath::GetSiderealTime(DateConstRef gmtDate, Float64 lambda)
{
	Float64Const T(gmtDate.GetJ2000Year());
	Float64Const d(gmtDate.Get() - floor(gmtDate.Get()));
	Float64 theta;

	theta = 100.4606 + 360.007700536 * T + 0.00000003879 * T * T + 360.0 * d + lambda;
	theta = fmod(theta, 360.0);
	return (theta >= 0.0) ? theta : (theta + 360.0);
}

// 寧楊傪庢摼偡傞
Float64 AstronomicalMath::GetMoonAge(DateConstRef gmtDate)
{
	Sun sun(gmtDate);
	Moon moon(gmtDate);
	Date t0_date;
	Float64 d_lambda, t0, t;

	t0 = gmtDate.Get();
	t = t0;
	d_lambda = moon._eclipticPosition.Longitude - sun._eclipticPosition.Longitude;
	if (d_lambda < 0.0) {
		d_lambda += 360.0;
	}
	while (d_lambda >= ASTRONOMICALMATH_SUN_MOON_ANGLE_EPS) {
		t0 -= d_lambda / 12.1818;
		t0_date.Set(t0);
		sun.Set(t0_date);
		moon.Set(t0_date);
		d_lambda = moon._eclipticPosition.Longitude - sun._eclipticPosition.Longitude;
	}

	return t - t0;
}