Skip to content

Instantly share code, notes, and snippets.

@dpoggi
Last active August 29, 2015 14:17
Show Gist options
  • Save dpoggi/c7704a3b6bf9b47171e6 to your computer and use it in GitHub Desktop.
Save dpoggi/c7704a3b6bf9b47171e6 to your computer and use it in GitHub Desktop.
AstroBuild Go port
package main
import (
"math"
"time"
"fmt"
"os"
)
var planets map[float64][]string
type Planet struct {
N, i, w, a, e, M float64
}
type OrbitalElements struct {
v, r, xh, yh, w, lonecl, latecl float64
}
func radians(degrees float64) float64 {
return degrees * (math.Pi / 180.0)
}
func degrees(radians float64) float64 {
return radians / (math.Pi / 180.0)
}
func getday(t time.Time) float64 {
year := t.Year()
month := int(t.Month())
day := t.Day()
hour := t.Hour()
return float64((367 * year) - (7 * (year + (month + 9) / 12) / 4) + (275 * month / 9) + day - 730530) + (float64(hour) / 24.0)
}
func getplanet(planet string, day float64) Planet {
var p Planet
switch planet {
case "Sun":
p.N = radians(0.0)
p.i = radians(0.0)
p.w = radians(282.9404 + 4.70935e-5 * day)
p.a = 1.000000
p.e = 0.016709 - 1.151e-9 * day
p.M = radians(356.0470 + 0.9856002585 * day)
case "Mercury":
p.N = radians(48.3313 + 3.24587e-5 * day)
p.i = radians(7.0047 + 5.00e-8 * day)
p.w = radians(29.1241 + 1.01444e-5 * day)
p.a = 0.387098
p.e = 0.205635 + 5.59e-10 * day
p.M = radians(168.6562 + 4.0923344368 * day)
case "Venus":
p.N = radians(76.6799 + 2.46590e-5 * day)
p.i = radians(3.3946 + 2.75e-8 * day)
p.w = radians(54.8910 + 1.38374e-5 * day)
p.a = 0.723330
p.e = 0.006773 - 1.302e-9 * day
p.M = radians(48.0052 + 1.6021302244 * day)
case "Mars":
p.N = radians(49.5574 + 2.11081e-5 * day)
p.i = radians(1.8497 - 1.78e-8 * day)
p.w = radians(286.5016 + 2.92961e-5 * day)
p.a = 1.523688
p.e = 0.093405 + 2.516e-9 * day
p.M = radians(18.6021 + 0.5240207766 * day)
case "Jupiter":
p.N = radians(100.4542 + 2.76854e-5 * day)
p.i = radians(1.3030 - 1.557e-7 * day)
p.w = radians(273.8777 + 1.64505e-5 * day)
p.a = 5.20256
p.e = 0.048498 + 4.469e-9 * day
p.M = radians(19.8950 + 0.0830853001 * day)
case "Saturn":
p.N = radians(113.6634 + 2.38980e-5 * day)
p.i = radians(2.4886 - 1.081e-7 * day)
p.w = radians(339.3939 + 2.97661e-5 * day)
p.a = 9.55475
p.e = 0.055546 - 9.499e-9 * day
p.M = radians(316.9670 + 0.0334442282 * day)
case "Uranus":
p.N = radians(74.0005 + 1.3978e-5 * day)
p.i = radians( 0.7733 + 1.9e-8 * day)
p.w = radians(96.6612 + 3.0565e-5 * day)
p.a = 19.18171 - 1.55e-8 * day
p.e = 0.047318 + 7.45e-9 * day
p.M = radians(142.5905 + 0.011725806 * day)
case "Neptune":
p.N = radians(131.7806 + 3.0173e-5 * day)
p.i = radians(1.7700 - 2.55e-7 * day)
p.w = radians(272.8461 - 6.027e-6 * day)
p.a = 30.05826 + 3.313e-8 * day
p.e = 0.008606 + 2.15e-9 * day
p.M = radians(260.2471 + 0.005995147 * day)
}
return p
}
func getorbitalelements(p Planet) OrbitalElements {
var o OrbitalElements
E := p.M + p.e * math.Sin(p.M) * (1.0 + p.e * math.Cos(p.M))
xv := p.a * (math.Cos(E) - p.e)
yv := p.a * (math.Sqrt(1.0 - p.e * p.e) * math.Sin(E))
o.v = math.Atan2(yv, xv)
o.r = math.Sqrt(xv * xv + yv * yv)
o.xh = o.r * (math.Cos(p.N) * math.Cos(o.v + p.w) - math.Sin(p.N) * math.Sin(o.v + p.w) * math.Cos(p.i))
o.yh = o.r * (math.Sin(p.N) * math.Cos(o.v + p.w) + math.Cos(p.N) * math.Sin(o.v + p.w) * math.Cos(p.i))
o.w = p.w
zh := o.r * (math.Sin(o.v + p.w) * math.Sin(p.i))
o.lonecl = math.Atan2(o.yh, o.xh)
o.latecl = math.Atan2(zh, math.Sqrt(o.xh * o.xh + o.yh * o.yh))
return o
}
func getgeocentricalignment(planet string, day float64) float64 {
s := getorbitalelements(getplanet("Sun", day))
p := getorbitalelements(getplanet(planet, day))
lonsun := s.v + s.w
xs := s.r * math.Cos(lonsun)
ys := s.r * math.Sin(lonsun)
xg := p.xh + xs
yg := p.yh + ys
degrees := 90 - degrees(math.Atan2(xg, yg))
if degrees < 0 {
degrees += 360
}
return degrees
}
func init() {
planetnames := []string{"Sun", "Mercury", "Venus", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"}
day := getday(time.Now().UTC())
planets = make(map[float64][]string)
var alignment float64
for _, planet := range planetnames {
alignment = math.Floor(getgeocentricalignment(planet, day) + 0.5)
if _, ok := planets[alignment]; ok {
planets[alignment] = append(planets[alignment], planet)
} else {
planets[alignment] = []string{planet}
}
}
}
func main() {
flag := 1
for alignment, bodies := range planets {
if len(bodies) > 1 {
fmt.Printf("%fdeg\t", alignment)
for _, body := range bodies {
fmt.Printf("\t%s", body)
}
fmt.Printf("\n")
flag = 0
}
}
os.Exit(flag)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment