MajorDoMo

<?php
/**
* Moon phase calculation class
* Adapted for PHP from Moontool for Windows (http://www.fourmilab.ch/moontoolw/)
* by Samir Shah (http://rayofsolaris.net)
* License: MIT
**/
//namespace Solaris;
class MoonPhase {
private $timestamp;
private $phase;
private $illum;
private $age;
private $dist;
private $angdia;
private $sundist;
private $sunangdia;
private $synmonth;
private $quarters = null;
function __construct( $pdate = null ) {
if( is_null( $pdate ) )
$pdate = time();
/* Astronomical constants */
$epoch = 2444238.5; // 1980 January 0.0
/* Constants defining the Sun's apparent orbit */
$elonge = 278.833540; // Ecliptic longitude of the Sun at epoch 1980.0
$elongp = 282.596403; // Ecliptic longitude of the Sun at perigee
$eccent = 0.016718; // Eccentricity of Earth's orbit
$sunsmax = 1.495985e8; // Semi-major axis of Earth's orbit, km
$sunangsiz = 0.533128; // Sun's angular size, degrees, at semi-major axis distance
/* Elements of the Moon's orbit, epoch 1980.0 */
$mmlong = 64.975464; // Moon's mean longitude at the epoch
$mmlongp = 349.383063; // Mean longitude of the perigee at the epoch
$mlnode = 151.950429; // Mean longitude of the node at the epoch
$minc = 5.145396; // Inclination of the Moon's orbit
$mecc = 0.054900; // Eccentricity of the Moon's orbit
$mangsiz = 0.5181; // Moon's angular size at distance a from Earth
$msmax = 384401; // Semi-major axis of Moon's orbit in km
$mparallax = 0.9507; // Parallax at distance a from Earth
$synmonth = 29.53058868; // Synodic month (new Moon to new Moon)
$this->synmonth = $synmonth;
$lunatbase = 2423436.0; // Base date for E. W. Brown's numbered series of lunations (1923 January 16)
/* Properties of the Earth */
// $earthrad = 6378.16; // Radius of Earth in kilometres
// $PI = 3.14159265358979323846; // Assume not near black hole
$this->timestamp = $pdate;
// pdate is coming in as a UNIX timstamp, so convert it to Julian
$pdate = $pdate / 86400 + 2440587.5;
/* Calculation of the Sun's position */
$Day = $pdate - $epoch; // Date within epoch
$N = $this->fixangle((360 / 365.2422) * $Day); // Mean anomaly of the Sun
$M = $this->fixangle($N + $elonge - $elongp); // Convert from perigee co-ordinates to epoch 1980.0
$Ec = $this->kepler($M, $eccent); // Solve equation of Kepler
$Ec = sqrt((1 + $eccent) / (1 - $eccent)) * tan($Ec / 2);
$Ec = 2 * rad2deg(atan($Ec)); // True anomaly
$Lambdasun = $this->fixangle($Ec + $elongp); // Sun's geocentric ecliptic longitude
$F = ((1 + $eccent * cos(deg2rad($Ec))) / (1 - $eccent * $eccent)); // Orbital distance factor
$SunDist = $sunsmax / $F; // Distance to Sun in km
$SunAng = $F * $sunangsiz; // Sun's angular size in degrees
/* Calculation of the Moon's position */
$ml = $this->fixangle(13.1763966 * $Day + $mmlong); // Moon's mean longitude
$MM = $this->fixangle($ml - 0.1114041 * $Day - $mmlongp); // Moon's mean anomaly
$MN = $this->fixangle($mlnode - 0.0529539 * $Day); // Moon's ascending node mean longitude
$Ev = 1.2739 * sin(deg2rad(2 * ($ml - $Lambdasun) - $MM)); // Evection
$Ae = 0.1858 * sin(deg2rad($M)); // Annual equation
$A3 = 0.37 * sin(deg2rad($M)); // Correction term
$MmP = $MM + $Ev - $Ae - $A3; // Corrected anomaly
$mEc = 6.2886 * sin(deg2rad($MmP)); // Correction for the equation of the centre
$A4 = 0.214 * sin(deg2rad(2 * $MmP)); // Another correction term
$lP = $ml + $Ev + $mEc - $Ae + $A4; // Corrected longitude
$V = 0.6583 * sin(deg2rad(2 * ($lP - $Lambdasun))); // Variation
$lPP = $lP + $V; // True longitude
$NP = $MN - 0.16 * sin(deg2rad($M)); // Corrected longitude of the node
$y = sin(deg2rad($lPP - $NP)) * cos(deg2rad($minc)); // Y inclination coordinate
$x = cos(deg2rad($lPP - $NP)); // X inclination coordinate
$Lambdamoon = rad2deg(atan2($y, $x)) + $NP; // Ecliptic longitude
$BetaM = rad2deg(asin(sin(deg2rad($lPP - $NP)) * sin(deg2rad($minc)))); // Ecliptic latitude
/* Calculation of the phase of the Moon */
$MoonAge = $lPP - $Lambdasun; // Age of the Moon in degrees
$MoonPhase = (1 - cos(deg2rad($MoonAge))) / 2;
// Phase of the Moon
// Distance of moon from the centre of the Earth
$MoonDist = ($msmax * (1 - $mecc * $mecc)) / (1 + $mecc * cos(deg2rad($MmP + $mEc)));
$MoonDFrac = $MoonDist / $msmax;
$MoonAng = $mangsiz / $MoonDFrac;
// Moon's angular diameter
// $MoonPar = $mparallax / $MoonDFrac; // Moon's parallax
// store results
$this->phase = $this->fixangle($MoonAge) / 360; // Phase (0 to 1)
$this->illum = $MoonPhase; // Illuminated fraction (0 to 1)
$this->age = $synmonth * $this->phase; // Age of moon (days)
$this->dist = $MoonDist; // Distance (kilometres)
$this->angdia = $MoonAng; // Angular diameter (degrees)
$this->sundist = $SunDist; // Distance to Sun (kilometres)
$this->sunangdia = $SunAng; // Sun's angular diameter (degrees)
}
private function fixangle($a) {
return ( $a - 360 * floor($a / 360) );
}
// KEPLER -- Solve the equation of Kepler.
private function kepler($m, $ecc) {
$epsilon = 0.000001; // 1E-6
$e = $m = deg2rad($m);
do {
$delta = $e - $ecc * sin($e) - $m;
$e -= $delta / ( 1 - $ecc * cos($e) );
}
while ( abs($delta) > $epsilon );
return $e;
}
/* Calculates time of the mean new Moon for a given
base date. This argument K to this function is the
precomputed synodic month index, given by:
K = (year - 1900) * 12.3685
where year is expressed as a year and fractional year.
*/
private function meanphase($sdate, $k){
// Time in Julian centuries from 1900 January 0.5
$t = ( $sdate - 2415020.0 ) / 36525;
$t2 = $t * $t;
$t3 = $t2 * $t;
$nt1 = 2415020.75933 + $this->synmonth * $k
+ 0.0001178 * $t2
- 0.000000155 * $t3
+ 0.00033 * sin( deg2rad( 166.56 + 132.87 * $t - 0.009173 * $t2 ) );
return $nt1;
}
/* Given a K value used to determine the mean phase of
the new moon, and a phase selector (0.0, 0.25, 0.5,
0.75), obtain the true, corrected phase time.
*/
private function truephase($k, $phase){
$apcor = false;
$k += $phase; // Add phase to new moon time
$t = $k / 1236.85; // Time in Julian centuries from 1900 January 0.5
$t2 = $t * $t; // Square for frequent use
$t3 = $t2 * $t; // Cube for frequent use
$pt = 2415020.75933 // Mean time of phase
+ $this->synmonth * $k
+ 0.0001178 * $t2
- 0.000000155 * $t3
+ 0.00033 * sin( deg2rad( 166.56 + 132.87 * $t - 0.009173 * $t2 ) );
$m = 359.2242 + 29.10535608 * $k - 0.0000333 * $t2 - 0.00000347 * $t3; // Sun's mean anomaly
$mprime = 306.0253 + 385.81691806 * $k + 0.0107306 * $t2 + 0.00001236 * $t3; // Moon's mean anomaly
$f = 21.2964 + 390.67050646 * $k - 0.0016528 * $t2 - 0.00000239 * $t3; // Moon's argument of latitude
if ( $phase < 0.01 || abs( $phase - 0.5 ) < 0.01 ) {
// Corrections for New and Full Moon
$pt += (0.1734 - 0.000393 * $t) * sin( deg2rad( $m ) )
+ 0.0021 * sin( deg2rad( 2 * $m ) )
- 0.4068 * sin( deg2rad( $mprime ) )
+ 0.0161 * sin( deg2rad( 2 * $mprime) )
- 0.0004 * sin( deg2rad( 3 * $mprime ) )
+ 0.0104 * sin( deg2rad( 2 * $f ) )
- 0.0051 * sin( deg2rad( $m + $mprime ) )
- 0.0074 * sin( deg2rad( $m - $mprime ) )
+ 0.0004 * sin( deg2rad( 2 * $f + $m ) )
- 0.0004 * sin( deg2rad( 2 * $f - $m ) )
- 0.0006 * sin( deg2rad( 2 * $f + $mprime ) )
+ 0.0010 * sin( deg2rad( 2 * $f - $mprime ) )
+ 0.0005 * sin( deg2rad( $m + 2 * $mprime ) );
$apcor = true;
} else if ( abs( $phase - 0.25 ) < 0.01 || abs( $phase - 0.75 ) < 0.01 ) {
$pt += (0.1721 - 0.0004 * $t) * sin( deg2rad( $m ) )
+ 0.0021 * sin( deg2rad( 2 * $m ) )
- 0.6280 * sin( deg2rad( $mprime ) )
+ 0.0089 * sin( deg2rad( 2 * $mprime) )
- 0.0004 * sin( deg2rad( 3 * $mprime ) )
+ 0.0079 * sin( deg2rad( 2 * $f ) )
- 0.0119 * sin( deg2rad( $m + $mprime ) )
- 0.0047 * sin( deg2rad ( $m - $mprime ) )
+ 0.0003 * sin( deg2rad( 2 * $f + $m ) )
- 0.0004 * sin( deg2rad( 2 * $f - $m ) )
- 0.0006 * sin( deg2rad( 2 * $f + $mprime ) )
+ 0.0021 * sin( deg2rad( 2 * $f - $mprime ) )
+ 0.0003 * sin( deg2rad( $m + 2 * $mprime ) )
+ 0.0004 * sin( deg2rad( $m - 2 * $mprime ) )
- 0.0003 * sin( deg2rad( 2 * $m + $mprime ) );
if ( $phase < 0.5 ) // First quarter correction
$pt += 0.0028 - 0.0004 * cos( deg2rad( $m ) ) + 0.0003 * cos( deg2rad( $mprime ) );
else // Last quarter correction
$pt += -0.0028 + 0.0004 * cos( deg2rad( $m ) ) - 0.0003 * cos( deg2rad( $mprime ) );
$apcor = true;
}
if (!$apcor) // function was called with an invalid phase selector
return false;
return $pt;
}
/* Find time of phases of the moon which surround the current date.
Five phases are found, starting and
ending with the new moons which bound the current lunation.
*/
private function phasehunt() {
$sdate = $this->utctojulian( $this->timestamp );
$adate = $sdate - 45;
$ats = $this->timestamp - 86400 * 45;
$yy = (int) gmdate( 'Y', $ats );
$mm = (int) gmdate( 'n', $ats );
$k1 = floor( ( $yy + ( ( $mm - 1 ) * ( 1 / 12 ) ) - 1900 ) * 12.3685 );
$adate = $nt1 = $this->meanphase( $adate, $k1 );
while (true) {
$adate += $this->synmonth;
$k2 = $k1 + 1;
$nt2 = $this->meanphase( $adate, $k2 );
// if nt2 is close to sdate, then mean phase isn't good enough, we have to be more accurate
if( abs( $nt2 - $sdate ) < 0.75 )
$nt2 = $this->truephase( $k2, 0.0 );
if ( $nt1 <= $sdate && $nt2 > $sdate )
break;
$nt1 = $nt2;
$k1 = $k2;
}
// results in Julian dates
$data = array(
$this->truephase( $k1, 0.0 ),
$this->truephase( $k1, 0.25 ),
$this->truephase( $k1, 0.5 ),
$this->truephase( $k1, 0.75 ),
$this->truephase( $k2, 0.0 ),
$this->truephase( $k2, 0.25 ),
$this->truephase( $k2, 0.5 ),
$this->truephase( $k2, 0.75 )
);
$this->quarters = array();
foreach( $data as $v )
$this->quarters[] = ( $v - 2440587.5 ) * 86400; // convert to UNIX time
}
/* Convert UNIX timestamp to astronomical Julian time (i.e. Julian date plus day fraction). */
private function utctojulian( $ts ) {
return $ts / 86400 + 2440587.5;
}
private function get_phase( $n ) {
if( is_null( $this->quarters ) )
$this->phasehunt();
return $this->quarters[$n];
}
/* Public functions for accessing results */
function phase(){
return $this->phase;
}
function illumination(){
return $this->illum;
}
function age(){
return $this->age;
}
function distance(){
return $this->dist;
}
function diameter(){
return $this->angdia;
}
function sundistance(){
return $this->sundist;
}
function sundiameter(){
return $this->sunangdia;
}
function new_moon(){
return $this->get_phase( 0 );
}
function first_quarter(){
return $this->get_phase( 1 );
}
function full_moon(){
return $this->get_phase( 2 );
}
function last_quarter(){
return $this->get_phase( 3 );
}
function next_new_moon(){
return $this->get_phase( 4 );
}
function next_first_quarter(){
return $this->get_phase( 5 );
}
function next_full_moon(){
return $this->get_phase( 6 );
}
function next_last_quarter(){
return $this->get_phase( 7 );
}
function phase_name() {
//$names = array( 'New Moon', 'Waxing Crescent', 'First Quarter', 'Waxing Gibbous', 'Full Moon', 'Waning Gibbous', 'Third Quarter', 'Waning Crescent', 'New Moon' );
$names = array( 'Новолуние', 'Растущая Луна', 'Первый квартал', 'Растущая Луна', 'Полнолуние', 'Убывающая луна', 'Третий квартал', 'Убывающий полумесяц', 'Новолуние');
// There are eight phases, evenly split. A "New Moon" occupies the 1/16th phases either side of phase = 0, and the rest follow from that.
return $names[ floor( ( $this->phase + 0.0625 ) * 8 ) ];
}
}

