// 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' );
    // 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.
    //FR
    $names = array( 'Nouvelle lune', 'Premier croissant', 'Premier quartier', 'Gibbeuse croissante', 'Pleine lune', 'Gibbeuse décroissante', 'Dernier quartier', 'Dernier croissant', 'Nouvelle lune' );
    return $names[ floor( ( $this->phase + 0.0625 ) * 8 ) ];
  }
}

// Récupérations des données voulues pour affichage :
// create an instance of the class, and use the current time
setlocale(LC_TIME, 'fr','fr_FR','fr_FR@euro','fr_FR.utf8','fr-FR','fra');
$moon = new MoonPhase();
$quelle_phase = $moon->phase_name();
$age = round( $moon->age(), 1 );
$stage = $moon->phase() < 0.5 ? 'Montante' : 'Descendante';
$distance = round( $moon->distance(), 0 );
$illumination = round($moon->illumination()*100, 0);

echo round($moon->phase()*100,0) . '% du cycle';
echo '<br>';
echo $quelle_phase;			// Phase (quartier, pleine, croissant, etc.)
echo "<br>";				// Age de la lune dans le cycle :
if ($age < 2) { echo $age." jour"; } else { echo $age." jours"; }
echo "<br>";
echo $stage;				// Etat : montante ou descendante
echo "<br>";
echo "Visibilité : ".$illumination." %";
echo "<br>";
echo "Distance : ".$distance." kms";
echo '<br>';

// Recherche des prochaines pleines et nouvelles lunes :
$next_full =  strftime( '%a %d/%m %R', $moon->full_moon() );
$next_full2 = strftime( '%a %d/%m %R', $moon->next_full_moon() );
$next_new =   strftime( '%a %d/%m %R', $moon->new_moon() );
$next_new2 =  strftime( '%a %d/%m %R', $moon->next_new_moon() );

//date de la prochaine pleine lune - format linux
$pleine_lune = strtotime(strftime( '%d.%m.%Y', $moon->full_moon() ));
//nouveau code pour gérer correctement la date de prochaine pleine lune renvoyée.
if ( $pleine_lune <= strtotime('now') ) //si la date de la prochaine pleine lune calculée est dépassée
{
  $date_pleine_lune = $next_full2; //alors on utilise la date de pleine lune suivante
}
else {
  $date_pleine_lune = $next_full; //sinon on prend la date de pleine lune la plus proche
}
echo "Pleine : ". $date_pleine_lune;
echo '<br>';

//date de la prochaine nouvelle lune - format linux
$nouvelle_lune = strtotime(strftime( '%d.%m.%Y', $moon->new_moon() ));
//nouveau code pour gérer correctement la date de prochaine nouvelle lune renvoyée.
if ( $nouvelle_lune <= strtotime('now') ) //si la date de la prochaine nouvelle lune calculée est dépassée
{
  $date_nouvelle_lune = $next_new2; //alors on utilise la date de nouvelle lune suivante
}
else {
  $date_nouvelle_lune = $next_new; //sinon on prend la date de nouvelle lune la plus proche
}
echo "Nouvelle : ". $date_nouvelle_lune;

// On affecte le contenu de $xxxxxxx à la variable « xxxxxxxx » que l’on pourra afficher par exemple dans un virtuel
$scenario->setData('v_lune_age',$age);								// Age de la lune (le jour du cycle)
$scenario->setData('v_lune_pourcent',round($moon->phase()*100,0));	// Pourcentage du cycle (~29.53 jours)
$scenario->setData('v_lune_phase',$quelle_phase);					// Phase (quartier, pleine, croissant, etc.)
$scenario->setData('v_lune_etat',$stage);							// Etat : montante ou descendante
$scenario->setData('v_lune_visibilite',$illumination);				// Pourcentage de visibilité
$scenario->setData('v_lune_distance',$distance);					// Distance depuis la Terre
$scenario->setData('v_lune_next_pleine',$date_pleine_lune);			// Prochaine pleine lune
$scenario->setData('v_lune_next_nouv',$date_nouvelle_lune);			// Prochaine nouvelle lune