Class Calcurate Sunrise / Sunset time

Suntime.php

<?php
/**
 * Class Calcurate Sunrise / Sunset time
 *
 * PHP version 5
 *
 * Licence:
 *
 * Copyright (c) 2006-2009, Ippei Suzuki <[email protected]>
 * All rights reserved.
 *
 * Sample:
 *
 *   $s = new Suntime(135.7717, 35.008, time());
 *   $sunRiseTime = date('Y-m-d H:i:s', $s->getSunRise());
 *   $sunSetTime  = date('Y-m-d H:i:s', $s->getSunSet());
 *
 * This Sample, you can get Sunrise time / Sunset time at Kyoto, Japan.
 *
 * Option:
 *
 * When you use option with getSunRise/getSunSet, you can get "Twilight Time".
 * Twilight Time is defined as the Sun's height in -7h21m40s
 *
 * Option Sample:
 *
 *   $s = new Suntime(135.7717, 35.008, time());
 *   $twilightSunRiseTime = date('Y-m-d H:i:s', $s->getSunRise(1));
 *   $twilightSunSetTime  = date('Y-m-d H:i:s', $s->getSunSet(1));
 *
 * Special Thanks:
 *
 * Math Logic of this program is based on "日の出・日の入りの計算" by Ko Nagasawa, Japan.
 * (see also: http://ja.wikipedia.org/wiki/%E9%95%B7%E6%B2%A2%E5%B7%A5)
 *
 * @version 1.0.0
 * @author Ippei Suzuki
 * @link http://www.nobdoyplace.com
 */
class Suntime
{
  /**
   * Circle ratio(radian)
   */
  const PAI = 0.017453278;
  /**
   * 太陽視半径
   */
  const SO  = 0.266994;
  /**
   * 大気差
   */
  const R   = 0.585556;
  /**
   * 視差
   */
  const PO  = 0.0024428;
  /**
   * 地球自転遅れの補正(day)
   */
  const DT  = 0.00074074074;

  /**
   * target time
   *
   * @var unixtime
   */
  protected $targettime;
  /**
   * longitude
   *
   * @var float
   */
  protected $longitude;
  /**
   * latitude
   *
   * @var float
   */
  protected $latitude;
  /**
   * Request needs Twilight Time or Not
   *
   * @var integer
   */
  protected $isTwilight;

  /**
   * factory
   * default: Today at Kyoto, Japan
   *
   * @param unixtime $unixtime
   * @param float $x
   * @param float $y
   * @return object
   */
  public function __construct($longitude=null, $latitude=null, $targettime=null)
  {
    $this->setLongitude($longitude);
    $this->setLatitude($latitude);
    $this->setTargetTime($targettime);
  }
  /**
   * Sets the Time to target
   *
   * @param unixtime $targettime
   * @return Suntime
   */
  public function setTargetTime($targettime)
  {
    $this->targettime = (is_null($targettime) === true || is_numeric($targettime) === false) ?
      time() :
      (int) $targettime;
    return $this;
  }
  /**
   * Returns the Time to target
   *
   * @return integer unixtime
   */
  public function getTargetTime()
  {
    return $this->targettime;
  }
  /**
   * Sets the Longitude('keido' in Japanese)
   *
   * @param float $longitude
   * @return Suntime
   */
  public function setLongitude($longitude)
  {
    $this->longitude = (is_null($longitude) === true || is_float($longitude) === false) ?
      135.7717 :
      (float) $longitude;
    return $this;
  }
  /**
   * Returns the Longitude('keido' in Japanese)
   *
   * @return float
   */
  public function getLongitude()
  {
    return $this->longitude;
  }
  /**
   * Sets the Latitude('ido' in Japanese)
   *
   * @param float $latitude
   * @return Suntime
   */
  public function setLatitude($latitude)
  {
    $this->latitude = (is_null($latitude) === true || is_float($latitude) === false) ?
      35.008 :
      $latitude;
    return $this;
  }
  /**
   * Returnd the Latitude('ido' in Japanese)
   *
   * @return float
   */
  public function getLatitude()
  {
    return $this->latitude;
  }
  /**
   * Sets request needs Twilight('hakumei' in Japanese) or not
   *
   * @param integer $isTwilight 1|0
   * @return Suntime
   */
  public function setIsTwilight($isTwilight)
  {
    $this->isTwilight = ($isTwilight === 1) ? 1 : 0;
    return $this;
  }
  /**
   * Returns request needs Twilight('hakumei' in Japanese) or not
   *
   * @return integer
   */
  public function getIsTwilight()
  {
    return $this->isTwilight;
  }



