13
\$\begingroup\$

Photography is one of my hobbies, and, living in Canada, I have the possibility of seeing, and photographing, the Aurora Borealis (Northern Lights).

Now, these happen when there's ionization in the ionosphere typically resulting from strong solar winds which create a 'tear', resulting in the plasma visible as the 'Northern Lights'. The more intense the solar wind, the further south you can possibly see the 'tear'.

Since I live relatively far South in Canada, it will take a very strong solar event to cause the Northern lights to appear South of my position, so, it will inevitably be best seen looking North, when it happens.

Additionally, they are best seen when there are no clouds, and when the sky is as dark as possible.

Finally, the solar wind 'pushes' the Northern lights away from the sun, so they stretch further south on the 'dark' side of the earth, than the light side (they form an egg-shape, not a circle).

Due to the nature of the phenomenon, they can best be seen (for me):

  • close to midnight (when my part of the world is furthest away from the sun - and the pointy-side of the egg-shaped aurora is closest to me)
  • when there is no cloud cover
  • when there is lots of solar activity (kp (a measure of solar activity) >= 8)
  • I am in 'the country' to avoid light-pollution (instead of a city)
  • there is no moon creating light.

I have written a system that checks various factors, and e-mails me when the conditions are good for viewing the Aurora. Hopefully, when all the celestial factors align, I will be able to wake up, get out, and take pictures of the majestic events....

The solar phase calculation is the part that is relevant to this question. I need to know how dark it is likely to be, and as a result, I need to know what phase the moon is in (and in a different process, when the moon rises and sets). If the moon is dark (new moon), then it is a good thing. The closer I am to new moon the better.

The following code is an implementation of the function described in Wikipedia for calculating the Lunation phase. I only need an approximate formula. If I am off by an hour or two, it does not matter.

$$ d=5.597661 + 29.5305888610 \times N + (102.026\times10^{-12}) \times N^2 $$

where

$$ \begin{align} N & => \text{The n'th new moon since 1 Jan 2000}\\ 29.5305888610 & => \text{Lunation period - mean time between new moons (days)}\\ 102.026\times10^{-12} & => \text{Cyclical correction factor (from tides, etc)}\\ d & => \text{The number of days since 1 Jan 2000 to the Nth new moon}\\ \end{align} $$

The above formula is related to terrestrial time (not wall-clock time), and an adjustment needs to be made to convert to a wall-clock time (using UTC in this case)

$$ adjustment = 0.000739 - (235 \times 10^{-12}) \times N^2 $$

Now, this calculates the day offset of a new moon relative to 1 Jan 2000. I need to know relative to 'now'. The following Perl program prints out the number of days before, or after the new moon (depending on which is smaller).

#!/usr/bin/perl -w

use strict;

# see http://en.wikipedia.org/wiki/New_moon#Determining_new_moons:_an_approximate_formula

my $date = defined($ARGV[0]) ? $ARGV[0] : time();

my $y2kref = 5.597661;
my $lunation = 29.5305888610;
my $kfactor = 102.026 * 10**-12;
my $utos = 0.000739;
my $utk = 235 * 10**-12;
my $y2kepoch = 946702800;
my $secsinday = 60 * 60 * 24.0;

## we have had at least 84 lunation cycles since 1 jan 2000
my $moons = 84;
## days to/from new moon
my $days = 0;
## days since 1 Jan 2000
my $now = ($date - $y2kepoch) / $secsinday;

while ($days < $now) {

    $moons++;

    # calculate the new moon in terrestrial time
    $days = $y2kref + ($lunation * $moons) + ($kfactor * $moons * $moons);

    # Convert the terrestrial time to universal time.
    $days = $days - $utos - ($utk * $moons * $moons);

}

#$days is the number of days to the next new moon.


my $closest = $days - $now;
my $previous = $days - $lunation;
$closest = ($previous - $now) if ($now - $previous < $closest);

printf "%.3f\n", $closest;

If I run this using the following hack:

perl -e 'foreach (1..30) {$t = (time() + 86400 * $_); $d = localtime($t); print "echo -n $d \"   \" && bin/newmoon $t\n";};' | bash

I get the results:

Mon Dec 15 22:05:52 2014    5.835
Tue Dec 16 22:05:52 2014    4.835
Wed Dec 17 22:05:52 2014    3.835
Thu Dec 18 22:05:52 2014    2.835
Fri Dec 19 22:05:52 2014    1.835
Sat Dec 20 22:05:52 2014    0.835
Sun Dec 21 22:05:52 2014    -0.165
Mon Dec 22 22:05:52 2014    -1.165
Tue Dec 23 22:05:52 2014    -2.165
Wed Dec 24 22:05:52 2014    -3.165
Thu Dec 25 22:05:52 2014    -4.165
Fri Dec 26 22:05:52 2014    -5.165
Sat Dec 27 22:05:52 2014    -6.165
Sun Dec 28 22:05:52 2014    -7.165
Mon Dec 29 22:05:52 2014    -8.165
Tue Dec 30 22:05:52 2014    -9.165
Wed Dec 31 22:05:52 2014    -10.165
Thu Jan 1 22:05:52 2015    -11.165
Fri Jan 2 22:05:52 2015    -12.165
Sat Jan 3 22:05:52 2015    -13.165
Sun Jan 4 22:05:52 2015    -14.165
Mon Jan 5 22:05:52 2015    14.366
Tue Jan 6 22:05:52 2015    13.366
Wed Jan 7 22:05:52 2015    12.366
Thu Jan 8 22:05:52 2015    11.366
Fri Jan 9 22:05:52 2015    10.366
Sat Jan 10 22:05:52 2015    9.366
Sun Jan 11 22:05:52 2015    8.366
Mon Jan 12 22:05:52 2015    7.366
Tue Jan 13 22:05:52 2015    6.366

which shows a new moon between:

Sat Dec 20 22:05:52 2014    0.835
Sun Dec 21 22:05:52 2014    -0.165
\$\endgroup\$
  • 1
    \$\begingroup\$ The results from the formula differ by up to several hours relative to the US Naval Observatory calculations. \$\endgroup\$ – 200_success Dec 26 '14 at 10:44
  • \$\begingroup\$ You don't appear to have asked a question. But I suggest that Astro::MoonPhase will help you. \$\endgroup\$ – Borodin Feb 19 '15 at 12:00
  • \$\begingroup\$ @Borodin - yes, absolutely, it will. I did some searches for lunar cycle calculators, and so on. Perhaps I should have posed a question on Software Recommendations ;-) ? \$\endgroup\$ – rolfl Feb 19 '15 at 12:22
  • \$\begingroup\$ I don't think Software Recommendations is sufficiently granular to suggest language modules, but it would be wonderful if it did. The biggest shortcoming of CPAN is a lack of curation \$\endgroup\$ – Borodin Feb 20 '15 at 20:43