class Moon
{
function phase($Year, $Month, $Day, $Hour, $Minutes, $Seconds)
{
$DateSec = time($Hour, $Minutes, $Seconds, $Month, $Day, $Year, 0);

ini_set(precision, "20"); //Defini la precision des calcules

# Astronomical constants.
$Epoch = 2444238.5; # 1980 January 0.0

# Constants defining the Sun's apparent orbit.
$Elonge = 278.833540; # ecliptic longitude of the Sun at epoch 1980.0
$Elongp = 282.596403; # ecliptic longitude of the Sun at perigee
$Eccent = 0.016718; # eccentricity of Earth's orbit
$Sunsmax = 1.495985e8; # semi-major axis of Earth's orbit, km
$Sunangsiz = 0.533128; # sun's angular size, degrees, at semi-major axis distance

# Elements of the Moon's orbit, epoch 1980.0.
$Mmlong = 64.975464; # moon's mean longitude at the epoch
$Mmlongp = 349.383063; # mean longitude of the perigee at the epoch
$Mlnode = 151.950429; # mean longitude of the node at the epoch
$Minc = 5.145396; # inclination of the Moon's orbit
$Mecc = 0.054900; # eccentricity of the Moon's orbit
$Mangsiz = 0.5181; # moon's angular size at distance a from Earth
$Msmax = 384401.0; # semi-major axis of Moon's orbit in km
$Mparallax = 0.9507; # parallax at distance a from Earth
$Synmonth = 29.53058868; # synodic month (new Moon to new Moon)

$pdate = Moon::jtime($DateSec);

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


# Calculation of the Sun's position.

$Day = $pdate - $Epoch; # date within epoch
$N = Moon::fixangle((360 / 365.2422) * $Day); # mean anomaly of the Sun
$M = Moon::fixangle($N + $Elonge - $Elongp); # convert from perigee
# co-ordinates to epoch 1980.0
$Ec = Moon::kepler($M, $Eccent); # solve equation of Kepler
$Ec = sqrt((1 + $Eccent) / (1 - $Eccent)) * tan($Ec / 2);
$Ec = 2 * Moon::todeg(atan($Ec)); # true anomaly
$Lambdasun = Moon::fixangle($Ec + $Elongp); # Sun's geocentric ecliptic
# longitude
# Orbital distance factor.
$F = ((1 + $Eccent * cos(Moon::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 = Moon::fixangle(13.1763966 * $Day + $Mmlong);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# Ecliptic latitude.
$BetaM = Moon::todeg(asin(sin(Moon::torad($lPP - $NP)) * sin(Moon::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(Moon::torad($MoonAge))) / 2;

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

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

# Calculate Moon's angular diameter.

$MoonDFrac = $MoonDist / $Msmax;
$MoonAng = $Mangsiz / $MoonDFrac;

# Calculate Moon's parallax.

$MoonPar = $Mparallax / $MoonDFrac;

$pphase = $MoonPhase; # illuminated fraction
$mage = $Synmonth * (Moon::fixangle($MoonAge) / 360.0); # age of moon in days
$dist = $MoonDist; # distance in kilometres
$angdia = $MoonAng; # angular diameter in degrees
$sudist = $SunDist; # distance to Sun
$suangdia = $SunAng; # sun's angular diameter
$mpfrac = Moon::fixangle($MoonAge) / 360.0;
return array( $pphase, $mage, $dist, $angdia, $sudist, $suangdia, $mpfrac, $mpfrac );
}

function fixangle($x) { return ($x - 360.0 * (floor($x / 360.0))); } # fix angle
function torad($x) { return ($x * (M_PI / 180.0)); } # deg->rad
function todeg($x) { return ($x * (180.0 / M_PI)); } # rad->deg

function jtime($t)
{
$julian = ($t / 86400) + 2440587.5; # (seconds /(seconds per day)) + julian date of epoch 2440587.5 / 86400 = 28,24753472222 Days
return ($julian);
}

function kepler($m, $ecc)
{
$EPSILON = 1e-6;

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

}


//Exemple d'utilisation :

//Pour le 11 Avril 2009 à 00h00
list($MoonPhase, $MoonAge, $MoonDist, $MoonAng, $SunDist, $SunAng, $mpfrac,$names) = Moon::phase(2009, 04, 11, 00, 00, 01);
echo "Видимость луны ".number_format($MoonPhase*100, 2, ',', '')."%"."<br>";
echo "Ее возраст ".number_format($MoonAge, 0, ',', '')." дней"."<br>";
echo "Луна находится на расстоянии ".number_format($MoonDist, 0, ',', '')." km по отношению к Земле."."<br>";
echo "Ее возраст ".number_format($MoonAng, 0, ',', '')." дней"."<br>";
echo "Расстояние до Солнца ".number_format($SunDist, 0, ',', '')." km"."<br>";
echo "Ее возраст ".number_format($SunAng, 0, ',', '')." дней"."<br>";
echo "Ее возраст ".number_format($mpfrac, 0, ',', '')." дней"."<br>";
//sg("MoonPhase.MoonPhase",number_format($MoonPhase*100, 2, ',', ''));
//sg("MoonPhase.MoonAge",number_format($MoonAge, 0, ',', ''));
//sg("MoonPhase.MoonDist",number_format($MoonDist, 0, ',', ''));
//sg("MoonPhase.names",$names);

<SCRIPT LANGUAGE="JavaScript">
var i, b = new Date();b.setTime(b.getTime() + (b.getTimezoneOffset()*60000));var c = new Date(96, 1, 3, 16, 15, 0), d = 2551443000, e = (b.getTime() - c.getTime()) % d, f = Math.round((d-e)/(24*3600*1000)), sklon='дней';if ((f==1)||(f==21)||(f==31))sklon='день';if ((f>1&&f<5)||(f>21&&f<25))sklon='дня';document.write("<br> ",f," ",sklon,"");</script>