  /**
   * Returns Time of Sun Rise
   *
   * @param integer $isTwilight
   * @return integer|unixtime
   */
  public function getSunRise($isTwilight=0)
  {
    $this->setIsTwilight($isTwilight);
    return $this->calc(true);
  }
  /**
   * Returns Time of Sun Set
   *
   * @param integer $isTwilight
   * @return integer|unixtime
   */
  public function getSunSet($isTwilight=0)
  {
    $this->setIsTwilight($isTwilight);
    return $this->calc(false);
  }





  /**
   * Returns the Sun Time
   *
   * @param boolean $rise Calculates Sun-Rise or Not
   * @return integer|unixtime
   */
  protected function calc($rise=true)
  {
  //check for Time to target
    if (is_null($this->targettime) === true):
      throw new Exception('Target Time is null');
    endif;
  //check for Time to target
    if (is_null($this->longitude) === true):
      throw new Exception('Longitude is null');
    endif;
  //check for Time to target
    if (is_null($this->latitude) === true):
      throw new Exception('Latitude is null');
    endif;

    $Ko = $this->getKo();
  //temporary result(day)
    $d = (float) (($rise === true) ? 6 / 24 : 17 / 24);
  //use the method of approximation
    $dd = 1;
    while($dd > 0.00005):
      $T = (float) (($Ko + $d + self::DT) / 365.25);
    //太陽黄経（単位：度）
      $ramda = $this->getRamda($T);
    //太陽距離
      $r = $this->getR($T);
    //黄道傾角
      $e = $this->getE($T);
    //太陽赤経
      $alpha = $this->getAlpha($ramda, $e);
    //太陽赤緯
      $delta = $this->getDelta($ramda, $e);
    //恒星時
      $theta = $this->getTheta($T, $d);
    //太陽出没高度[k]
      $k = $this->getK($r);
    //日の出入高度の時角
      $tk = $this->getTk($k, $delta, $rise);
    //仮定時刻の時角
      $t_tmp = $theta - $alpha;
    //仮定時刻dに対する補正値
      $dd = ($tk - $t_tmp) / 360;
      $dd -= floor($d);
    //補正
      $d += $dd;
    endwhile;
  //追加補正
    if($d >= 1)
      $d -= floor($d);
  //整形
    $H = $d * 24;
    $M = ($H - floor($H)) * 60;
    $hour = floor($H);
    $min = floor($M);
    $second = floor(($M - floor($M)) * 60);
    if($second >= 60){$second -= 60;$min++;}
    if($min >= 60){$min -= 60;$hour++;}

    return mktime($hour, $min, $second, date('n', $this->targettime), date('j', $this->targettime), date('Y', $this->targettime));
  }