6
\$\begingroup\$

To accept user input for dates, DateTime::Format::Natural is a nice CPAN module to use. It spares you from having to input dates using the Unix epoch format and the tortured use of Perl-Bash-Perl to generate the table. (Defining a subroutine would also have helped you generate the whole table in Perl.)

Perl has some features for writing numeric literals that you should take advantage of. Instead of 102.026 * 10**-12, you could write 102.026E-12. You can also insert underscores, as in 29.530_588_861. Constants would be better written as

use constant {
    Y2K_EPOCH    => 946_702_800,
    Y2K_REF      => 5.597_661,

    LUNATION     => 29.530_588_861,
    K_FACTOR     => 102.026E-12,

    U_TO_S       => 0.000_739,
    UTK          => 235E-12,

    SECS_PER_DAY => 60 * 60 * 24,
};

Considering that you only want an approximate answer, I am perplexed by your concern for tiny numbers. \$(102.026 \times 10^{-12}) \times N\$ amounts to slowdown of an additional 9 millionths of a second per month. The \$(235 \times 10^{-12}) \times N^2\$ term works out to 20 millionths of a second per month per month. At this rate, it would take millions of years before the time of the new moon is delayed by a whole night due to drag forces!

Once you drop the insignificant \$N^2\$ terms, the quadratic expressions become linear. You no longer need to loop (nor do you need to solve a quadratic equation).

$now is a misnomer, considering that the time can be overridden on the command line.

To decide between the previous and the following month, just round off $moons to the nearest integer.

sub days_until_new_moon {
    my $date = shift // time;
    my $days = ($date - Y2K_EPOCH) / SECS_PER_DAY;

    # Use sprintf() for rounding: http://stackoverflow.com/q/178539
    my $moons = sprintf('%.0f', ($days - Y2K_REF) / LUNATION);
    my $newmoon_days = Y2K_REF + LUNATION * $moons;

    # Convert terrestrial time to universal time
    $newmoon_days -= U_TO_S;

    return $newmoon_days - $days;
}

say days_until_new_moon(@ARGV);
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.