diff --git a/README.md b/README.md index 9679458..cd4f3b6 100644 --- a/README.md +++ b/README.md @@ -12,19 +12,15 @@ I decided to port the SGP4 library to GoLang as one of my first projects with th #### Constants ```go -const DEG2RAD float64 = math.Pi / 180.0 +const deg2rad float64 = math.Pi / 180.0 ``` ```go -const RAD2DEG float64 = 180.0 / math.Pi +const twopi float64 = math.Pi * 2.0 ``` ```go -const TWOPI float64 = math.Pi * 2.0 -``` - -```go -const XPDOTP float64 = 1440.0 / (2.0 * math.Pi) +const xPDOTP float64 = 1440.0 / (2.0 * math.Pi) ``` #### func ECIToLLA @@ -57,10 +53,10 @@ func Propagate(sat Satellite, year int, month int, day, hours, minutes, seconds ``` Calculates position and velocity vectors for given time -#### func ThetaG_JD +#### func ThetaGFromJulianDate ```go -func ThetaG_JD(jday float64) (ret float64) +func ThetaGFromJulianDate(jday float64) (ret float64) ``` Calculate GMST from Julian date. Reference: The 1992 Astronomical Almanac, page B6. diff --git a/conversions.go b/conversions.go index eefb2c8..752ccc3 100644 --- a/conversions.go +++ b/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,7 +37,7 @@ func days2mdhms(year int64, epochDays float64) (mon, day, hr, min, sec float64) return } -// Calc julian date given year, month, day, hour, minute and second +// JDay 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. 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) @@ -47,22 +47,22 @@ func JDay(year, mon, day, hr, min, sec int) float64 { func gstime(jdut1 float64) (temp float64) { tut1 := (jdut1 - 2451545.0) / 36525.0 temp = -6.2e-6*tut1*tut1*tut1 + 0.093104*tut1*tut1 + (876600.0*3600+8640184.812866)*tut1 + 67310.54841 - temp = math.Mod((temp * DEG2RAD / 240.0), TWOPI) + temp = math.Mod((temp * deg2rad / 240.0), twoPi) if temp < 0.0 { - temp += TWOPI + temp += twoPi } return } -// Calc GST given year, month, day, hour, minute and second +// GSTimeFromDate calc 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. +// ThetaGFromJulianDate calculates GMST from Julian date. // Reference: The 1992 Astronomical Almanac, page B6. -func ThetaG_JD(jday float64) (ret float64) { +func ThetaGFromJulianDate(jday float64) (ret float64) { _, UT := math.Modf(jday + 0.5) 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(ThetaGFromJulianDate(jday)+obsCoords.Longitude, twoPi) r := (re + alt) * math.Cos(obsCoords.Latitude) eciObs.X = r * math.Cos(theta) eciObs.Y = r * math.Sin(theta) @@ -135,7 +135,7 @@ 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) @@ -144,30 +144,30 @@ 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) + theta := math.Mod(ThetaGFromJulianDate(jday)+obsCoords.Longitude, 2*math.Pi) obsPos := LLAToECI(obsCoords, obsAlt, jday) rx := eciSat.X - obsPos.X ry := eciSat.Y - obsPos.Y rz := eciSat.Z - obsPos.Z - top_s := math.Sin(obsCoords.Latitude)*math.Cos(theta)*rx + math.Sin(obsCoords.Latitude)*math.Sin(theta)*ry - math.Cos(obsCoords.Latitude)*rz - 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 + topS := math.Sin(obsCoords.Latitude)*math.Cos(theta)*rx + math.Sin(obsCoords.Latitude)*math.Sin(theta)*ry - math.Cos(obsCoords.Latitude)*rz + topE := -math.Sin(theta)*rx + math.Cos(theta)*ry + topZ := 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 = math.Atan(-topE / topS) + if topS > 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) + lookAngles.El = math.Asin(topZ / lookAngles.Rg) return } diff --git a/dspace.go b/dspace.go index 8a73fc3..8a80d79 100644 --- a/dspace.go +++ b/dspace.go @@ -4,12 +4,12 @@ import ( "math" ) -// A struct returned from the dsinit function +// DeepSpaceInitResult a 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 a struct returned from the dspace function type DeepSpaceResult struct { atime, em, argpm, inclm, xli, mm, xni, nodem, dndt, nm float64 } @@ -34,7 +34,6 @@ func dsinit(whichconst GravConst, cosim, emsq, argpo, s1, s2, s3, s4, s5, sinim, xke := whichconst.xke - irez = 0 if 0.0034906585 < nm && nm < 0.0052359877 { irez = 1 } @@ -74,7 +73,7 @@ func dsinit(whichconst GravConst, cosim, emsq, argpo, s1, s2, s3, s4, s5, sinim, } dndt := 0.0 - theta = math.Mod(gsto+tc*rptim, TWOPI) + theta = math.Mod(gsto+tc*rptim, twoPi) em = em + dedt*t inclm = inclm + didt*t argpm = argpm + domdt*t @@ -154,7 +153,7 @@ func dsinit(whichconst GravConst, cosim, emsq, argpo, s1, s2, s3, s4, s5, sinim, temp = 2.0 * temp1 * root54 d5421 = temp * f542 * g521 d5433 = temp * f543 * g533 - xlamo = math.Mod(mo+nodeo+nodeo-theta-theta, TWOPI) + xlamo = math.Mod(mo+nodeo+nodeo-theta-theta, twoPi) xfact = mdot + dmdt + 2.0*(nodedot+dnodt-rptim) - no em = emo emsq = emsqo @@ -171,7 +170,7 @@ func dsinit(whichconst GravConst, cosim, emsq, argpo, s1, s2, s3, s4, s5, sinim, del2 = 2.0 * del1 * f220 * g200 * q22 del3 = 3.0 * del1 * f330 * g300 * q33 * aonv del1 = del1 * f311 * g310 * q31 * aonv - xlamo = math.Mod(mo+nodeo+argpo-theta, TWOPI) + xlamo = math.Mod(mo+nodeo+argpo-theta, twoPi) xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - no } xli = xlamo @@ -233,7 +232,7 @@ func dspace(irez, d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, step2 := 259200.0 dndt := 0.0 - theta = math.Mod((gsto + tc*rptim), TWOPI) + theta = math.Mod((gsto + tc*rptim), twoPi) em = em + dedt*t inclm = inclm + didt*t @@ -314,7 +313,7 @@ func dspace(irez, d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, return } -// A struct returned from the dpper function +// DpperResult a struct returned from the dpper function type DpperResult struct { ep, inclp, nodep, argpp, mp float64 } @@ -423,21 +422,21 @@ func dpper(satrec *Satellite, inclo float64, init string, ep, inclp, nodep, argp alfdp = alfdp + dalf betdp = betdp + dbet - nodep = math.Mod(nodep, TWOPI) + nodep = math.Mod(nodep, twoPi) if nodep < 0.0 && opsmode == "a" { - nodep = nodep + TWOPI + nodep = nodep + twoPi } xls := mp + argpp + pl + pgh + (cosip-pinc*sinip)*nodep xnoh := nodep nodep = math.Atan2(alfdp, betdp) if nodep < 0.0 && opsmode == "a" { - nodep = nodep + TWOPI + nodep = nodep + twoPi } if math.Abs(xnoh-nodep) > math.Pi { if nodep < xnoh { - nodep = nodep + TWOPI + nodep = nodep + twoPi } else { - nodep = nodep - TWOPI + nodep = nodep - twoPi } } mp += pl @@ -454,7 +453,7 @@ func dpper(satrec *Satellite, inclo float64, init string, ep, inclp, nodep, argp return } -// A struct returned from the dscom function +// DSComResults a 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 } @@ -490,7 +489,7 @@ func dscom(epoch, ep, argpp, tc, inclp, nodep, np, e3, ee2, peo, pgho, pho, pinc pgho = 0.0 pho = 0.0 day := epoch + 18261.5 + tc/1440.0 - xnodce = math.Mod(4.5236020-9.2422029e-4*day, TWOPI) + xnodce = math.Mod(4.5236020-9.2422029e-4*day, twoPi) stem = math.Sin(xnodce) ctem = math.Cos(xnodce) zcosil = 0.91375164 - 0.03568096*ctem @@ -588,8 +587,8 @@ func dscom(epoch, ep, argpp, tc, inclp, nodep, np, e3, ee2, peo, pgho, pho, pinc } } - zmol = math.Mod(4.7199672+0.22997150*day-gam, TWOPI) - zmos = math.Mod(6.2565837+0.017201977*day, TWOPI) + zmol = math.Mod(4.7199672+0.22997150*day-gam, twoPi) + zmos = math.Mod(6.2565837+0.017201977*day, twoPi) se2 = 2.0 * ss1 * ss6 se3 = 2.0 * ss1 * ss7 diff --git a/helpers.go b/helpers.go index 504656f..bc8db4f 100644 --- a/helpers.go +++ b/helpers.go @@ -7,28 +7,28 @@ import ( "strings" ) -// 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 ( + twoPi float64 = math.Pi * 2.0 + deg2rad float64 = math.Pi / 180.0 + 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,23 +58,23 @@ 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) opsmode := "i" - sat.no = sat.no / XPDOTP - sat.ndot = sat.ndot / (XPDOTP * 1440.0) - sat.nddot = sat.nddot / (XPDOTP * 1440.0 * 1440) + sat.no = sat.no / xPDOTP + sat.ndot = sat.ndot / (xPDOTP * 1440.0) + sat.nddot = sat.nddot / (xPDOTP * 1440.0 * 1440) - sat.inclo = sat.inclo * DEG2RAD - sat.nodeo = sat.nodeo * DEG2RAD - sat.argpo = sat.argpo * DEG2RAD - sat.mo = sat.mo * DEG2RAD + sat.inclo = sat.inclo * deg2rad + sat.nodeo = sat.nodeo * deg2rad + 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..e96b545 100644 --- a/satellite.go +++ b/satellite.go @@ -24,8 +24,6 @@ type Satellite struct { argpo float64 mo float64 no float64 - alta float64 - altp float64 method string operationmode string diff --git a/sgp4.go b/sgp4.go index d7cbca8..ffb6513 100644 --- a/sgp4.go +++ b/sgp4.go @@ -194,12 +194,6 @@ func sgp4init(opsmode *string, epoch float64, satrec *Satellite) (position, velo dsinitResults := dsinit(satrec.whichconst, cosim, emsq, satrec.argpo, s1, s2, s3, s4, s5, sinim, ss1, ss2, ss3, ss4, ss5, sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33, satrec.t, tc, satrec.gsto, satrec.mo, satrec.mdot, satrec.no, satrec.nodeo, satrec.nodedot, xpidot, z1, z3, z11, z13, z21, z23, z31, z33, satrec.ecco, eccsq, em, argpm, inclm, mm, nm, nodem, satrec.irez, satrec.atime, satrec.d2201, satrec.d2211, satrec.d3210, satrec.d3222, satrec.d4410, satrec.d4422, satrec.d5220, satrec.d5232, satrec.d5421, satrec.d5433, satrec.dedt, satrec.didt, satrec.dmdt, satrec.dnodt, satrec.domdt, satrec.del1, satrec.del2, satrec.del3, satrec.xfact, satrec.xlamo, satrec.xli, satrec.xni) - em = dsinitResults.em - argpm = dsinitResults.argpm - inclm = dsinitResults.inclm - mm = dsinitResults.mm - nm = dsinitResults.nm - nodem = dsinitResults.nodem satrec.irez = dsinitResults.irez satrec.atime = dsinitResults.atime satrec.d2201 = dsinitResults.d2201 @@ -279,10 +273,10 @@ func initl(satn int64, grav GravConst, ecco, epoch, inclo, noIn float64, methodI c1 := 1.72027916940703639e-2 thgr70 := 1.7321343856509374 fk5r := 5.07551419432269442e-15 - c1p2p := c1 + TWOPI - gsto = math.Mod((thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r), TWOPI) + c1p2p := c1 + twoPi + gsto = math.Mod((thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r), twoPi) if gsto < 0.0 { - gsto = gsto + TWOPI + gsto = gsto + twoPi } } else { gsto = gstime(epoch + 2433281.5) @@ -302,7 +296,7 @@ func Propagate(sat Satellite, year int, month int, day, hours, minutes, seconds // satrec - initialized Satellite struct from sgp4init // tsince - time since epoch in minutes func sgp4(satrec *Satellite, tsince float64) (position, velocity Vector3) { - var am, axnl, aynl, betal, cosim, sinim, cnod, snod, cos2u, sin2u, coseo1, sineo1, cosi, sini, cosip, sinip, cosisq, cossu, sinsu, cosu, sinu, delm, delomg, emsq, ecose, el2, eo1, esine, argpm, argpp, pl, rdotl, rl, rvdot, rvdotl, su, t2, t3, t4, tc, tem5, temp, temp1, temp2, tempa, tempe, templ, u, ux, uy, uz, vx, vy, vz, inclm, mm, nm, nodem, xinc, xincp, xl, xlm, mp, xmdf, xmx, xmy, nodedf, xnode, nodep, mrt float64 + var am, axnl, aynl, betal, cosim, sinim, cnod, snod, cos2u, sin2u, coseo1, sineo1, cosi, sini, cosip, sinip, cosisq, cossu, sinsu, cosu, sinu, delm, delomg, ecose, el2, eo1, esine, argpm, argpp, pl, rdotl, rl, rvdot, rvdotl, su, t2, t3, t4, tc, tem5, temp, temp1, temp2, tempa, tempe, templ, u, ux, uy, uz, vx, vy, vz, inclm, mm, nm, nodem, xinc, xincp, xl, xlm, mp, xmdf, xmx, xmy, nodedf, xnode, nodep, mrt float64 mrt = 0.0 temp4 := 1.5e-12 @@ -380,13 +374,11 @@ func sgp4(satrec *Satellite, tsince float64) (position, velocity Vector3) { } mm = mm + satrec.no*templ xlm = mm + argpm + nodem - emsq = em * em - temp = 1.0 - emsq - nodem = math.Mod(nodem, TWOPI) - argpm = math.Mod(argpm, TWOPI) - xlm = math.Mod(xlm, TWOPI) - mm = math.Mod((xlm - argpm - nodem), TWOPI) + nodem = math.Mod(nodem, twoPi) + argpm = math.Mod(argpm, twoPi) + xlm = math.Mod(xlm, twoPi) + mm = math.Mod((xlm - argpm - nodem), twoPi) sinim = math.Sin(inclm) cosim = math.Cos(inclm) @@ -436,7 +428,7 @@ func sgp4(satrec *Satellite, tsince float64) (position, velocity Vector3) { aynl = ep*math.Sin(argpp) + temp*satrec.aycof xl = mp + argpp + nodep + temp*satrec.xlcof*axnl - u = math.Mod((xl - nodep), TWOPI) + u = math.Mod((xl - nodep), twoPi) eo1 = u tem5 = 9999.9 ktr := 1