Update from main

conversions.go

@@ -20,7 +20,7 @@ func days2mdhms(year int64, epochDays float64) (mon, day, hr, min, sec float64)
 
 	for dayofyr > inttemp+float64(lmonth[int(i-1)]) && i < 22 {
 		inttemp = inttemp + float64(lmonth[int(i-1)])
-		i += 1
+		i++
 	}
 
 	mon = i
@@ -37,8 +37,8 @@ func days2mdhms(year int64, epochDays float64) (mon, day, hr, min, sec float64)
 	return
 }
 
-// Calc julian date given year, month, day, hour, minute and second
-// the julian date is defined by each elapsed day since noon, jan 1, 4713 bc.
+// JDay calculates julian date given year, month, day, hour, minute and second.
+// The julian date is defined by each elapsed day since noon, jan 1, 4713 bc. No leap second is taken into account.
 func JDay(year, mon, day, hr, min, sec int) float64 {
 	return (367.0*float64(year) - math.Floor((7*(float64(year)+math.Floor((float64(mon)+9)/12.0)))*0.25) + math.Floor(275*float64(mon)/9.0) + float64(day) + 1721013.5 + ((float64(sec)/60.0+float64(min))/60.0+float64(hr))/24.0)
 }
@@ -56,13 +56,13 @@ func gstime(jdut1 float64) (temp float64) {
 	return
 }
 
-// Calc GST given year, month, day, hour, minute and second
+// GSTimeFromDate computes GST given year, month, day, hour, minute and second
 func GSTimeFromDate(year, mon, day, hr, min, sec int) float64 {
 	jDay := JDay(year, mon, day, hr, min, sec)
 	return gstime(jDay)
 }
 
