Notice. New forum software under development. It's going to miss a few functions and look a bit ugly for a while, but I'm working on it full time now as the old forum was too unstable. Couple days, all good. If you notice any issues, please contact me.
|
Forum Index : Microcontroller and PC projects : MMX/RPI/DOS - repeating groundtrack orbit
Author | Message | ||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
This post describes a MMBASIC computer program called repeat3.bas that can be used for the preliminary design of repeating groundtrack orbits. These types of orbits are useful for remote sensing and other special applications since they overfly the same points on the Earth’s surface every repeat cycle. Repeating ground track orbits are usually specified by an integer number of days N and integer number of orbits K in the repeat cycle. This program calculates the mean semimajor axis required for a repeating ground track orbit using an algorithm devised by Carl Wagner. This iterative numerical method is described in “A Prograde Geosat Exact Repeat Mission?”, The Journal of the Astronautical Sciences, Vol. 39, No. 3, July-September 1991, pp. 313-326. Here's a typical user interaction with the software and the displayed results. program repeat3 repeating ground track - Wagner's method ---------------------------------------- please input the mean orbital eccentricity (non-dimensional) (0 <= eccentricity < 1) ? 0 please input the mean orbital inclination (degrees) (0 <= inclination < 180) ? 108 please input the number of orbits in the repeat cycle ? 271 please input the number of nodal days in the repeat cycle ? 19 program repeat3.bas repeating ground track - Wagner's method ---------------------------------------- mean semimajor axis 7192.2303 kilometers mean eccentricity 0.0000000000 mean inclination 108.0000 degrees orbits to repeat 271 solar days to repeat 19.0548 Keplerian period 101.1708 minutes nodal period 101.2507 minutes length of nodal day 1444.1545 minutes fundamental interval 25.2399 degrees Here's the MMBASIC source code ' program repeat3.bas August 3, 2017 ' determine mean semimajor axis required ' for a repeating ground track orbit ' Wagner's algorithm ' MMX/RasPi/DOS versions of MMBASIC ''''''''''''''''''''''''''''''''''' option default float option base 1 ' astrodynamic and utility constants const pi2 = 2.0 * pi ' conversion factor - degrees to radians const dtr = pi / 180.0 ' conversion factor - radians to degrees const rtd = 180.0 / pi ' Earth gravitational constant (km**3/sec**2) const xmu = 398600.4415 ' equatorial radius of the Earth (kilometers) const req = 6378.14 ' first zonal coefficient (non-dimensional) const xj2 = 1.08262668355e-3 ' Earth rotation rate (radians/second) const omega = 7.29211564186e-05 print " " print "program repeat3" print " " print "repeating ground track - Wagner's method" print "----------------------------------------" ' request mean orbital elements do print " " print "please input the mean orbital eccentricity (non-dimensional)" print "(0 <= eccentricity < 1)" input ecc if (ecc >= 0.0 and ecc < 1.0) then exit do loop do print " " print "please input the mean orbital inclination (degrees)" print "(0 <= inclination < 180)" input xdata if (xdata >= 0.0 and xdata < 180.0) then exit do loop xinc = dtr * xdata do print " " print "please input the number of orbits in the repeat cycle" input xnorbits if (xnorbits > 0) then exit do loop do print " " print "please input the number of nodal days in the repeat cycle" input xndays if (xndays > 0) then exit do loop c20 = -xj2 ' calculate initial guess for semimajor axis (kilometers) a0 = xmu^(1.0 / 3.0) * ((xnorbits / xndays) * omega)^(-2.0 / 3.0) aold = a0 * (1.0 - c20 * (req / a0)^2 * (4.0 * cos(xinc)^2 - (xnorbits / xndays) * cos(xinc) - 1.0)) do slr = aold * (1.0 - ecc * ecc) tmp1 = xmu^(1.0 / 3.0) * (xnorbits * omega / xndays)^(-2.0 / 3.0) tmp2 = (1.0 - 1.5 * c20 * (req / slr)^2 * (1.0 - 1.5 * sin(xinc)^2))^(2.0 / 3.0) tmp3 = (1.0 + c20 * (req / slr)^2 * (1.5 * (xnorbits / xndays) * cos(xinc) - 0.75 * (5.0 * cos(xinc)^2 - 1.0)))^(2.0 / 3.0) anew = tmp1 * tmp2 * tmp3 if (abs(anew - aold) < 1.0e-6) then ' converged to within tolerance exit do else ' continue iteration aold = anew end if loop sma = anew ' Keplerian period (seconds) tkepler = pi2 * sma * sqr(sma / xmu) ' keplerian mean motion (rad/sec) xmm = sqr(xmu / (sma * sma * sma)) ' orbital semiparameter (km) slr = sma * (1.0 - ecc * ecc) b = sqr(1.0 - ecc * ecc) c = req / slr d = c * c e = sin(xinc) * sin(xinc) ' perturbed mean motion (rad/sec) pmm = xmm * (1.0 + 1.5 * xj2 * d * b * (1.0 - 1.5 * e)) ' argument of perigee perturbation (rad/sec) apdot = 1.5 * xj2 * pmm * d * (2.0 - 2.5 * e) ' raan perturbation (rad/sec) raandot = -1.5 * xj2 * pmm * d * cos(xinc) ' nodal period - time between nodal crossings tnode = pi2 / (apdot + pmm) ' delta longitude per nodal period (radians) dlong = tnode * (omega - raandot) ' length of nodal day (seconds) tnday = pi2 / (omega - raandot) ' number of solar days to repeat xnsdays = tnode * xnorbits / 86400.0 ' print results print " " print "program repeat3.bas" print " " print "repeating ground track - Wagner's method" print "----------------------------------------" print " " print "mean semimajor axis ", str$(sma, 0, 4) + " kilometers" print " " print "mean eccentricity ", str$(ecc, 0, 10) print " " print "mean inclination ", str$(rtd * xinc, 0, 4) + " degrees" print " " print " " print "orbits to repeat ", str$(xnorbits) print " " print "solar days to repeat ", str$(xnsdays) print " " print " " print "Keplerian period ", str$(tkepler / 60.0, 0, 4) + " minutes" print " " print "nodal period ", str$(tnode / 60.0, 0, 4) + " minutes" print " " print " " print "length of nodal day ", str$(tnday / 60.0, 0, 4) + " minutes" print " " print "fundamental interval ", str$(dlong * rtd, 0, 4) + " degrees" print " " end Here's a PDF document that describes the computer programs that comprise the repeating groundtrack series. 2017-08-06_104742_repeat.zip |
||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
' repeat1.bas August 7, 2017 ' time to repeat ground track ' Kozai analytic orbit propagation ' MMBASIC for MMX, RasPi and DOS Program example. program repeat1 - time to repeat ground track - analytic solution please input the semimajor axis (kilometers) (semimajor axis > 0) ? 8000 please input the orbital eccentricity (non-dimensional) (0 <= eccentricity < 1) ? 0.025 please input the orbital inclination (degrees) (0 <= inclination <= 180) ? 28.5 please input the argument of perigee (degrees) (0 <= argument of perigee <= 360) ? 120 please input the closure tolerance (degrees) (a value between 0.1 and 0.5 degrees is recommended) ? 0.1 program repeat1 - time to repeat ground track - analytic solution mean semimajor axis 8000.000000 kilometers mean eccentricity (nd) 0.02500000 mean inclination 28.500000 degrees mean argument of perigee 120.000000 degrees number of orbits to repeat 2027 number of solar days to repeat 166.705086 Keplerian period 118.684693 minutes nodal period 118.428872 minutes length of nodal day 1420.446851 minutes fundamental interval 30.014776 degrees closure tolerance 0.100000 degrees actual closure 0.049141 degrees MMBASIC source code listing. ' repeat1.bas August 7, 2017 ' time to repeat ground track ' Kozai analytic orbit propagation ' MMBASIC for MMX, RasPi and DOS '''''''''''''''''''''''''''''''' option default float option base 1 dim ioev%(6), oev1(7), oev2(7), rs(3), vs(3) ' global constants and variables dim mm, pmm, apdot, raandot, sinc, cinc, slr ' zonal gravity constant (nd) const xj2 = 0.00108263 ' Earth equatorial radius (kilometers) const req = 6378.14 ' Earth gravitational constant (km^3/sec^2) const xmu = 398600.4415 ' Earth inertial rotation rate (rad/sec) const omega = 7.2921151467e-5 const pi2 = 2.0 * pi const dtr = pi / 180.0 const rtd = 180.0 / pi print " " print "program repeat1 - time to repeat ground track - analytic solution" ''''''''''''''''''''''''''''''''''''''''' ' request mean classical orbital elements ''''''''''''''''''''''''''''''''''''''''' ioev%(1) = 1 ioev%(2) = 1 ioev%(3) = 1 ioev%(4) = 1 ioev%(5) = 0 ioev%(6) = 0 getoe(ioev%(), oev1()) do print " " print "please input the closure tolerance (degrees)" print "(a value between 0.1 and 0.5 degrees is recommended)" input tol if (tol > 0.0) then exit do loop ' set true anomaly to ascending node oev1(6) = -oev1(4) ' calculate initial mean anomaly (radians) a = sqr(1.0 - oev1(2) * oev1(2)) * sin(oev1(6)) b = cos(oev1(6)) + oev1(2) eanom = atan3(a, b) oev1(6) = modulo(eanom - oev1(2) * sin(eanom)) ' compute Kozai mean motion and perturbations t = 0.0 kozai1(1, t, oev1(), oev2(), rs(), vs()) ' Keplerian period (seconds) tkepler = pi2 / mm ' nodal period => time between nodal crossings tnode = pi2 / (apdot + pmm) ' delta longitude per nodal period (radians) dlong = tnode * (omega - raandot) ' length of nodal day (seconds) tnday = pi2 / (omega - raandot) ' initialize xlong = 0.0 norbits = 0 ' increment nodal crossing and look for closure do ' increment current longitude xlong = xlong + dlong ' modulo longitude if (xlong >= pi2) then xlong = xlong - pi2 ' increment number of orbits norbits = norbits + 1 ' check for ground track closure if ((abs(xlong - pi2) <= dtr * tol) or (abs(xlong) <= dtr * tol)) then exit do loop ' number of days to repeat ground track ndays = tnode * norbits / 86400.0 ' actual closure delta longitude (radians) delta = abs(xlong - pi2) if (abs(xlong) < delta) then delta = abs(xlong) ' print results print " " print "program repeat1 - time to repeat ground track - analytic solution" print " " print "mean semimajor axis ", str$(oev1(1), 0, 6) + " kilometers" print " " print "mean eccentricity (nd) ", str$(oev1(2), 0, 8) print " " print "mean inclination ", str$(rtd * oev1(3), 0, 6) + " degrees" print " " print "mean argument of perigee ", str$(rtd * oev1(4), 0, 6) + " degrees" print " " print "number of orbits to repeat ", norbits print " " print "number of solar days to repeat ", str$(ndays, 0, 6) print " " print "Keplerian period ", str$(tkepler / 60.0, 0, 6) + " minutes" print " " print "nodal period ", str$(tnode / 60.0, 0, 6) + " minutes" print " " print "length of nodal day ", str$(tnday / 60.0, 0, 6) + " minutes" print " " print "fundamental interval ", str$(dlong * rtd, 0, 6) + " degrees" print " " print "closure tolerance ", str$(tol, 0, 6) + " degrees" print " " print "actual closure ", str$(delta * rtd, 0, 6) + " degrees" print " " end '''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''' sub kozai1(iniz%, t, oev1(), oev2(), r(), v()) ' analytic orbit propagation ' Kozai's method - ECI version ' input ' iniz = initialization flag (1 = initialize, 0 = bypass) ' t = elapsed simulation time (seconds) ' initial orbital elements ' oev1(1) = semimajor axis (kilometers) ' oev1(2) = orbital eccentricity (non-dimensional) ' (0 <= eccentricity < 1) ' oev1(3) = orbital inclination (radians) ' (0 <= oev1(3) <= pi) ' oev1(4) = argument of perigee (radians) ' (0 <= oev1(4) <= 2 pi) ' oev1(5) = right ascension of ascending node (radians) ' (0 <= oev1(5) <= 2 pi) ' oev1(6) = mean anomaly (radians) ' (0 <= oev1(6) <= 2 pi) ' output ' final orbital elements and state vector at time = t ' oev2(1) = semimajor axis (kilometers) ' oev2(2) = orbital eccentricity (non-dimensional) ' (0 <= eccentricity < 1) ' oev2(3) = orbital inclination (radians) ' (0 <= oev2(3) <= pi) ' oev2(4) = updated argument of perigee (radians) ' (0 <= oev2(4) <= 2 pi) ' oev2(5) = updated right ascension of ascending node (radians) ' (0 <= oev2(5) <= 2 pi) ' oev2(6) = updated mean anomaly (radians) ' (0 <= oev2(6) <= 2 pi) ' oev2(7) = updated true anomaly (radians) ' (0 <= oev2(6) <= 2 pi) ' r = eci position vector (kilometers) ' v = eci velocity vector (kilometers/second) '''''''''''''''''''''''''''''''''''''''''''''''''''' local b, c, d, e, c1, c2, c3, rmag local argperi, raani, manomi ' set final orbital elements equal to initial for i% =1 to 6 oev2(i%) = oev1(i%) next i% if (iniz% = 1) then sinc = sin(oev1(3)) cinc = cos(oev1(3)) ' keplerian mean motion mm = sqr(xmu / (oev1(1) * oev1(1) * oev1(1))) ' orbital semiparameter slr = oev1(1) * (1.0 - oev1(2) * oev1(2)) b = sqr(1.0 - oev1(2) * oev1(2)) c = req / slr d = c * c e = sinc * sinc ' perturbed mean motion pmm = mm * (1.0 + 1.5 * xj2 * d * b * (1.0 - 1.5 * e)) ' argument of perigee perturbation apdot = 1.5 * xj2 * pmm * d * (2.0 - 2.5 * e) ' raan perturbation raandot = -1.5 * xj2 * pmm * d * cinc end if ' set initial values of argument of perigee, ' raan and mean anomaly argperi = oev1(4) raani = oev1(5) manomi = oev1(6) ' propagate orbit argper = modulo(argperi + apdot * t) raan = modulo(raani + raandot * t) manom = modulo(manomi + pmm * t) oev2(4) = argper oev2(5) = raan oev2(6) = manom ' solve Kepler's equation using Danby's method ecc = oev1(2) kepler1(manom, ecc, ea, tanom) oev2(7) = tanom ' position magnitude rmag = slr / (1.0 + ecc * cos(tanom)) sargper = sin(argper) cargper = cos(argper) ' argument of latitude arglat = modulo(argper + tanom) sarglat = sin(arglat) carglat = cos(arglat) sraan = sin(raan) craan = cos(raan) ' eci position vector r(1) = rmag * (craan * carglat - cinc * sarglat * sraan) r(2) = rmag * (sraan * carglat + cinc * sarglat * craan) r(3) = rmag * sinc * sarglat ' eci velocity vector c1 = sqr(xmu / slr) c2 = ecc * cargper + carglat c3 = ecc * sargper + sarglat v(1) = -c1 * (craan * c3 + sraan * cinc * c2) v(2) = -c1 * (sraan * c3 - craan * cinc * c2) v(3) = c1 * c2 * sinc end sub '''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''' sub kepler1 (manom, ecc, eanom, tanom, niter%) ' solve Kepler's equation for circular, ' elliptic and hyperbolic orbits ' Danby's method ' input ' manom = mean anomaly (radians) ' ecc = orbital eccentricity (non-dimensional) ' output ' eanom = eccentric anomaly (radians) ' tanom = true anomaly (radians) ' niter% = number of algorithm iterations '''''''''''''''''''''''''''''''''''''''''' local ktol, xma, s, c, f, fp, fpp, fppp local delta, deltastar, deltak, sta, cta ' define convergence criterion ktol = 1.0e-10 xma = manom - pi2 * fix(manom / pi2) ' initial guess if (ecc = 0.0) then ' circular orbit tanom = xma eanom = xma return elseif (ecc < 1.0) then ' elliptic orbit eanom = xma + 0.85 * sgn(sin(xma)) * ecc else ' hyperbolic orbit eanom = log(2.0 * xma / ecc + 1.8) end if ' perform iterations niter% = 0 do if (ecc < 1.0) then ' elliptic orbit s = ecc * sin(eanom) c = ecc * cos(eanom) f = eanom - s - xma fp = 1.0 - c fpp = s fppp = c else ' hyperbolic orbit s = ecc * sinh(eanom) c = ecc * cosh(eanom) f = s - eanom - xma fp = c - 1.0 fpp = s fppp = c end if niter% = niter% + 1 ' check for convergence if (abs(f) <= ktol or niter% > 20) then exit do ' update eccentric anomaly delta = -f / fp deltastar = -f / (fp + 0.5 * delta * fpp) deltak = -f / (fp + 0.5 * deltastar * fpp + deltastar * deltastar * fppp / 6.0) eanom = eanom + deltak loop if (niter% > 20) then print " " print "more than 20 iterations in kepler1" stop end if ' compute true anomaly if (ecc < 1.0) then ' elliptic orbit sta = sqr(1.0 - ecc * ecc) * sin(eanom) cta = cos(eanom) - ecc else ' hyperbolic orbit sta = sqr(ecc * ecc - 1.0) * sinh(eanom) cta = ecc - cosh(eanom) end if tanom = atan3(sta, cta) end sub ''''''''''''''''''''''' ''''''''''''''''''''''' sub getoe(ioev%(), oev()) ' interactive request of classical orbital elements ' input ' ioev%(1) = request semimajor axis (1 = yes, 0 = no) ' ioev%(2) = request orbital eccentricity (1 = yes, 0 = no) ' ioev%(3) = request orbital inclination (1 = yes, 0 = no) ' ioev%(4) = request argument of perigee (1 = yes, 0 = no) ' ioev%(5) = request right ascension of the ascending node (1 = yes, 0 = no) ' ioev%(6) = request true anomaly (1 = yes, 0 = no) ' output ' oev(1) = semimajor axis (kilometers) ' oev(2) = orbital eccentricity ' oev(3) = orbital inclination (radians) ' oev(4) = argument of perigee (radians) ' oev(5) = right ascension of the ascending node (radians) ' oev(6) = true anomaly (radians) '''''''''''''''''''''''''''''''''' local sma, ecc, orbinc, argper, raan, tanom ' initialize all orbital elements for i% = 1 to 6 oev(i%) = 0.0 next i% if (ioev%(1) = 1) then ' semimajor axis (kilometers) ' (semimajor axis > req) do print " " print "please input the semimajor axis (kilometers)" print "(semimajor axis > 0)" input sma if (sma > 0.0) then exit do loop oev(1) = sma end if if (ioev%(2) = 1) then ' orbital eccentricity (non-dimensional) ' (0 <= eccentricity < 1) do print " " print "please input the orbital eccentricity (non-dimensional)" print "(0 <= eccentricity < 1)" input ecc if (ecc >= 0.0 and ecc < 1.0) then exit do loop oev(2) = ecc end if if (ioev%(3) = 1) then ' orbital inclination (degrees) ' (0 <= inclination <= 180) do print " " print "please input the orbital inclination (degrees)" print "(0 <= inclination <= 180)" input orbinc if (orbinc >= 0.0 and orbinc <= 180.0) then exit do loop oev(3) = dtr * orbinc end if if (ioev%(4) = 1) then ' argument of perigee (degrees) ' (0 <= argument of perigee <= 360) if (ecc > 0.0) then do print " " print "please input the argument of perigee (degrees)" print "(0 <= argument of perigee <= 360)" input argper if (argper >= 0.0 and argper <= 360.0) then exit do loop else argper = 0.0 end if oev(4) = dtr * argper end if if (ioev%(5) = 1) then ' raan (degrees) ' (0 <= raan <= 360) do print " " print "please input the right ascension of the ascending node (degrees" print "(0 <= raan <= 360)" input raan if (raan >= 0.0 and raan <= 360.0) then exit do loop oev(5) = dtr * raan end if if (ioev%(6) = 1) then ' true anomaly (degrees) ' (0 <= true anomaly <= 360) do print " " print "please input the true anomaly (degrees)" print "(0 <= true anomaly <= 360)" input tanom if (tanom >= 0.0 and tanom <= 360.0) then exit do loop oev(6) = dtr * tanom end if end sub ''''''''''''''''''''''''''' ''''''''''''''''''''''''''' function atan3(a, b) as float ' four quadrant inverse tangent function ' input ' a = sine of angle ' b = cosine of angle ' output ' atan3 = angle (0 =< atan3 <= 2 * pi; radians) '''''''''''''''''''''''''''''''''''''''''''''''' local c if (abs(a) < 1.0e-10) then atan3 = (1.0 - sgn(b)) * pidiv2 exit function else c = (2.0 - sgn(a)) * pidiv2 endif if (abs(b) < 1.0e-10) then atan3 = c exit function else atan3 = c + sgn(a) * sgn(b) * (abs(atn(a / b)) - pidiv2) endif end function ''''''''''''''''''''''''' ''''''''''''''''''''''''' function modulo(x) as float ' modulo 2 pi function '''''''''''''''''''''' local a a = x - pi2 * fix(x / pi2) if (a < 0.0) then a = a + pi2 end if modulo = a end function |
||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
Updated links to the software on my Google Drive. Thanks to chronic for noticing that they were broken. MMBASIC software Octave software |
||||
WhiteWizzard Guru Joined: 05/04/2013 Location: United KingdomPosts: 2794 |
Hi David, Off Topic (sorry!) I meant to contact you last week to check everything is ok with you out there! I really hope you weren't too affected with things, and that you were indeed able to stay safe. Great to see more mind-blowing stuff from you . . . WW For everything Micromite visit micromite.org Direct Email: whitewizzard@micromite.o |
||||
Print this page |