  /**
   * Returns the number of days since J2000.0(2000/01/01 in Japan)
   *
   * @return float
   */
  protected function getKo()
  {
    $year  = ((int) date('n', $this->targettime) < 3) ? (int) date('Y', $this->targettime) - 2001 : (int) date('Y', $this->targettime) - 2000;
    $month = ((int) date('n', $this->targettime) < 3) ? (int) date('n', $this->targettime) +12 : (int) date('n', $this->targettime);
    $day   = (int) date('j', $this->targettime);
    return (float) ($year * 365 + $month * 30 + $day - 33.875 + floor(3 * ($month + 1) / 5) + floor($year / 4));
  }
  /**
   * 太陽黄経
   * 0 <= $ramda < 360
   *
   * @param float $T
   * @return float
   */
  protected function getRamda($T)
  {
    $ramda = 280.4603 + 360.00769 * $T
      + (1.9146 - 0.00005 * $T) * sin((357.538 + 359.991 * $T) * self::PAI)
      + 0.0200 * sin((355.05 + 719.981 * $T) * self::PAI)
      + 0.0048 * sin((234.95 + 19.341 * $T) * self::PAI)
      + 0.0020 * sin((247.1 + 329.64 * $T) * self::PAI)
      + 0.0018 * sin((297.8 + 4452.67 * $T) * self::PAI)
      + 0.0018 * sin((251.3 + 0.20 * $T) * self::PAI)
      + 0.0015 * sin((343.2 + 450.37 * $T) * self::PAI)
      + 0.0013 * sin((81.4 + 225.18 * $T) * self::PAI)
      + 0.0008 * sin((132.5 + 659.29 * $T) * self::PAI)
      + 0.0007 * sin((153.3 + 90.38 * $T) * self::PAI)
      + 0.0007 * sin((206.8 + 30.35 * $T) * self::PAI)
      + 0.0006 * sin((29.8 + 337.18 * $T) * self::PAI)
      + 0.0005 * sin((207.4 + 1.50 * $T) * self::PAI)
      + 0.0005 * sin((291.2 + 22.81 * $T) * self::PAI)
      + 0.0004 * sin((234.9 + 315.56 * $T) * self::PAI)
      + 0.0004 * sin((157.3 + 299.30 * $T) * self::PAI)
      + 0.0004 * sin((21.1 + 720.02 * $T) * self::PAI)
      + 0.0003 * sin((352.5 + 1079.97 * $T) * self::PAI)
      + 0.0003 * sin((329.7 + 44.43 * $T) * self::PAI);
    if($ramda > 0):
      while($ramda >= 360):
        $ramda -= 360;
      endwhile;
    else:
      while($ramda < 0):
        $ramda += 360;
      endwhile;
    endif;
    return $ramda;
  }
  /**
   * 太陽距離
   *
   * @param float $T
   * @return integer
   */
  protected function getR($T)
  {
    $q = (0.007256 - 0.0000002 * $T) * sin((267.54 + 359.991 * $T) * self::PAI)
      + 0.000091 * sin((265.1 + 719.98 * $T) * self::PAI)
      + 0.000030 * sin(90.0 * self::PAI)
      + 0.000013 * sin((27.8 + 4452.67 * $T) * self::PAI)
      + 0.000007 * sin((254 + 450.4 * $T) * self::PAI)
      + 0.000007 * sin((156 + 329.6 * $T) * self::PAI);
    return pow(10, $q);
  }
  /**
   * 黄道傾角
   *
   * @param float $T
   * @return float
   */
  protected function getE($T)
  {
    return (float) (23.439291 - 0.000130042 * $T);
  }

// 太陽赤経[alpha]の計算
  /**
   *  太陽赤経
   *
   * @param float $ramda
   * @param float $e
   * @return float
   */
  protected function getAlpha($ramda, $e)
  {
    $alpha = atan(tan($ramda * self::PAI) * cos($e * self::PAI)) / self::PAI;
    if ($ramda < 180):
      if ($alpha < 0):
        while ($alpha < 0):
          $alpha += 180;
        endwhile;
      else:
        while ($alpha >= 180):
          $alpha -= 180;
        endwhile;
      endif;
    else:
      if ($alpha < 180):
        while ($alpha < 180):
          $alpha += 180;
        endwhile;
      else:
        while ($alpha >= 360):
          $alpha -= 180;
        endwhile;
      endif;
    endif;
    return $alpha;
  }
  /**
   * 太陽赤緯
   *
   * @param float $ramda
   * @param float $e
   * @return float
   */
  protected function getDelta($ramda, $e)
  {
    return asin(sin($ramda * self::PAI) * sin($e * self::PAI)) / self::PAI;
  }
  /**
   * 恒星時
   * 0 <= $this->theta < 360
   *
   * @param float $T
   * @param float $time
   * @return float
   */
  protected function getTheta($T, $time)
  {
    $theta = 325.4606 + 360.007700536 * $T + 0.00000003879 * $T * $T + 360 * $time + $this->longitude;
    if ($theta < 0):
      while ($theta < 0):
        $theta += 360;
      endwhile;
    elseif ($theta >= 360):
      while ($theta >= 360):
        $theta -= 360;
      endwhile;
    endif;
    return $theta;
  }
  /**
   * 太陽出没高度
   *
   * @param integer $r
   * @return float
   */
  protected function getK($r)
  {
    $S = self::SO / $r;
    $P = self::PO / $r;
    return (float) (($this->isTwilight) ? $P - 7.36111 : $P - $S - self::R);
  }
  /**
   * 日の出高度の時角
   *
   * @param float $k
   * @return float
   */
  protected function getTk($k, $delta, $rise)
  {
    $costk = (sin($k * self::PAI) - sin($delta * self::PAI) * sin($this->latitude * self::PAI)) / (cos($delta * self::PAI) * cos($this->latitude * self::PAI));
    $tk = acos($costk) / self::PAI;
    return ($rise === true) ? (($tk > 0) ? $tk * -1 : $tk) : (($tk < 0) ? $tk * -1 : $tk);
  }
}
?>