-// Convert Earth Centered Inertial coordinated into equivalent latitude, longitude, altitude and velocity.
+// ECIToLLA converts Earth Centered Inertial coordinated into equivalent latitude, longitude, altitude and velocity.
 // Reference: http://celestrak.com/columns/v02n03/
 func ECIToLLA(eciCoords Vector3, gmst float64) (altitude, velocity float64, ret LatLong) {
 	a := 6378.137     // Semi-major Axis
@@ -95,7 +95,7 @@ func ECIToLLA(eciCoords Vector3, gmst float64) (altitude, velocity float64, ret
 	return
 }
 
-// Convert LatLong in radians to LatLong in degrees
+// LatLongDeg converts LatLong in radians to LatLong in degrees
 func LatLongDeg(rad LatLong) (deg LatLong) {
 	deg.Longitude = math.Mod(rad.Longitude/math.Pi*180, 360)
 	if deg.Longitude > 180 {
@@ -111,10 +111,10 @@ func LatLongDeg(rad LatLong) (deg LatLong) {
 	return
 }
 
-// Calculate GMST from Julian date.
+// ThetaGJD calculates GMST from Julian date.
 // Reference: The 1992 Astronomical Almanac, page B6.
-func ThetaG_JD(jday float64) (ret float64) {
-	_, UT := math.Modf(jday + 0.5)
+func ThetaGJD(jday float64) (ret float64) {
+	UT := jday + 0.5 - jday
 	jday = jday - UT
 	TU := (jday - 2451545.0) / 36525.0
 	GMST := 24110.54841 + TU*(8640184.812866+TU*(0.093104-TU*6.2e-6))
@@ -123,19 +123,19 @@ func ThetaG_JD(jday float64) (ret float64) {
 	return
 }
 
-// Convert latitude, longitude and altitude(km) into equivalent Earth Centered Intertial coordinates(km)
+// LLAToECI converts latitude, longitude and altitude into equivalent Earth Centered Intertial coordinates.
 // Reference: The 1992 Astronomical Almanac, page K11.
 func LLAToECI(obsCoords LatLong, alt, jday float64) (eciObs Vector3) {
 	re := 6378.137
-	theta := math.Mod(ThetaG_JD(jday)+obsCoords.Longitude, TWOPI)
+	theta := math.Mod(ThetaGJD(jday)+obsCoords.Longitude, TWOPI)
 	r := (re + alt) * math.Cos(obsCoords.Latitude)
 	eciObs.X = r * math.Cos(theta)
 	eciObs.Y = r * math.Sin(theta)
 	eciObs.Z = (re + alt) * math.Sin(obsCoords.Latitude)
 	return
 }
 
-// Convert Earth Centered Intertial coordinates into Earth Cenetered Earth Final coordinates
+// ECIToECEF converts Earth Centered Intertial coordinates into Earth Cenetered Earth Final coordinates.
 // Reference: http://ccar.colorado.edu/ASEN5070/handouts/coordsys.doc
 func ECIToECEF(eciCoords Vector3, gmst float64) (ecfCoords Vector3) {
 	ecfCoords.X = eciCoords.X*math.Cos(gmst) + eciCoords.Y*math.Sin(gmst)
@@ -144,12 +144,11 @@ func ECIToECEF(eciCoords Vector3, gmst float64) (ecfCoords Vector3) {
 	return
 }
 
-// Calculate look angles for given satellite position and observer position
+// ECIToLookAngles calculates look angles for given satellite position and observer position
 // obsAlt in km
 // Reference: http://celestrak.com/columns/v02n02/
 func ECIToLookAngles(eciSat Vector3, obsCoords LatLong, obsAlt, jday float64) (lookAngles LookAngles) {
-	theta := math.Mod(ThetaG_JD(jday)+obsCoords.Longitude, 2*math.Pi)
-	obsPos := LLAToECI(obsCoords, obsAlt, jday)
+	gst := ThetaGJD(jday)
 
 	rx := eciSat.X - obsPos.X
 	ry := eciSat.Y - obsPos.Y
@@ -159,15 +158,17 @@ func ECIToLookAngles(eciSat Vector3, obsCoords LatLong, obsAlt, jday float64) (l
 	top_e := -math.Sin(theta)*rx + math.Cos(theta)*ry
 	top_z := math.Cos(obsCoords.Latitude)*math.Cos(theta)*rx + math.Cos(obsCoords.Latitude)*math.Sin(theta)*ry + math.Sin(obsCoords.Latitude)*rz
 
-	lookAngles.Az = math.Atan(-top_e / top_s)
-	if top_s > 0 {
-		lookAngles.Az = lookAngles.Az + math.Pi
-	}
-	if lookAngles.Az < 0 {
-		lookAngles.Az = lookAngles.Az + 2*math.Pi
-	}
-	lookAngles.Rg = math.Sqrt(rx*rx + ry*ry + rz*rz)
-	lookAngles.El = math.Asin(top_z / lookAngles.Rg)
+	rx := ecefSat.X - ecefObs.X
+	ry := ecefSat.Y - ecefObs.Y
+	rz := ecefSat.Z - ecefObs.Z
+
+	topS := math.Sin(obsCoords.Latitude)*math.Cos(obsCoords.Longitude)*rx + math.Sin(obsCoords.Latitude)*math.Sin(obsCoords.Longitude)*ry - math.Cos(obsCoords.Latitude)*rz
+	topE := -math.Sin(obsCoords.Longitude)*rx + math.Cos(obsCoords.Longitude)*ry
+	topZ := math.Cos(obsCoords.Latitude)*math.Cos(obsCoords.Longitude)*rx + math.Cos(obsCoords.Latitude)*math.Sin(obsCoords.Longitude)*ry + math.Sin(obsCoords.Latitude)*rz
+
+	lookAngles.Rg = math.Sqrt(topS*topS + topE*topE + topZ*topZ)
+	lookAngles.Az = math.Atan2(-topE, topS) + math.Pi
+	lookAngles.El = math.Asin(topZ / lookAngles.Rg)
 
 	return
 }

dspace.go

@@ -4,12 +4,12 @@ import (
 	"math"
 )
 
-// A struct returned from the dsinit function
+// DeepSpaceInitResult is the struct returned from the dsinit function
 type DeepSpaceInitResult struct {
 	em, argpm, inclm, mm, nm, nodem, irez, atime, d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433, dedt, didt, dmdt, dndt, dnodt, domdt, del1, del2, del3, xfact, xlamo, xli, xni float64
 }
 
-// A struct returned from the dspace function
+// DeepSpaceResult is the struct returned from the dspace function
 type DeepSpaceResult struct {
 	atime, em, argpm, inclm, xli, mm, xni, nodem, dndt, nm float64
 }
@@ -314,7 +314,7 @@ func dspace(irez, d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421,
 	return
 }
 
-// A struct returned from the dpper function
+// DpperResult is the struct returned from the dpper function
 type DpperResult struct {
 	ep, inclp, nodep, argpp, mp float64
 }
@@ -454,7 +454,7 @@ func dpper(satrec *Satellite, inclo float64, init string, ep, inclp, nodep, argp
 	return
 }
 
-// A struct returned from the dscom function
+// DSComResults is the struct returned from the dscom function
 type DSComResults struct {
 	snodm, cnodm, sinim, cosim, sinomm, cosomm, day, e3, ee2, em, emsq, gam, peo, pgho, pho, pinco, plo, rtemsq, se2, se3, sgh2, sgh3, sgh4, sh2, sh3, si2, si3, sl2, sl3, sl4, s1, s2, s3, s4, s5, s6, s7, ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3, sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33, xgh2, xgh3, xgh4, xh2, xh3, xi2, xi3, xl2, xl3, xl4, nm, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33, zmol, zmos float64
 }

gravity.go

@@ -5,12 +5,12 @@ import (
 	"math"
 )
 
-// Holds variables that are dependent upon selected gravity model
+// GravConst holds variables that are dependent upon selected gravity model.
 type GravConst struct {
 	mu, radiusearthkm, xke, tumin, j2, j3, j4, j3oj2 float64
 }
 
-// Returns a GravConst with correct information on requested model provided through the name parameter
+// Returns a GravConst with correct information on requested model provided through the name parameter.
 func getGravConst(name string) (grav GravConst) {
 	switch name {
 	case "wgs72old":

helpers.go

@@ -8,27 +8,33 @@ import (
 )
 
 // Constants
-const TWOPI float64 = math.Pi * 2.0
-const DEG2RAD float64 = math.Pi / 180.0
-const RAD2DEG float64 = 180.0 / math.Pi
-const XPDOTP float64 = 1440.0 / (2.0 * math.Pi)
+const (
+	// Pi*2
+	TWOPI float64 = math.Pi * 2.0
+	// Degrees to radians
+	DEG2RAD float64 = math.Pi / 180.0
+	// Radians to degrees
+	RAD2DEG float64 = 180.0 / math.Pi
+	// Xp'p
+	XPDOTP float64 = 1440.0 / (2.0 * math.Pi)
+)
 
-// Holds latitude and Longitude in either degrees or radians
+// LatLong holds latitude and Longitude in either degrees or radians
 type LatLong struct {
 	Latitude, Longitude float64
 }
 
-// Holds X, Y, Z position
+// Vector3 holds X, Y, Z position
 type Vector3 struct {
 	X, Y, Z float64
 }
 
-// Holds an azimuth, elevation and range
+// LookAngles holds an azimuth, elevation and range
 type LookAngles struct {
 	Az, El, Rg float64
 }
 
-// Parses a two line element dataset into a Satellite struct
+// ParseTLE parses a two line element dataset into a Satellite struct
 func ParseTLE(line1, line2, gravconst string) (sat Satellite) {
 	sat.Line1 = line1
 	sat.Line2 = line2
@@ -58,7 +64,7 @@ func ParseTLE(line1, line2, gravconst string) (sat Satellite) {
 	return
 }
 
-// Converts a two line element data set into a Satellite struct and runs sgp4init
+// TLEToSat converts a two line element data set into a Satellite struct and runs sgp4init
 func TLEToSat(line1, line2 string, gravconst string) Satellite {
 	//sat := Satellite{Line1: line1, Line2: line2}
 	sat := ParseTLE(line1, line2, gravconst)
@@ -74,7 +80,7 @@ func TLEToSat(line1, line2 string, gravconst string) Satellite {
 	sat.argpo = sat.argpo * DEG2RAD
 	sat.mo = sat.mo * DEG2RAD
 
-	var year int64 = 0
+	var year int64
 	if sat.epochyr < 57 {
 		year = sat.epochyr + 2000
 	} else {

satellite.go

@@ -1,6 +1,6 @@
 package satellite
 
-// Struct for holding satellite information during and before propagation
+// Satellite is the struct for holding satellite information during and before propagation
 type Satellite struct {
 	Line1 string
 	Line2 string

satellite_suite_test.go

@@ -1,6 +1,8 @@
 package satellite
 
 import (
+	"time"
+
 	. "github.com/onsi/ginkgo"
 	. "github.com/onsi/gomega"
 
@@ -221,7 +223,7 @@ var _ = Describe("go-satellite", func() {
 			PropagationTestCase{
 				line1: "1 23599U 95029B   06171.76535463  .00085586  12891-6  12956-2 0  2905",
 				line2: "2 23599   6.9327   0.2849 5782022 274.4436  25.2425  4.47796565123555",
-				grav: "wgs72",
+				grav:  "wgs72",
 				testData: `0.00000000 9892.63794341 35.76144969 -1.08228838 3.556643237 6.456009375 0.783610890       
 20.00000000 11931.95642997 7340.74973750 886.46365987 0.308329116 5.532328972 0.672887281
 40.00000000 11321.71039205 13222.84749156 1602.40119049 -1.151973982 4.285810871 0.521919425
@@ -296,4 +298,19 @@ func propagationTest(testCase PropagationTestCase) {
 			})
 		})
 	}
-}
+}
+
+func TestPropagationDate(t *testing.T) {
+	line1 := "1 25730U 99025A   18054.54151885  .00000600  00000-0  33327-3 0  9994"
+	line2 := "2 25730  99.0366  41.4381 0012822 245.3407 236.6366 14.15015322967876"
+
+	sat := TLEToSat(line1, line2, "wgs84")
+	pos, vel := Propagate(sat, 2018, 01, 10, 20, 30, 40)
+
+	date := time.Date(2018, 1, 10, 20, 30, 40, 0, time.UTC)
+	pos2, vel2 := sat.PropagateDate(&date)
+
+	if pos.X != pos2.X || pos.Y != pos2.Y || pos.Z != pos2.Z || vel.X != vel2.X || vel.Y != vel2.Y || vel.Z != vel2.Z {
+		t.Fatal("Returned values between Propagate and PropagateDate differ.")
+	}
+}

sgp4.go

@@ -2,6 +2,7 @@ package satellite
 
 import (
 	"math"
+	"time"
 )
 
 // this procedure initializes variables for sgp4.
@@ -258,10 +259,10 @@ func initl(satn int64, grav GravConst, ecco, epoch, inclo, noIn float64, methodI
 
 	ak = math.Pow(grav.xke/noIn, x2o3)
 	d1 = 0.75 * grav.j2 * (3.0*cosio2 - 1.0) / (rteosq * omeosq)
-	del_ := d1 / (ak * ak)
-	adel = ak * (1.0 - del_*del_ - del_*(1.0/3.0+134.0*del_*del_/81.0))
-	del_ = d1 / (adel * adel)
-	no = noIn / (1.0 + del_)
+	del := d1 / (ak * ak)
+	adel = ak * (1.0 - del*del - del*(1.0/3.0+134.0*del*del/81.0))
+	del = d1 / (adel * adel)
+	no = noIn / (1.0 + del)
 
 	ao = math.Pow(grav.xke/no, x2o3)
 	sinio = math.Sin(inclo)
@@ -291,13 +292,22 @@ func initl(satn int64, grav GravConst, ecco, epoch, inclo, noIn float64, methodI
 	return
 }
 
-// Calculates position and velocity vectors for given time
+// Propagate computes position and velocity vectors for given time.
 func Propagate(sat Satellite, year int, month int, day, hours, minutes, seconds int) (position, velocity Vector3) {
 	j := JDay(year, month, day, hours, minutes, seconds)
 	m := (j - sat.jdsatepoch) * 1440
 	return sgp4(&sat, m)
 }
 
+// PropagateDate computes position and velocity vectors for given time.Time.
+func (sat *Satellite) PropagateDate(date *time.Time) (position, velocity Vector3) {
+	date2 := date.UTC()
+	j := JDay(date2.Year(), int(date2.Month()), date2.Day(), date2.Hour(), date2.Minute(), date2.Second())
+	j += float64(date2.Nanosecond()) / 1.e9 / 86400.
+	m := (j - sat.jdsatepoch) * 1440
+	return sgp4(sat, m)
+}
+
 // this procedure is the sgp4 prediction model from space command. this is an updated and combined version of sgp4 and sdp4, which were originally published separately in spacetrack report #3. this version follows the methodology from the aiaa paper (2006) describing the history and development of the code.
 // satrec - initialized Satellite struct from sgp4init
 // tsince - time since epoch in minutes