From 056c8caf1719846026922d8979afdadba41aed29 Mon Sep 17 00:00:00 2001 From: siritinga Date: Sat, 24 Feb 2018 10:57:06 +0100 Subject: [PATCH 1/2] Fixed some comments and variable names. --- conversions.go | 44 ++++++++++++++++++++++---------------------- dspace.go | 8 ++++---- gravity.go | 4 ++-- helpers.go | 26 ++++++++++++++++---------- satellite.go | 2 +- sgp4.go | 10 +++++----- 6 files changed, 50 insertions(+), 44 deletions(-) diff --git a/conversions.go b/conversions.go index 1ea2213..b2171fc 100644 --- a/conversions.go +++ b/conversions.go @@ -9,7 +9,7 @@ import ( func days2mdhms(year int64, epochDays float64) (mon, day, hr, min, sec float64) { lmonth := [12]int{31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31} - if year % 4 == 0 { + if year%4 == 0 { lmonth = [12]int{31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31} } @@ -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,9 +111,9 @@ 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) { +func ThetaGJD(jday float64) (ret float64) { UT := jday + 0.5 - jday jday = jday - UT TU := (jday - 2451545.0) / 36525.0 @@ -123,11 +123,11 @@ func ThetaG_JD(jday float64) (ret float64) { return } -// Convert latitude, longitude and altitude into equivalent Earth Centered Intertial coordinates +// 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) @@ -135,20 +135,20 @@ func LLAToECI(obsCoords LatLong, alt, jday float64) (eciObs Vector3) { 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) - ecfCoords.Y = eciCoords.X *-math.Sin(gmst) + eciCoords.Y * math.Cos(gmst) + ecfCoords.X = eciCoords.X*math.Cos(gmst) + eciCoords.Y*math.Sin(gmst) + ecfCoords.Y = eciCoords.X*-math.Sin(gmst) + eciCoords.Y*math.Cos(gmst) ecfCoords.Z = eciCoords.Z 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) { - gst := ThetaG_JD(jday) + gst := ThetaGJD(jday) // Convert eciSat to ecf ecefSat := ECIToECEF(eciSat, gst) @@ -161,13 +161,13 @@ func ECIToLookAngles(eciSat Vector3, obsCoords LatLong, obsAlt, jday float64) (l ry := ecefSat.Y - ecefObs.Y rz := ecefSat.Z - ecefObs.Z - top_s := math.Sin(obsCoords.Latitude) * math.Cos(obsCoords.Longitude) * rx + math.Sin(obsCoords.Latitude) * math.Sin(obsCoords.Longitude) * ry - math.Cos(obsCoords.Latitude) * rz - top_e := -math.Sin(obsCoords.Longitude) * rx + math.Cos(obsCoords.Longitude) * ry - top_z := math.Cos(obsCoords.Latitude) * math.Cos(obsCoords.Longitude) * rx + math.Cos(obsCoords.Latitude) * math.Sin(obsCoords.Longitude) * ry + math.Sin(obsCoords.Latitude) * rz + 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(top_s*top_s + top_e*top_e + top_z*top_z) - lookAngles.Az = math.Atan2(-top_e, top_s) + math.Pi - lookAngles.El = math.Asin(top_z / lookAngles.Rg) + 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 } diff --git a/dspace.go b/dspace.go index 8a73fc3..454cf8c 100644 --- a/dspace.go +++ b/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 } diff --git a/gravity.go b/gravity.go index e87b33c..d6f5295 100644 --- a/gravity.go +++ b/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": diff --git a/helpers.go b/helpers.go index 504656f..568aa3e 100644 --- a/helpers.go +++ b/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 { diff --git a/satellite.go b/satellite.go index f0e17c8..3cd9720 100644 --- a/satellite.go +++ b/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 diff --git a/sgp4.go b/sgp4.go index d7cbca8..9e1afbb 100644 --- a/sgp4.go +++ b/sgp4.go @@ -258,10 +258,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,7 +291,7 @@ 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 From a8b25910a81b340d605af2e424c2601ae04fc21f Mon Sep 17 00:00:00 2001 From: siritinga Date: Sat, 24 Feb 2018 20:13:56 +0100 Subject: [PATCH 2/2] Additional func to propagate using time.Time. It also works with fractions of second. --- satellite_suite_test.go | 21 +++++++++++++++++++-- sgp4.go | 10 ++++++++++ 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/satellite_suite_test.go b/satellite_suite_test.go index cd0c6fa..cee32a3 100644 --- a/satellite_suite_test.go +++ b/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) { }) }) } -} \ No newline at end of file +} + +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.") + } +} diff --git a/sgp4.go b/sgp4.go index 9e1afbb..924c32d 100644 --- a/sgp4.go +++ b/sgp4.go @@ -2,6 +2,7 @@ package satellite import ( "math" + "time" ) // this procedure initializes variables for sgp4. @@ -298,6 +299,15 @@ func Propagate(sat Satellite, year int, month int, day, hours, minutes, seconds 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