1package Astro::MoonPhase; 2 3use strict; 4use vars qw($VERSION @ISA @EXPORT @EXPORT_OK); 5 6require Exporter; 7 8@ISA = qw(Exporter); 9@EXPORT = qw(phase phasehunt phaselist); 10$VERSION = '0.60'; 11 12use vars qw ( 13 $Epoch 14 $Elonge $Elongp $Eccent $Sunsmax $Sunangsiz 15 $Mmlong $Mmlongp $Mlnode $Minc $Mecc $Mangsiz $Msmax $Mparallax $Synmonth 16 $Pi 17 ); 18 19# Astronomical constants. 20 21$Epoch = 2444238.5; # 1980 January 0.0 22 23# Constants defining the Sun's apparent orbit. 24 25$Elonge = 278.833540; # ecliptic longitude of the Sun at epoch 1980.0 26$Elongp = 282.596403; # ecliptic longitude of the Sun at perigee 27$Eccent = 0.016718; # eccentricity of Earth's orbit 28$Sunsmax = 1.495985e8; # semi-major axis of Earth's orbit, km 29$Sunangsiz = 0.533128; # sun's angular size, degrees, at semi-major axis distance 30 31# Elements of the Moon's orbit, epoch 1980.0. 32 33$Mmlong = 64.975464; # moon's mean longitude at the epoch 34$Mmlongp = 349.383063; # mean longitude of the perigee at the epoch 35$Mlnode = 151.950429; # mean longitude of the node at the epoch 36$Minc = 5.145396; # inclination of the Moon's orbit 37$Mecc = 0.054900; # eccentricity of the Moon's orbit 38$Mangsiz = 0.5181; # moon's angular size at distance a from Earth 39$Msmax = 384401.0; # semi-major axis of Moon's orbit in km 40$Mparallax = 0.9507; # parallax at distance a from Earth 41$Synmonth = 29.53058868; # synodic month (new Moon to new Moon) 42 43# Properties of the Earth. 44 45$Pi = 3.14159265358979323846; # assume not near black hole nor in Tennessee 46 47# Handy mathematical functions. 48 49sub sgn { return (($_[0] < 0) ? -1 : ($_[0] > 0 ? 1 : 0)); } # extract sign 50sub fixangle { return ($_[0] - 360.0 * (floor($_[0] / 360.0))); } # fix angle 51sub torad { return ($_[0] * ($Pi / 180.0)); } # deg->rad 52sub todeg { return ($_[0] * (180.0 / $Pi)); } # rad->deg 53sub dsin { return (sin(torad($_[0]))); } # sin from deg 54sub dcos { return (cos(torad($_[0]))); } # cos from deg 55 56sub tan { return sin($_[0])/cos($_[0]); } 57sub asin { return ($_[0]<-1 or $_[0]>1) ? undef : atan2($_[0],sqrt(1-$_[0]*$_[0])); } 58sub atan { 59 if ($_[0]==0) { return 0; } 60 elsif ($_[0]>0) { return atan2(sqrt(1+$_[0]*$_[0]),sqrt(1+1/($_[0]*$_[0]))); } 61 else { return -atan2(sqrt(1+$_[0]*$_[0]),sqrt(1+1/($_[0]*$_[0]))); } 62} 63 64sub floor { 65 my $val = shift; 66 my $neg = $val < 0; 67 my $asint = int($val); 68 my $exact = $val == $asint; 69 70 return ($exact ? $asint : $neg ? $asint - 1 : $asint); 71} 72 73 74 75# jtime - convert internal date and time to astronomical Julian 76# time (i.e. Julian date plus day fraction) 77 78 79sub jtime { 80 81 my $t = shift; 82 my ($julian); 83 $julian = ($t / 86400) + 2440587.5; # (seconds /(seconds per day)) + julian date of epoch 84 return ($julian); 85 86} 87 88 89# jdaytosecs - convert Julian date to a UNIX epoch 90 91sub jdaytosecs { 92 my $jday = shift; 93 my $stamp; 94 95 $stamp = ($jday - 2440587.5)*86400; # (juliandate - jdate of unix epoch)*(seconds per julian day) 96 return($stamp); 97 98} 99 100 101 102# jyear - convert Julian date to year, month, day, which are 103# returned via integer pointers to integers 104sub jyear { 105 106 my ($td, $yy, $mm, $dd) = @_; 107 my ($z, $f, $a, $alpha, $b, $c, $d, $e); 108 109 $td += 0.5; # astronomical to civil 110 $z = floor($td); 111 $f = $td - $z; 112 113 if ($z < 2299161.0) { 114 $a = $z; 115 } else { 116 $alpha = floor(($z - 1867216.25) / 36524.25); 117 $a = $z + 1 + $alpha - floor($alpha / 4); 118 } 119 120 121 $b = $a + 1524; 122 $c = floor(($b - 122.1) / 365.25); 123 $d = floor(365.25 * $c); 124 $e = floor(($b - $d) / 30.6001); 125 126 $$dd = $b - $d - floor(30.6001 * $e) + $f; 127 $$mm = $e < 14 ? $e - 1 : $e - 13; 128 $$yy = $$mm > 2 ? $c - 4716 : $c - 4715; 129 130} 131 132## meanphase -- Calculates time of the mean new Moon for a given 133## base date. This argument K to this function is the 134## precomputed synodic month index, given by: 135## 136## K = (year - 1900) * 12.3685 137## 138## where year is expressed as a year and fractional year. 139 140 141sub meanphase { 142 my ($sdate, $k) = @_; 143 my ($t, $t2, $t3, $nt1); 144 145 ## Time in Julian centuries from 1900 January 0.5 146 $t = ($sdate - 2415020.0) / 36525; 147 $t2 = $t * $t; ## Square for frequent use 148 $t3 = $t2 * $t; ## Cube for frequent use 149 150 $nt1 = 2415020.75933 + $Synmonth * $k 151 + 0.0001178 * $t2 152 - 0.000000155 * $t3 153 + 0.00033 * dsin(166.56 + 132.87 * $t - 0.009173 * $t2); 154 155 return ($nt1); 156} 157 158 159# truephase - given a K value used to determine the mean phase of the 160# new moon, and a phase selector (0.0, 0.25, 0.5, 0.75), 161# obtain the true, corrected phase time 162 163sub truephase { 164 my ($k, $phase) = @_; 165 my ($t, $t2, $t3, $pt, $m, $mprime, $f); 166 my $apcor = 0; 167 168 $k += $phase; # add phase to new moon time 169 $t = $k / 1236.85; # time in Julian centuries from 170 # 1900 January 0.5 171 $t2 = $t * $t; # square for frequent use 172 $t3 = $t2 * $t; # cube for frequent use 173 174 # mean time of phase */ 175 $pt = 2415020.75933 176 + $Synmonth * $k 177 + 0.0001178 * $t2 178 - 0.000000155 * $t3 179 + 0.00033 * dsin(166.56 + 132.87 * $t - 0.009173 * $t2); 180 181 # Sun's mean anomaly 182 $m = 359.2242 183 + 29.10535608 * $k 184 - 0.0000333 * $t2 185 - 0.00000347 * $t3; 186 187 # Moon's mean anomaly 188 $mprime = 306.0253 189 + 385.81691806 * $k 190 + 0.0107306 * $t2 191 + 0.00001236 * $t3; 192 193 # Moon's argument of latitude 194 $f = 21.2964 195 + 390.67050646 * $k 196 - 0.0016528 * $t2 197 - 0.00000239 * $t3; 198 199 if (($phase < 0.01) || (abs($phase - 0.5) < 0.01)) { 200 # Corrections for New and Full Moon. 201 202 $pt += (0.1734 - 0.000393 * $t) * dsin($m) 203 + 0.0021 * dsin(2 * $m) 204 - 0.4068 * dsin($mprime) 205 + 0.0161 * dsin(2 * $mprime) 206 - 0.0004 * dsin(3 * $mprime) 207 + 0.0104 * dsin(2 * $f) 208 - 0.0051 * dsin($m + $mprime) 209 - 0.0074 * dsin($m - $mprime) 210 + 0.0004 * dsin(2 * $f + $m) 211 - 0.0004 * dsin(2 * $f - $m) 212 - 0.0006 * dsin(2 * $f + $mprime) 213 + 0.0010 * dsin(2 * $f - $mprime) 214 + 0.0005 * dsin($m + 2 * $mprime); 215 $apcor = 1; 216 } 217 elsif ((abs($phase - 0.25) < 0.01 || (abs($phase - 0.75) < 0.01))) { 218 $pt += (0.1721 - 0.0004 * $t) * dsin($m) 219 + 0.0021 * dsin(2 * $m) 220 - 0.6280 * dsin($mprime) 221 + 0.0089 * dsin(2 * $mprime) 222 - 0.0004 * dsin(3 * $mprime) 223 + 0.0079 * dsin(2 * $f) 224 - 0.0119 * dsin($m + $mprime) 225 - 0.0047 * dsin($m - $mprime) 226 + 0.0003 * dsin(2 * $f + $m) 227 - 0.0004 * dsin(2 * $f - $m) 228 - 0.0006 * dsin(2 * $f + $mprime) 229 + 0.0021 * dsin(2 * $f - $mprime) 230 + 0.0003 * dsin($m + 2 * $mprime) 231 + 0.0004 * dsin($m - 2 * $mprime) 232 - 0.0003 * dsin(2 * $m + $mprime); 233 if ($phase < 0.5) { 234 # First quarter correction. 235 $pt += 0.0028 - 0.0004 * dcos($m) + 0.0003 * dcos($mprime); 236 } 237 else { 238 # Last quarter correction. 239 $pt += -0.0028 + 0.0004 * dcos($m) - 0.0003 * dcos($mprime); 240 } 241 $apcor = 1; 242 } 243 if (!$apcor) { 244 die "truephase() called with invalid phase selector ($phase).\n"; 245 } 246 return ($pt); 247} 248 249# phasehunt - find time of phases of the moon which surround the current 250# date. Five phases are found, starting and ending with the 251# new moons which bound the current lunation 252 253sub phasehunt { 254 my $sdate = jtime(shift || time()); 255 my ($adate, $k1, $k2, $nt1, $nt2); 256 my ($yy, $mm, $dd); 257 258 $adate = $sdate - 45; 259 260 jyear($adate, \$yy, \$mm, \$dd); 261 $k1 = floor(($yy + (($mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685); 262 263 $adate = $nt1 = meanphase($adate, $k1); 264 265 while (1) { 266 $adate += $Synmonth; 267 $k2 = $k1 + 1; 268 $nt2 = meanphase($adate, $k2); 269 if (($nt1 <= $sdate) && ($nt2 > $sdate)) { 270 last; 271 } 272 $nt1 = $nt2; 273 $k1 = $k2; 274 275 } 276 277 278 279 return ( 280 jdaytosecs(truephase($k1, 0.0)), 281 jdaytosecs(truephase($k1, 0.25)), 282 jdaytosecs(truephase($k1, 0.5)), 283 jdaytosecs(truephase($k1, 0.75)), 284 jdaytosecs(truephase($k2, 0.0)) 285 ); 286} 287 288 289 290# phaselist - find time of phases of the moon between two dates 291# times (in & out) are seconds_since_1970 292 293sub phaselist 294{ 295 my ($sdate, $edate) = map { jtime($_) } @_; 296 297 my (@phases, $d, $k, $yy, $mm); 298 299 jyear($sdate, \$yy, \$mm, \$d); 300 $k = floor(($yy + (($mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685) - 2; 301 302 while (1) { 303 ++$k; 304 for my $phase (0.0, 0.25, 0.5, 0.75) { 305 $d = truephase($k, $phase); 306 307 return @phases if $d >= $edate; 308 309 if ($d >= $sdate) { 310 push @phases, int(4 * $phase) unless @phases; 311 push @phases, jdaytosecs($d); 312 } # end if date should be listed 313 } # end for each $phase 314 } # end while 1 315} # end phaselist 316 317 318 319# kepler - solve the equation of Kepler 320 321sub kepler { 322 my ($m, $ecc) = @_; 323 my ($e, $delta); 324 my $EPSILON = 1e-6; 325 326 $m = torad($m); 327 $e = $m; 328 do { 329 $delta = $e - $ecc * sin($e) - $m; 330 $e -= $delta / (1 - $ecc * cos($e)); 331 } while (abs($delta) > $EPSILON); 332 return ($e); 333} 334 335 336 337# phase - calculate phase of moon as a fraction: 338# 339# The argument is the time for which the phase is requested, 340# expressed as a Julian date and fraction. Returns the terminator 341# phase angle as a percentage of a full circle (i.e., 0 to 1), 342# and stores into pointer arguments the illuminated fraction of 343# the Moon's disc, the Moon's age in days and fraction, the 344# distance of the Moon from the centre of the Earth, and the 345# angular diameter subtended by the Moon as seen by an observer 346# at the centre of the Earth. 347 348sub phase { 349 my $pdate = jtime(shift || time()); 350 351 my $pphase; # illuminated fraction 352 my $mage; # age of moon in days 353 my $dist; # distance in kilometres 354 my $angdia; # angular diameter in degrees 355 my $sudist; # distance to Sun 356 my $suangdia; # sun's angular diameter 357 358 my ($Day, $N, $M, $Ec, $Lambdasun, $ml, $MM, $MN, $Ev, $Ae, $A3, $MmP, 359 $mEc, $A4, $lP, $V, $lPP, $NP, $y, $x, $Lambdamoon, $BetaM, 360 $MoonAge, $MoonPhase, 361 $MoonDist, $MoonDFrac, $MoonAng, $MoonPar, 362 $F, $SunDist, $SunAng, 363 $mpfrac); 364 365 # Calculation of the Sun's position. 366 367 $Day = $pdate - $Epoch; # date within epoch 368 $N = fixangle((360 / 365.2422) * $Day); # mean anomaly of the Sun 369 $M = fixangle($N + $Elonge - $Elongp); # convert from perigee 370 # co-ordinates to epoch 1980.0 371 $Ec = kepler($M, $Eccent); # solve equation of Kepler 372 $Ec = sqrt((1 + $Eccent) / (1 - $Eccent)) * tan($Ec / 2); 373 $Ec = 2 * todeg(atan($Ec)); # true anomaly 374 $Lambdasun = fixangle($Ec + $Elongp); # Sun's geocentric ecliptic 375 # longitude 376 # Orbital distance factor. 377 $F = ((1 + $Eccent * cos(torad($Ec))) / (1 - $Eccent * $Eccent)); 378 $SunDist = $Sunsmax / $F; # distance to Sun in km 379 $SunAng = $F * $Sunangsiz; # Sun's angular size in degrees 380 381 382 # Calculation of the Moon's position. 383 384 # Moon's mean longitude. 385 $ml = fixangle(13.1763966 * $Day + $Mmlong); 386 387 # Moon's mean anomaly. 388 $MM = fixangle($ml - 0.1114041 * $Day - $Mmlongp); 389 390 # Moon's ascending node mean longitude. 391 $MN = fixangle($Mlnode - 0.0529539 * $Day); 392 393 # Evection. 394 $Ev = 1.2739 * sin(torad(2 * ($ml - $Lambdasun) - $MM)); 395 396 # Annual equation. 397 $Ae = 0.1858 * sin(torad($M)); 398 399 # Correction term. 400 $A3 = 0.37 * sin(torad($M)); 401 402 # Corrected anomaly. 403 $MmP = $MM + $Ev - $Ae - $A3; 404 405 # Correction for the equation of the centre. 406 $mEc = 6.2886 * sin(torad($MmP)); 407 408 # Another correction term. 409 $A4 = 0.214 * sin(torad(2 * $MmP)); 410 411 # Corrected longitude. 412 $lP = $ml + $Ev + $mEc - $Ae + $A4; 413 414 # Variation. 415 $V = 0.6583 * sin(torad(2 * ($lP - $Lambdasun))); 416 417 # True longitude. 418 $lPP = $lP + $V; 419 420 # Corrected longitude of the node. 421 $NP = $MN - 0.16 * sin(torad($M)); 422 423 # Y inclination coordinate. 424 $y = sin(torad($lPP - $NP)) * cos(torad($Minc)); 425 426 # X inclination coordinate. 427 $x = cos(torad($lPP - $NP)); 428 429 # Ecliptic longitude. 430 $Lambdamoon = todeg(atan2($y, $x)); 431 $Lambdamoon += $NP; 432 433 # Ecliptic latitude. 434 $BetaM = todeg(asin(sin(torad($lPP - $NP)) * sin(torad($Minc)))); 435 436 # Calculation of the phase of the Moon. 437 438 # Age of the Moon in degrees. 439 $MoonAge = $lPP - $Lambdasun; 440 441 # Phase of the Moon. 442 $MoonPhase = (1 - cos(torad($MoonAge))) / 2; 443 444 # Calculate distance of moon from the centre of the Earth. 445 446 $MoonDist = ($Msmax * (1 - $Mecc * $Mecc)) / 447 (1 + $Mecc * cos(torad($MmP + $mEc))); 448 449 # Calculate Moon's angular diameter. 450 451 $MoonDFrac = $MoonDist / $Msmax; 452 $MoonAng = $Mangsiz / $MoonDFrac; 453 454 # Calculate Moon's parallax. 455 456 $MoonPar = $Mparallax / $MoonDFrac; 457 458 $pphase = $MoonPhase; 459 $mage = $Synmonth * (fixangle($MoonAge) / 360.0); 460 $dist = $MoonDist; 461 $angdia = $MoonAng; 462 $sudist = $SunDist; 463 $suangdia = $SunAng; 464 $mpfrac = fixangle($MoonAge) / 360.0; 465 return wantarray ? ( $mpfrac, $pphase, $mage, $dist, $angdia, $sudist,$suangdia ) : $mpfrac; 466} 467 4681; 469__END__ 470 471=head1 NAME 472 473Astro::MoonPhase - Information about the phase of the Moon 474 475=head1 SYNOPSIS 476 477use Astro::MoonPhase; 478 479 ( $MoonPhase, 480 $MoonIllum, 481 $MoonAge, 482 $MoonDist, 483 $MoonAng, 484 $SunDist, 485 $SunAng ) = phase($seconds_since_1970); 486 487 @phases = phasehunt($seconds_since_1970); 488 489 ($phase, @times) = phaselist($start, $stop); 490 491=head1 DESCRIPTION 492 493MoonPhase calculates information about the phase of the moon 494at a given time. 495 496=head1 FUNCTIONS 497 498=head2 phase() 499 500 ( $MoonPhase, 501 $MoonIllum, 502 $MoonAge, 503 $MoonDist, 504 $MoonAng, 505 $SunDist, 506 $SunAng ) = phase($seconds_since_1970); 507 508 $MoonPhase = phase($seconds_since_1970); 509 510The argument is the time for which the phase is requested, 511expressed as a time returned by the C<time> function. If C<$seconds_since_1970> 512is omitted, it does C<phase(time)>. 513 514Return value in scalar context is $MoonPhase, 515the terminator phase angle as a percentage of a full circle (i.e., 0 to 1). 516 517=over 4 518 519=item B<Return values in array context:> 520 521=item $MoonPhase: 522 523the terminator phase angle as a percentage of a full circle (i.e., 0 to 1) 524 525=item $MoonIllum: 526 527the illuminated fraction of the Moon's disc 528 529=item $MoonAge: 530 531the Moon's age in days and fraction 532 533=item $MoonDist: 534 535the distance of the Moon from the centre of the Earth 536 537=item $MoonAng: 538 539the angular diameter subtended by the Moon as seen by 540an observer at the centre of the Earth. 541 542=item $SunDist: 543 544the distance from the Sun in km 545 546=item $SunAng: 547 548the angular size of Sun in degrees 549 550=back 551 552Example: 553 554 ( $MoonPhase, 555 $MoonIllum, 556 $MoonAge, 557 $MoonDist, 558 $MoonAng, 559 $SunDist, 560 $SunAng ) = phase(); 561 562 print "MoonPhase = $MoonPhase\n"; 563 print "MoonIllum = $MoonIllum\n"; 564 print "MoonAge = $MoonAge\n"; 565 print "MoonDist = $MoonDist\n"; 566 print "MoonAng = $MoonAng\n"; 567 print "SunDist = $SunDist\n"; 568 print "SunAng = $SunAng\n"; 569 570could print something like this: 571 572 MoonPhase = 0.598939375319023 573 MoonIllum = 0.906458030827876 574 MoonAge = 17.6870323368022 575 MoonDist = 372479.357420033 576 MoonAng = 0.534682403555093 577 SunDist = 152078368.820205 578 SunAng = 0.524434538105092 579 580=head2 phasehunt() 581 582 @phases = phasehunt($seconds_since_1970); 583 584Finds time of phases of the moon which surround the given 585date. Five phases are found, starting and ending with the 586new moons which bound the current lunation. 587 588The argument is the time, expressed as a time returned 589by the C<time> function. If C<$seconds_since_1970> 590is omitted, it does C<phasehunt(time)>. 591 592Example: 593 594 @phases = phasehunt(); 595 print "New Moon = ", scalar(localtime($phases[0])), "\n"; 596 print "First quarter = ", scalar(localtime($phases[1])), "\n"; 597 print "Full moon = ", scalar(localtime($phases[2])), "\n"; 598 print "Last quarter = ", scalar(localtime($phases[3])), "\n"; 599 print "New Moon = ", scalar(localtime($phases[4])), "\n"; 600 601could print something like this: 602 603 New Moon = Wed Jun 24 06:51:47 1998 604 First quarter = Wed Jul 1 21:42:19 1998 605 Full moon = Thu Jul 9 19:02:47 1998 606 Last quarter = Thu Jul 16 18:15:18 1998 607 New Moon = Thu Jul 23 16:45:01 1998 608 609=head2 phaselist() 610 611 ($phase, @times) = phaselist($start, $stop); 612 613Finds times of all phases of the moon which occur on or after 614C<$start> but before C<$stop>. Both the arguments and the return 615values are expressed as seconds since 1970 (like the C<time> function 616returns). 617 618C<$phase> is an integer indicating the phase of the moon at 619C<$times[0]>, as shown in this table: 620 621 0 New Moon 622 1 First quarter 623 2 Full Moon 624 3 Last quarter 625 626The remaining values in C<@times> indicate subsequent phases of the 627moon (in ascending order by time). If there are no phases of the moon 628between C<$start> and C<$stop>, C<phaselist> returns the empty list. 629 630Example: 631 632 @name = ("New Moon", "First quarter", "Full moon", "Last quarter"); 633 ($phase, @times) = phaselist($start, $stop); 634 635 while (@times) { 636 printf "%-14s= %s\n", $name[$phase], scalar localtime shift @times; 637 $phase = ($phase + 1) % 4; 638 } 639 640could produce the same output as the C<phasehunt> example above (given 641the appropriate start & stop times). 642 643=head1 ABOUT THE ALGORITHMS 644 645The algorithms used in this program to calculate the positions of Sun and 646Moon as seen from the Earth are given in the book I<Practical Astronomy 647With Your Calculator> by B<Peter Duffett-Smith, Second Edition, 648Cambridge University Press, 1981>. Ignore the word "Calculator" in the 649title; this is an essential reference if you're interested in 650developing software which calculates planetary positions, orbits, 651eclipses, and the like. If you're interested in pursuing such 652programming, you should also obtain: 653 654I<Astronomical Formulae for Calculators> by B<Jean Meeus, Third Edition, 655Willmann-Bell, 1985>. A must-have. 656 657I<Planetary Programs and Tables from -4000 to +2800> by B<Pierre 658Bretagnon and Jean-Louis Simon, Willmann-Bell, 1986>. If you want the 659utmost (outside of JPL) accuracy for the planets, it's here. 660 661I<Celestial BASIC> by B<Eric Burgess, Revised Edition, Sybex, 1985>. Very 662cookbook oriented, and many of the algorithms are hard to dig out of 663the turgid BASIC code, but you'll probably want it anyway. 664 665Many of these references can be obtained from Willmann-Bell, P.O. Box 66635025, Richmond, VA 23235, USA. Phone: (804) 320-7016. In addition 667to their own publications, they stock most of the standard references 668for mathematical and positional astronomy. 669 670=head1 LICENCE 671 672This program is in the public domain: "Do what thou wilt shall be the 673whole of the law". 674 675=head1 AUTHORS 676 677The moontool.c Release 2.0: 678 679 A Moon for the Sun 680 Designed and implemented by John Walker in December 1987, 681 revised and updated in February of 1988. 682 683Initial Perl transcription: 684 685 Raino Pikkarainen, 1998 686 raino.pikkarainen@saunalahti.fi 687 688The moontool.c Release 2.4: 689 690 Major enhancements by Ron Hitchens, 1989 691 692Revisions: 693 694 Brett Hamilton http://simple.be/ 695 Bug fix, 2003 696 Second transcription and bugfixes, 2004 697 698 Christopher J. Madsen http://www.cjmweb.net/ 699 Added phaselist function, March 2007 700