(item C.11) describes John Horton Conway's approximation for the moon phase of a given date. Here it is in Tcl: proc moonphase date {
set y [clock format $date -format %Y]
set delta [expr {$y/100==19? 4: 8.3}]
scan [clock format $date -format %m] %d m
set m [expr {$m==1? 3: $m==2? 4: $m}]
scan [clock format $date -format %d] %d d
set t [expr ($y%100)%19]
if {$t>9} {incr t -19}
expr {(round(($t*11)%30+$m+$d-$delta)%30)}
}--- Testing:% moonphase [clock scan "June 6, 1944"] 14which matches the FAQ's example and stands for (almost) full moon.
Fred Limouzin - 2005-03-21, updated on 2005-03-24: Here's a little week-end project related to the above. It's still a draft, but it's getting there! It simply displays the moon phase in a canvas, based on a percentage value representing the progression within a full moon cycle. 0% is a new Moon, 50% a full Moon, 100% is the next new Moon (in fact 100% is clamped back to 0, etc.).I was going to create a new page, but after a quick search on wiki, I found out about this page, so I'm gonna append to it!Here's a screeshot:
You'll also need the background image (in fact the most important part!): (I haven't decided yet if the image should be rotated. Last night was overcast ;-))
The code and images can also be found at http://dire.straits.free.fr/vertigo ([1]). #!/bin/sh
# Frederic Limouzin ; Copyrights (c)2005 - all rights reserved \
exec tclsh "$0" ${1+"$@"}
package require Tk
set ::OFFSET 99
set ::RADIUS 100
set ::STEPY 2
set ::pi [expr {acos(-1.0)}]
####################################################
# ellipse-shaped Mask
# x^2/a^2 + y^2/b^2 = 1
####################################################
proc ElipseMask {rx1 rx2 ry clr} {
if {$ry == 0} { return 0 }
set lst [list]
foreach rx [list $rx1 $rx2] rysign {1.0 -1.0} {
if {$rx < 0} {
set rxsign -1.0
set rx [expr {abs($rx)}]
} else {
set rxsign 1.0
}
for {set y -$ry} {$y <= $ry} {incr y $::STEPY} {
set t [expr {1.0 - ((1.0 * $y * $y) / (1.0 * $ry * $ry))}]
set x [expr {round(1.0 * $rx * sqrt($t))}]
lappend lst [expr {round(($rxsign * $x) + $::OFFSET)}]
lappend lst [expr {round(($rysign * $y) + $::OFFSET)}]
}
}
return [.c create polygon $lst -fill $clr -outline $clr]
}
####################################################
# p in percent of a full cycle within [0.0;100.0]
# 0% and 100% : new moon; 50% full moon; etc.
####################################################
# This gives an approximation for a flat Moon rather
# than spherical Moon. But after all it is well known
# that Earth is flat, so the Moon must be too...
####################################################
proc ElipsePhase {p} {
while {$p >= 100.0} {
set p [expr {1.0 * ($p - 100.0)}]
}
set phasesLst { {New Moon} {Waxing Crescent} {First Quarter} \
{Waxing Gibbous} {Full Moon} {Waning Gibbous} \
{Last Quarter} {Waning Crescent} }
set idx [expr {int(8.0 * ($p + (100.0/16.0)) / 100.0) % 8}]
set phase [lindex $phasesLst $idx]
set quadrant [expr {int(1.0 * $p / 25.0)}]
set mod [expr {(1.0 * $p) - (25.0 * $quadrant)}]
if {$quadrant == 3} {
set rx1 [expr {-1.0 * (4.0 * $mod)}]
} elseif {$quadrant == 2} {
set rx1 [expr {100.0 - (4.0 * $mod)}]
} else {
set rx1 -$::RADIUS
}
if {$quadrant == 0} {
set rx2 [expr {100.0 - (4.0 * $mod)}]
} elseif {$quadrant == 1} {
set rx2 [expr {-1.0 * (4.0 * $mod)}]
} else {
set rx2 $::RADIUS
}
return [list $rx1 $rx2 $phase]
}
# normalize in [0.0;1.0]
proc normalize {{v 0.0}} {
set v [expr {1.0 * ($v - floor($v))}]
if {$v < 0.0} {
set v [expr {1.0 + $v}]
}
return $v
}
# based on "Lunar Phase Computation" by Stephen R. Schmitt
# based on "Sky & Telescope, Astronomical Computing", April 1994
proc ComputeLunarPhase {{tmr {}}} {
if {$tmr eq {}} {
set tmr [clock seconds] ;# today by default
}
set year [clock format $tmr -format %Y]
scan [clock format $tmr -format %m] %d month
scan [clock format $tmr -format %d] %d day
set date [format {%4d/%02d/%02d} $year $month $day]
# Julian date at 12h UT
set JulianYear [expr {$year - int((12.0 - $month) / 10.0)}]
set JulianMonth [expr {($month + 9) % 12}]
set K1 [expr {int(365.25 * ($JulianYear + 4712.0))}]
set K2 [expr {int((30.6 * $JulianMonth) + 0.5)}]
set K3 [expr {int(floor(($JulianYear / 100.0) + 49.0) * 3 / 4.0) - 38}]
set JulianDay [expr {$K1 + $K2 + $day + 59}]
set GregorianDay [expr {$JulianDay - $K3}]
set JGDay [expr {($JulianDay > 2299160) ? $GregorianDay : $JulianDay}]
#calculate moon's age in days
set MaxAge 29.530588853
set IP [normalize [expr {($JGDay - 2451550.1) / $MaxAge}]]
set MoonAge [expr {$IP * $MaxAge}]
set pc [expr {(100.0 * $MoonAge) / $MaxAge}]
set IPrad [expr {$IP * $::pi}] ;# radians
#calculate moon's distance
set tmp [normalize [expr {($JGDay - 2451562.2 ) / 27.55454988}]]
set DP [expr {2.0 * $::pi * $tmp}]
set dist [expr {60.4 - 3.3*cos($DP) - 0.6*cos((2*$IPrad) - $DP)}]
set dist [expr {$dist - 0.5*cos(2*$IPrad)}]
#calculate moon's ecliptic latitude
set tmp [normalize [expr {($JGDay - 2451565.2) / 27.212220817}]]
set NP [expr {2.0 * $::pi * $tmp}]
set EclipticLatitude [expr {5.1 * sin($NP)}]
# calculate moon's ecliptic longitude
set RP [normalize [expr {($JGDay - 2451555.8) / 27.321582241}]]
set tmp [expr {(360.0 * $RP) + (6.3*sin($DP))}]
set tmp [expr {$tmp + (1.3*sin(2.0*$IPrad - $DP))}]
set EclipticLongitude [expr {$tmp + (0.7*sin(2.0*$IPrad))}]
if {$EclipticLongitude < 33.18} { set Zodiac "Pisces"
} elseif {$EclipticLongitude < 51.16} { set Zodiac "Aries"
} elseif {$EclipticLongitude < 93.44} { set Zodiac "Taurus"
} elseif {$EclipticLongitude < 119.48} { set Zodiac "Gemini"
} elseif {$EclipticLongitude < 135.30} { set Zodiac "Cancer"
} elseif {$EclipticLongitude < 173.34} { set Zodiac "Leo"
} elseif {$EclipticLongitude < 224.17} { set Zodiac "Virgo"
} elseif {$EclipticLongitude < 242.57} { set Zodiac "Libra"
} elseif {$EclipticLongitude < 271.26} { set Zodiac "Scorpio"
} elseif {$EclipticLongitude < 302.49} { set Zodiac "Sagittarius"
} elseif {$EclipticLongitude < 311.72} { set Zodiac "Capricorn"
} elseif {$EclipticLongitude < 348.58} { set Zodiac "Aquarius"
} else { set Zodiac "Pisces"
}
return [list $pc $date $dist $EclipticLatitude $EclipticLongitude $Zodiac]
}
####################################################
proc UpdateMoon args {
global obj
global phaselbl
global infotxt
global zodiac
set ry $::RADIUS
set secinaday [expr {24 * 60 * 60}]
set co #080808
catch {.c delete $obj} res
set dayoffset [.s get]
set tmr [expr {[clock second] + ($dayoffset * $secinaday)}]
;#;#test: 28 fev 2004 new moon, 6 march 2004 full, 13 march last quarter
;#set tmr [expr {[clock scan {March 1, 2004}] + ($dayoffset * $secinaday)}]
foreach {pc date -- eclat eclon zod} [ComputeLunarPhase $tmr] {break;#}
set infotxt [format {(Ecliptic Lat=%f°;Lon=%f°)} $eclat $eclon]
set zodiac [format {%s: Constellation : %s} $date $zod]
foreach {rx1 rx2 phaselbl} [ElipsePhase $pc] {break;#} ;#assign
set obj [ElipseMask $rx1 $rx2 $ry $co]
update
}
####################################################
set phaselbl {}
set zodiac {}
set infotxt {}
canvas .c -width [expr {2 * $::RADIUS}] \
-height [expr {2 * $::RADIUS}] -background black
label .l -textvariable phaselbl
label .i -textvariable infotxt
label .z -textvariable zodiac
scale .s -from -30 -to 30 -length 300 -resolution 1 \
-label "Day offset from Today" -variable dayoffset -command {UpdateMoon} \
-orient horizontal -tickinterval 4 -showvalue true
pack .c -side top
pack .s -side bottom
pack .i -side bottom
pack .z -side bottom
pack .l -side bottom
# Full moon image as the 'background'
set fname [file join [file dirname [info script]] moon.gif]
set img [image create photo -file $fname]
.c create image 0 0 -image $img -anchor nw
set ry $::RADIUS
set co #111111
####################################################
.s set 0
UpdateMoonCheck the code for references of the 'ComputeLunarPhase' procedure. The rest is all mine :-). Note that the mask is only an approximation for a flat Moon rather than a spherical Moon. But after all it is well know that Earth is flat, therefore the Moon must be too... I haven't done much testing, the result may be 1 day behind or ahead of schedule. Yet again, it could be due to the flat approximation. I'll need to add some Sine function in the ElipsePhase procedure at some stage!Cheers,--Fred
Arts and crafts of tcl-Tk programming - Category Science
