view release on metacpan or  search on metacpan

MoonPhase.pm  view on Meta::CPAN

	do {
		$delta = $e - $ecc * sin($e) - $m;
		$e -= $delta / (1 - $ecc * cos($e));
	} while (abs($delta) > $EPSILON);
	return ($e);
}



# phase - calculate phase of moon as a fraction:
#
# The argument is the time for which the phase is requested,
# expressed as a Julian date and fraction.  Returns the terminator
# phase angle as a percentage of a full circle (i.e., 0 to 1),
# and stores into pointer arguments the illuminated fraction of
# the Moon's disc, the Moon's age in days and fraction, the
# distance of the Moon from the centre of the Earth, and the
# angular diameter subtended by the Moon as seen by an observer
# at the centre of the Earth.

sub phase {
	my $pdate = jtime(shift || time());

	my $pphase;				# illuminated fraction
	my $mage;				# age of moon in days
	my $dist;				# distance in kilometres
	my $angdia;				# angular diameter in degrees
	my $sudist;				# distance to Sun
	my $suangdia;				# sun's angular diameter

	my ($Day, $N, $M, $Ec, $Lambdasun, $ml, $MM, $MN, $Ev, $Ae, $A3, $MmP,
	   $mEc, $A4, $lP, $V, $lPP, $NP, $y, $x, $Lambdamoon, $BetaM,
	   $MoonAge, $MoonPhase,
	   $MoonDist, $MoonDFrac, $MoonAng, $MoonPar,
	   $F, $SunDist, $SunAng,
	   $mpfrac);

	# Calculation of the Sun's position.

	$Day = $pdate - $Epoch;						# date within epoch
	$N = fixangle((360 / 365.2422) * $Day);				# mean anomaly of the Sun
	$M = fixangle($N + $Elonge - $Elongp);				# convert from perigee
									# co-ordinates to epoch 1980.0
	$Ec = kepler($M, $Eccent);					# solve equation of Kepler
	$Ec = sqrt((1 + $Eccent) / (1 - $Eccent)) * tan($Ec / 2);
	$Ec = 2 * todeg(atan($Ec));					# true anomaly
	$Lambdasun = fixangle($Ec + $Elongp);				# Sun's geocentric ecliptic
									# longitude
	# Orbital distance factor.
	$F = ((1 + $Eccent * cos(torad($Ec))) / (1 - $Eccent * $Eccent));
	$SunDist = $Sunsmax / $F;					# distance to Sun in km
	$SunAng = $F * $Sunangsiz;					# Sun's angular size in degrees


	# Calculation of the Moon's position.

	# Moon's mean longitude.
	$ml = fixangle(13.1763966 * $Day + $Mmlong);

	# Moon's mean anomaly.
	$MM = fixangle($ml - 0.1114041 * $Day - $Mmlongp);

	# Moon's ascending node mean longitude.
	$MN = fixangle($Mlnode - 0.0529539 * $Day);

	# Evection.
	$Ev = 1.2739 * sin(torad(2 * ($ml - $Lambdasun) - $MM));

	# Annual equation.
	$Ae = 0.1858 * sin(torad($M));

	# Correction term.
	$A3 = 0.37 * sin(torad($M));

	# Corrected anomaly.
	$MmP = $MM + $Ev - $Ae - $A3;

	# Correction for the equation of the centre.
	$mEc = 6.2886 * sin(torad($MmP));

	# Another correction term.
	$A4 = 0.214 * sin(torad(2 * $MmP));

	# Corrected longitude.
	$lP = $ml + $Ev + $mEc - $Ae + $A4;

	# Variation.
	$V = 0.6583 * sin(torad(2 * ($lP - $Lambdasun)));

	# True longitude.
	$lPP = $lP + $V;

	# Corrected longitude of the node.
	$NP = $MN - 0.16 * sin(torad($M));

	# Y inclination coordinate.
	$y = sin(torad($lPP - $NP)) * cos(torad($Minc));

	# X inclination coordinate.
	$x = cos(torad($lPP - $NP));

	# Ecliptic longitude.
	$Lambdamoon = todeg(atan2($y, $x));
	$Lambdamoon += $NP;

	# Ecliptic latitude.
	$BetaM = todeg(asin(sin(torad($lPP - $NP)) * sin(torad($Minc))));

	# Calculation of the phase of the Moon.

	# Age of the Moon in degrees.
	$MoonAge = $lPP - $Lambdasun;

	# Phase of the Moon.
	$MoonPhase = (1 - cos(torad($MoonAge))) / 2;

	# Calculate distance of moon from the centre of the Earth.

	$MoonDist = ($Msmax * (1 - $Mecc * $Mecc)) /
		(1 + $Mecc * cos(torad($MmP + $mEc)));