view release on metacpan or  search on metacpan

MoonPhase.pm  view on Meta::CPAN

                        jdaytosecs(truephase($k1, 0.25)),
                        jdaytosecs(truephase($k1, 0.5)),
                        jdaytosecs(truephase($k1, 0.75)),
                        jdaytosecs(truephase($k2, 0.0))
                        );
}



# phaselist - find time of phases of the moon between two dates
# times (in & out) are seconds_since_1970

sub phaselist
{
  my ($sdate, $edate) = map { jtime($_) } @_;

  my (@phases, $d, $k, $yy, $mm);

  jyear($sdate, \$yy, \$mm, \$d);
  $k = floor(($yy + (($mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685) - 2;

  while (1) {
    ++$k;
    for my $phase (0.0, 0.25, 0.5, 0.75) {
      $d = truephase($k, $phase);

      return @phases if $d >= $edate;

      if ($d >= $sdate) {
        push @phases, int(4 * $phase) unless @phases;
        push @phases, jdaytosecs($d);
      } # end if date should be listed
    } # end for each $phase
  } # end while 1
} # end phaselist



# kepler - solve the equation of Kepler

sub kepler {
	my ($m, $ecc) = @_;
	my ($e, $delta);
	my $EPSILON = 1e-6;

	$m = torad($m);
	$e = $m;
	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));

MoonPhase.pm  view on Meta::CPAN

	$MoonPar = $Mparallax / $MoonDFrac;

	$pphase = $MoonPhase;
	$mage = $Synmonth * (fixangle($MoonAge) / 360.0);
	$dist = $MoonDist;
	$angdia = $MoonAng;
	$sudist = $SunDist;
	$suangdia = $SunAng;
	$mpfrac = fixangle($MoonAge) / 360.0;
	return wantarray ? ( $mpfrac, $pphase, $mage, $dist, $angdia, $sudist,$suangdia ) : $mpfrac;
}

1;
__END__

=head1 NAME

Astro::MoonPhase - Information about the phase of the Moon

=head1 SYNOPSIS

use Astro::MoonPhase;

	( $MoonPhase,
	  $MoonIllum,
	  $MoonAge,
	  $MoonDist,
	  $MoonAng,
	  $SunDist,
	  $SunAng ) = phase($seconds_since_1970);

	@phases  = phasehunt($seconds_since_1970);

	($phase, @times) = phaselist($start, $stop);

=head1 DESCRIPTION

MoonPhase calculates information about the phase of the moon
at a given time.

=head1 FUNCTIONS

=head2 phase()

	( $MoonPhase,
	  $MoonIllum,
	  $MoonAge,
	  $MoonDist,
	  $MoonAng,
	  $SunDist,
	  $SunAng )  = phase($seconds_since_1970);

	  $MoonPhase = phase($seconds_since_1970);

The argument is the time for which the phase is requested,
expressed as a time returned by the C<time> function. If C<$seconds_since_1970>
is omitted, it does C<phase(time)>.

Return value in scalar context is $MoonPhase,
the terminator phase angle as a percentage of a full circle (i.e., 0 to 1).

=over 4

=item B<Return values in array context:>

=item $MoonPhase:

the terminator phase angle as a percentage of a full circle (i.e., 0 to 1)

=item $MoonIllum:

the illuminated fraction of the Moon's disc

=item $MoonAge:

the Moon's age in days and fraction

=item $MoonDist:

the distance of the Moon from the centre of the Earth

=item $MoonAng:

the angular diameter subtended by the Moon as seen by
an observer at the centre of the Earth.

=item $SunDist:

the distance from the Sun in km

=item $SunAng:

the angular size of Sun in degrees

=back

Example:

   ( $MoonPhase,
     $MoonIllum,
     $MoonAge,
     $MoonDist,
     $MoonAng,
     $SunDist,
     $SunAng ) = phase();

     print "MoonPhase  = $MoonPhase\n";
     print "MoonIllum  = $MoonIllum\n";
     print "MoonAge    = $MoonAge\n";
     print "MoonDist   = $MoonDist\n";
     print "MoonAng    = $MoonAng\n";
     print "SunDist    = $SunDist\n";
     print "SunAng     = $SunAng\n";

could print something like this:

     MoonPhase  = 0.598939375319023
     MoonIllum  = 0.906458030827876
     MoonAge    = 17.6870323368022
     MoonDist   = 372479.357420033
     MoonAng    = 0.534682403555093
     SunDist    = 152078368.820205
     SunAng     = 0.524434538105092

=head2 phasehunt()

     @phases = phasehunt($seconds_since_1970);