Astronomy
<< SpaceInvadersApplet | FinalProjectsTrailIndex | Perspective >>
moonset.php5
<?php
class MiniMoon
{
public $time=0;
public $RA=0;
public $dec=0;
public function __construct($t){
$this->setTime($t);
}
public function setTime($t){
$P2=2*pi();
$cosEPS=0.91748;// cos(obliquity ecliptic)
$sinEPS=0.39778;// sin(obliquity ecliptic)
$arc = 206264.8062;
$this->time = $t;
// mean elements of lunar orbit
$l0 = frac(0.606433+1336.855225*$t); // mean longitude of moon in rev
$l = $P2*frac(0.374897+1325.552410*$t); // mean anomaly of the Moon
$ls = $P2*frac(0.993133+99.997361*$t); // mean anomaly of the Sun
$d = $P2*frac(0.827361+1236.853086*$t); // diff longitude Moon-Sun
$f = $P2*frac(0.259086+1342.227825*$t); // mean argument of latitude
$dl = 22640*sin($l) - 4586*sin($l-2*$d) + 2370* sin(2*$d) + 769*sin(2*$l) -668*sin($ls) - 412*sin(2*$f) - 212*sin(2*$l-2*$d) - 206*sin($l+$ls-2*$d) + 192*sin($l+2*$d) - 165*sin($ls-2*$d) - 125*sin($d) - 110*sin($l+$ls) + 148*sin($l-$ls) -55*sin(2*$f-2*$d);
$s = $f + ($dl+412*sin(2*$f)+541*sin($ls))/$arc;
$h = $f-2*$d;
$n = -516*sin($h) + 44*sin($l+$h) - 31*sin($h-$l) - 23*sin($ls+$h) +11*sin($h-$ls) - 25*sin($f-2*$l) + 21*sin($f-$l);
$bMoon = (18520.0*sin($s) + $n) /$arc; // in radians
$lMoon = $P2 * frac( $l0 +$dl/1296E3) ; // in radians
//echo "bMoon=".$bMoon." lMoon=".$lMoon."<br>";
// equatorial coords
$cb = cos($bMoon);
$x = $cb*cos($lMoon);
$v = $cb*sin($lMoon);
$w = sin($bMoon);
$y = $cosEPS*$v-$sinEPS*$w;
$z = $sinEPS*$v+$cosEPS*$w;
$rho = sqrt(1.0-$z*$z);
//echo "rho=".$rho." z*z=".sqrt(1.0-$z*$z)."<br>";
$this->dec = (360.0/$P2)*atan($z/$rho);
$r = (48.0/$P2)*atan($y/($x+$rho));
//echo "sinEPS=".$sinEPS."<br>";
if ($r<0) $r=$r+24;
if ($r>24) $r=$r-24;
$this->RA = $r;
}
}
class quad{
public $nz;//number of roots
public $z1;
public $z2;//the roots
public $xe; // x extreme
public $ye; // y extreme
public function __construct($yMinus, $y0, $yPlus){
$this->nz=0;
$this->z1=0;
$this->z2=0;
$c = $y0;
$a = 0.5*($yPlus+$yMinus)-$c;
$b = 0.5*($yPlus-$yMinus);
$this->xe = -1*$b/(2.0*$a);
$this->ye = ($a*$this->xe*$this->xe)+ ($b * $this->xe) + $c;
$discriminant = $b*$b - (4.0*$a*$c);
if ($discriminant >=0) { //parabola above x axis
$dx = .5*sqrt($discriminant)/abs($a);
$this->z1=$this->xe-$dx;
$this->z2=$this->xe+$dx;
if (abs($this->z1)<= 1.0) $this->nz++;
if (abs($this->z2)<= 1.0) $this->nz++;
if ($this->z1 < -1.0) $this->z1=$this->z2;
}
}
}
class moonRiseSet
{
public $risesAt="-------";
public $setsAt="-------";
public function __construct($year, $month, $day, $longitude, $latitude, $timeZone)
{
$datum=mjd($year,$month,$day)-$timeZone/24.0;
$sphi=sinD($latitude);
$cphi=cosD($latitude);
$hour=1.0;
$angSize=(8.0/60.0)*pi()/180.0;
//echo "datum=".$datum."<br>";
$yMinus = sinAltMoon( $datum, $hour-1, $longitude, $cphi, $sphi) - sin($angSize);
//echo "yMinus=".$yMinus."<br>";
$rise = false;
$sset = false;
$utRise=0;
$utSet=0;
// loop over search intervals from [0h-2h] to [22h-24h]
// go thru the day, stop if you find both rise and set
while ( ($hour < 25) && !( $rise && $sset)){
$y0=sinAltMoon( $datum, $hour, $longitude, $cphi, $sphi) - sin($angSize);
$yPlus=sinAltMoon ( $datum, $hour+1, $longitude, $cphi, $sphi) - sin($angSize);
$x = new Quad($yMinus,$y0,$yPlus);
if (($x->nz==1)&&($yMinus < 0)){
$utRise= $hour+$x->z1;
$rise = true;
}
if (($x->nz==1)&&($yMinus >=0)) {
$utSet= $hour +$x->z1;
$sset = true;
}
if (($x->nz==2) &&($x->ye <0)) {
$utRise=$hour+$x->z2;
$utSet= $hour+$x->z1;
$rise = true;
$sset = true;
}
if (($x->nz==2) && ($x->ye>=0)){
$utRise=$hour+$x->z1;
$utSet= $hour+$x->z2;
$rise = true;
$sset = true;
}
$yMinus = $yPlus;
$hour = $hour + 2;
}//end of while loop
//Now report moon, phase, sun and twilight times for the day
if ($rise) {
$this->risesAt=ut2String($utRise);
//echo "utRise=".$utRise."<br>";
}
if ($sset) {
$this->setsAt=ut2String($utSet);
}
}
}
function sinAltMoon($mjd0, $hr, $lambda, $cphi, $sphi){
$mjd=$mjd0+$hr/24.0;
$t= ($mjd-51544.5)/36525.0;
//echo "lmst=".LMST($mjd, $lambda)."<br>";
$obj=new MiniMoon($t);
$tau = 15.0 * (LMST($mjd, $lambda) - $obj->RA);
//echo "RA=".$obj->RA." dec=".$obj->dec."<br>";
return $sphi*sinD($obj->dec) + $cphi*cosD($obj->dec)*cosD($tau);
}
function LMST($mjd, $lambda) {
$mjd0 = trunc($mjd);
$ut= ($mjd-$mjd0)*24;
$t = ($mjd0-51544.5)/36525.0;
$gmst = 6.697374558 + 1.0027379093*$ut + (8640184.812866+(0.093104-0.0000062*$t)*$t)*$t/3600.0;
return 24.0*frac(($gmst+$lambda/15.0)/24.0);
}
function frac($a){
$a1 = $a - trunc($a);
if ($a1<0) $a1=$a1+1;
return $a1;
}
function trunc($t){
if ($t<0) {
return ceil($t);
}
return floor($t);
}
function div($q, $d){
return floor($q/$d);
}
function modulo($q, $d){
return ($q-$d*floor($q/$d));
}
function sinD($a){
return sin(pi()*$a/180.0);
}
function cosD($a){
return cos(pi()*$a/180.0);
}
function ut2String($ut){
$result="am";
while($ut<=0){
$ut=$ut+24.0;
}
$h=floor($ut);
$min=round( ($ut-round($ut))*60);
while ($min<=0){
$min=$min+60;
}
if ($min==60){$min=0;}
if ($h>11){
$result="pm";
}
$h=($h)%12;
if($h==0){
$h=12;
}
if ($min>9){
$result=$h.":".$min.$result;
} else {
$result=$h.":0".$min.$result;
}
return $result;
}
function mjd2String($mjd) {
$result="am";
$min=round(24*60*($mjd-floor($mjd)));
$h= div($min,60);
$min=$min % 60;
if ( $h > 11)
$result="pm";
$h = $h%12;
if ( $h==0)
$h = 12;
if ($min>9) {
$result=$h.":".$min.$result;
} else {
$result=$h.":0".$min.$result;
}
return $result;
}
function mjd($year, $month, $day)
{
$julianDay=jd($year, $month, $day,0,0,0);
return $julianDay-2400000.5;
}
function jd($year, $month, $day, $hour, $min, $sec)
{
return jdn($year, $month, $day)+($hour-12)/24.0+$min/1440.0+$sec/86400.0;
}
function jdn($year, $month, $day){
$a1=div((14-$month),12);
$y=$year+4800-$a1;
$m=$month+12*$a1-3;
return $day+div((153*$m+2),5)+365*$y+div($y,4)-div($y,100)+div($y,400)-32045;
}
function jdUnix($unix){
$sinceNoon=($unix-43200)%86400;
$jdn= 2440588+($unix-43200)/86400;
return $jdn+$sinceNoon/86400.0;
}
function mjdUnix($unix){
return jdUnix($unix)-2400000.5;
}
function jdnUnixTime($unix){
$u=$unix-43200;
$days=div($u,86400);
return 2440588+$days;
}
// diff of 1/1/1970 and new moon of 1/4/2011 9:03 UTC
// average lunation is 2551443 seconds
function moonAgeDaysUnix($unix){
return round((($unix-1294131780) % 2551443)/86400,2);
}
function moonAgeDays($year, $month, $day){
$unix=43200+86400*(mjd($year,$month,$day)-mjd(1970,1,1) );
return round((($unix-1294131780) % 2551443)/86400,2);
}
function moonPhase($mjd){
$t=($mjd-51544.5)/36525.0;
$d = 297.85020420 + 445267.1115168*$t - .00163*$t*$t + $t*$t*$t/545868 - $t*$t*$t*$t/113065000; // Moon's mean elongation
$m = 357.5291092 + 35999.0502909*$t - .0001536*$t*$t + $t*$t*$t/24490000; // Sun's mean anomaly
$m1 = 134.9634114 + 477198.8676313*$t + .0089970*$t*$t + $t*$t*$t/69699 - $t*$t*$t*$t/14712000; //Moon's mean anomaly
$i= 180 - $d - 6.289*sinD($m1) + 2.1*sinD($m) - 1.274*sinD(2*$d-$m1) - 0.658*sinD(2*$d) - .214*sinD(2*$m1) - .11*sinD($d);
return round(50*(1+cosD($i)),2);
}
echo "jdn of April 9 , 2011 is (should be 2455661)".jdn(2011,4,9)."<br>";
echo "JD of Noon April 9, 2001 is (2455661.0)".jd(2011,4,9,12,0,0)."<br>";
echo "JDN For April 9, 2011 11:59:59 AM UTC (should be 2455660)".jdnUnixTime(1302350399)."<br>";
echo "MJD of 9/27/1995 should be 49987 is ".mjd(1995,9,27)."<br>";
echo "MJD of 4/9/1995 should be 55660 is ".mjd(2011,4,9)."<br>";
echo "Age of Moon now ".moonAgeDaysUnix(time())." days<br>";
echo "Age of 4/12/2011 ".moonAgeDays(2011,4,12)." days<br>";
echo "Illumination ".moonPhase(mjdUnix(time()))." %<br>";
echo "MJD from unix time ".mjdUnix(time())."<br>";
echo "JD from unix time ".jdUnix(time())."<br>";
echo "1. MJD2 string from unix time ".mjd2String(mjdUnix(time()))."<br>";
echo "2. MJD from unix time should be 8:27 is ".mjd2String(55664.18560)."<br>";
echo "3. 4/13/2011 4:27am utc testing mjd2String : ".mjd2String(2455664.685417-2400000.5)."<br><P>";
$m = new moonRiseSet(2011,1,1,-118.2,34.2,-8);
echo "1/1 2011 Moon rise (4:36am rise, 2:37pm set) ".$m->risesAt." set ".$m->setsAt."<br>";
$m = new moonRiseSet(2011,4,13,-118.2,34.2,-8);
echo "4/13 2011 Moon rise (1:43pm rise, 2:15pm set PST) ".$m->risesAt." set ".$m->setsAt."<br>";
$m = new moonRiseSet(2011,3,24,-118.2,34.2,-8);
echo "3/24 2011 Moon rise (no rise, 9:16am set PST) ".$m->risesAt." set ".$m->setsAt."<br>";
$m = new moonRiseSet(2011,1,11,-118.2,34.2,-8);
echo "1/11 2011 Moon rise (10:43pm rise, no set PST) ".$m->risesAt." set ".$m->setsAt."<br>";
for ($i=1;$i<32;$i++){
$m = new moonRiseSet(2011,1,$i,-118.2,34.2,-8);
echo "1/".$i."/2011 ".$m->risesAt." ".$m->setsAt."<br>";
}
?>