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 - interplanetary injection
Author | Message | ||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
This post describes a MMBASIC program that solves one of my favorite problems in geometry and orbital mechanics. The program is named hyper1.bas and tackles the problem of interplanetary injection from a circular Earth orbit. The numerical method is due to Professor Richard Battin. The solution of this problem provides an initial estimate of the delta-v and energy required to transfer a spacecraft onto an interplanetary trajectory (perhaps to Mars) from a circular orbit about the Earth. It also provides the approximate geometry of the possible Earth departure hyperbolas. Here's a typical user interaction and program output ... Interplanetary Injection from a Circular Earth Park Orbit --------------------------------------------------------- please input the altitude of the circular Earth park orbit (kilometers) ? 185.32 please input the orbital inclination of the Earth park orbit (degrees) (0 <= inclination <= 180) ? 28.5 please input the C3 of the departure hyperbola (kilometers^2/second^2) (C3 > 0) ? 9.28 please input the right ascension of the outgoing asymptote (degrees) (0 <= right ascension <= 360) ? 352.59 please input the declination of the outgoing asymptote (degrees) ? 2.27 --------------------------------------------------------- Interplanetary Injection from a Circular Earth Park Orbit --------------------------------------------------------- departure hyperbola characteristics ----------------------------------- c3 9.2800 kilometers^2/second^2 asymptote right ascension 352.590000 degrees asymptote declination 2.270000 degrees orbital elements and state vector of park orbit at injection - opportunity #1 ----------------------------------------------------------------------------- semimajor axis 6563.4600 kilometers eccentricity 0.0000000000 orbital inclination 28.5000 degrees argument of perigee 0.0000 degrees RAAN 176.7767 degrees true anomaly 25.0750 degrees position vector x-component -6.072921967203e+03 kilometers position vector y-component -2.106409915567e+03 kilometers position vector z-component 1.327276617537e+03 kilometers position vector magnitude 6.563460000000e+03 kilometers velocity vector x-component 2.948684715564e+00 kilometers/second velocity vector y-component -6.379019476975e+00 kilometers/second velocity vector z-component 3.368026111924e+00 kilometers/second velocity vector magnitude 7.792960344441e+00 kilometers/second orbital elements and state vector of hyperbola at injection - opportunity #1 ---------------------------------------------------------------------------- semimajor axis -42952.6338 kilometers eccentricity 1.1528069276 orbital inclination 28.5000 degrees argument of perigee 25.0750 degrees RAAN -3.2233 degrees true anomaly 0.0000 degrees position vector x-component -6.072921967203e+03 kilometers position vector y-component -2.106409915567e+03 kilometers position vector z-component 1.327276617537e+03 kilometers position vector magnitude 6.563460000000e+03 kilometers velocity vector x-component 4.326441938393e+00 kilometers/second velocity vector y-component -9.359582340336e+00 kilometers/second velocity vector z-component 4.941718367961e+00 kilometers/second velocity vector magnitude 1.143417954468e+01 kilometers/second injection delta-v vector and magnitude - opportunity #1 ------------------------------------------------------- x-component of delta-v 1377.757223 meters/second y-component of delta-v -2980.562863 meters/second z-component of delta-v 1573.692256 meters/second delta-v magnitude 3641.219200 meters/second orbital elements and state vector of park orbit at injection - opportunity #2 ----------------------------------------------------------------------------- semimajor axis 6563.4600 kilometers eccentricity 0.0000000000 orbital inclination 28.5000 degrees argument of perigee 0.0000 degrees RAAN 348.4033 degrees true anomaly 214.5981 degrees position vector x-component -5.950846004649e+03 kilometers position vector y-component -2.122286481036e+03 kilometers position vector z-component -1.778296683053e+03 kilometers position vector magnitude 6.563460000000e+03 kilometers velocity vector x-component 3.201396629997e+00 kilometers/second velocity vector y-component -6.411885875654e+00 kilometers/second velocity vector z-component -3.060883869907e+00 kilometers/second velocity vector magnitude 7.792960344441e+00 kilometers/second orbital elements and state vector of hyperbola at injection - opportunity #2 ---------------------------------------------------------------------------- semimajor axis -42952.6338 kilometers eccentricity 1.1528069276 orbital inclination 28.5000 degrees argument of perigee 34.5981 degrees RAAN -11.5967 degrees true anomaly 360.0000 degrees position vector x-component -5.950846004649e+03 kilometers position vector y-component -2.122286481036e+03 kilometers position vector z-component -1.778296683053e+03 kilometers position vector magnitude 6.563460000000e+03 kilometers velocity vector x-component 4.697232148402e+00 kilometers/second velocity vector y-component -9.407805388687e+00 kilometers/second velocity vector z-component -4.491065549808e+00 kilometers/second velocity vector magnitude 1.143417954468e+01 kilometers/second injection delta-v vector and magnitude - opportunity #2 ------------------------------------------------------- x-component of delta-v 1495.835518 meters/second y-component of delta-v -2995.919513 meters/second z-component of delta-v -1430.181680 meters/second delta-v magnitude 3641.219200 meters/second Here's a link to a PDF file of the MATLAB version with illustrations of the orbital geometry. 2017-07-20_185149_hyper1_matlab.pdf Here's the MMBASIC source code. ' hyper1.bas July 20, 2017 ' Battin's method for single impulse coplanar hyperbolic ' injection from a circular Earth park orbit ' MMBASIC for MMX, Raspberry Pi and DOS ''''''''''''''''''''''''''''''''''''''' option default float option base 1 dim pv(3), shat(3), rhat1(3), rhat2(3) dim rpo1(3), vpo1(3), oev_po1(6) dim rpo2(3), vpo2(3), oev_po2(6) dim rhyper1(3), vhyper1(3), oev_hyper1(6) dim rhyper2(3), vhyper2(3), oev_hyper2(6) dim deltav1(3), deltav2(3) ' angular conversion factors const pi2 = 2.0 * pi, rtd = 180.0 / pi, dtr = pi / 180.0 ' gravitational constant of the Earth (kilometers^3/second^2) const mu = 398600.4415 ' equatorial radius of the Earth (kilometers) const req = 6378.14 ' request user inputs print " " print "Interplanetary Injection from a Circular Earth Park Orbit" print "---------------------------------------------------------" print " " do print " " print "please input the altitude of the circular Earth park orbit (kilometers)" input alt if (alt > 0.0) then exit do loop do print " " print "please input the orbital inclination of the Earth park orbit (degrees)" print "(0 <= inclination <= 180)" input orbinc if (orbinc >= 0.0 and orbinc <= 180.0) then exit do loop do print " " print "please input the C3 of the departure hyperbola (kilometers^2/second^2)" print "(C3 > 0)" input c3 if (c3 > 0.0) then exit do loop do print " " print "please input the right ascension of the outgoing asymptote (degrees)" print "(0 <= right ascension <= 360)" input rasc if (rasc >= 0.0 and rasc <= 360.0) then exit do loop do print " " print "please input the declination of the outgoing asymptote (degrees)" input decl if (abs(decl) <= 90.0) then exit do loop if (abs(decl) >= orbinc) then print " " print "*********************************************************************" print "Please note: this program is valid for |DLA| < park orbit inclination" print "*********************************************************************" end end if ' geocentric radius of park orbit (kilometers) rpmag = req + alt ' convert angular elements to radians xinc = orbinc * dtr decl_asy = decl * dtr rasc_asy = rasc * dtr ' v-infinity (km/sec) vinf = sqr(c3) ' unit vector along the Earth's spin axis pv(1) = 0.0 pv(2) = 0.0 pv(3) = 1.0 cdecl_asy = cos(decl_asy) sdecl_asy = sin(decl_asy) crasc_asy = cos(rasc_asy) srasc_asy = sin(rasc_asy) ' compute unit asymptote vector shat(1) = cdecl_asy * crasc_asy shat(2) = cdecl_asy * srasc_asy shat(3) = sdecl_asy ' beta = angle between unit asymptote vector and spin axis beta = acos(vdot(pv(), shat())) sum1 = acos(cos(beta) / sin(xinc)) ' alpha = angle between unit asymptote vector and position vector alpha = asin(1.0 / (1.0 + rpmag * vinf * vinf / mu)) '''''''''''''''''''''' ' tangential injection '''''''''''''''''''''' ' injection argument of latitude #1 arglat1 = modulo(sum1 - alpha) ' injection argument of latitude #2 arglat2 = modulo(-sum1 - alpha) ' injection raan #1 raan1 = modulo(rasc_asy + pi + asin(cotan(beta) / tan(xinc))) ' injection raan #2 raan2 = modulo(rasc_asy + 2.0 * pi - asin(cotan(beta) / tan(xinc))) ' opportunity #1 - park orbit state vector and orbital elements oev_po1(1) = rpmag oev_po1(2) = 0.0 oev_po1(3) = xinc oev_po1(4) = 0.0 oev_po1(5) = raan1 oev_po1(6) = arglat1 orb2eci(oev_po1(), rpo1(), vpo1()) uvector(rpo1(), rhat1()) ctheta = vdot(shat(), rhat1()) d = sqr(mu / ((1.0 + ctheta) * rpmag) + vinf * vinf / 4.0) ' velocity vector of hyperbola at injection for i% = 1 to 3 vhyper1(i%) = (d + 0.5 * vinf) * shat(i%) + (d - 0.5 * vinf) * rhat1(i%) next i% ' injection delta-v vector for i% = 1 to 3 deltav1(i%) = vhyper1(i%) - vpo1(i%) rhyper1(i%) = rpo1(i%) next i% eci2orb(rhyper1(), vhyper1(), oev_hyper1()) print " " print "---------------------------------------------------------" print "Interplanetary Injection from a Circular Earth Park Orbit" print "---------------------------------------------------------" print " " print "departure hyperbola characteristics" print "-----------------------------------" print " " print "c3 ", str$(c3, 0, 4) + " kilometers^2/second^2" print " " print "asymptote right ascension ", str$(rtd * rasc_asy, 0, 6) + " degrees" print " " print "asymptote declination ", str$(rtd * decl_asy, 0, 6) + " degrees" print " " print "orbital elements and state vector of park orbit at injection - opportunity #1" print "-----------------------------------------------------------------------------" print " " oeprint(oev_po1()) svprint(rpo1(), vpo1()) print " " print "orbital elements and state vector of hyperbola at injection - opportunity #1" print "----------------------------------------------------------------------------" print " " oeprint(oev_hyper1()) svprint(rhyper1(), vhyper1()) print " " print "injection delta-v vector and magnitude - opportunity #1" print "-------------------------------------------------------" print " " print "x-component of delta-v ", str$(1000.0 * deltav1(1), 0, 6) + " meters/second" print " " print "y-component of delta-v ", str$(1000.0 * deltav1(2), 0, 6) + " meters/second" print " " print "z-component of delta-v ", str$(1000.0 * deltav1(3), 0, 6) + " meters/second" print " " print "delta-v magnitude ", str$(1000.0 * vecmag(deltav1()), 0, 6) + " meters/second" print " " ' opportunity #2 - park orbit state vector oev_po2(1) = rpmag oev_po2(2) = 0.0 oev_po2(3) = xinc oev_po2(4) = 0.0 oev_po2(5) = raan2 oev_po2(6) = arglat2 orb2eci(oev_po2(), rpo2(), vpo2()) uvector(rpo2(), rhat2()) ctheta = vdot(shat(), rhat2()) d = sqr(mu / ((1.0 + ctheta) * rpmag) + vinf * vinf / 4.0) ' velocity vector of hyperbola at injection for i% = 1 to 3 vhyper2(i%) = (d + 0.5 * vinf) * shat(i%) + (d - 0.5 * vinf) * rhat2(i%) next i% ' injection delta-v vector for i% = 1 to 3 deltav2(i%) = vhyper2(i%) - vpo2(i%) rhyper2(i%) = rpo2(i%) next i% eci2orb(rhyper2(), vhyper2(), oev_hyper2()) print " " print "orbital elements and state vector of park orbit at injection - opportunity #2" print "-----------------------------------------------------------------------------" print " " oeprint(oev_po2()) svprint(rpo2(), vpo2()) print " " print "orbital elements and state vector of hyperbola at injection - opportunity #2" print "----------------------------------------------------------------------------" print " " oeprint(oev_hyper2()) svprint(rhyper2(), vhyper2()) print " " print "injection delta-v vector and magnitude - opportunity #2" print "-------------------------------------------------------" print " " print "x-component of delta-v ", str$(1000.0 * deltav2(1), 0, 6) + " meters/second" print " " print "y-component of delta-v ", str$(1000.0 * deltav2(2), 0, 6) + " meters/second" print " " print "z-component of delta-v ", str$(1000.0 * deltav2(3), 0, 6) + " meters/second" print " " print "delta-v magnitude " str$(1000.0 * vecmag(deltav2()), 0, 6) + " meters/second" print " " end '''''''''''''''''''''''''' '''''''''''''''''''''''''' sub eci2orb(r(), v(), oev()) ' convert eci state vector to classical orbital elements ' input ' r() = eci position vector (kilometers) ' v() = eci velocity vector (kilometers/second) ' output ' oev(1) = semimajor axis (kilometers) ' oev(2) = orbital eccentricity (non-dimensional) ' (0 <= eccentricity < 1) ' oev(3) = orbital inclination (radians) ' (0 <= inclination <= pi) ' oev(4) = argument of perigee (radians) ' (0 <= argument of perigee <= 2 pi) ' oev(5) = right ascension of ascending node (radians) ' (0 <= raan <= 2 pi) ' oev(6) = true anomaly (radians) ' (0 <= true anomaly <= 2 pi) '''''''''''''''''''''''''''''''''''''''' local rmag, vmag, sma, p, q, const1 local h, xk, x1, y1, eccm, inc, xlambdat local raan, argper, tanom local rhat(3), vhat(3), hhat(3), vtmp(3), ecc(3) local hv(3), fhat(3), ghat(3) ' position and velocity magnitude rmag = vecmag(r()) vmag = vecmag(v()) ' position and velocity unit vectors uvector(r(), rhat()) uvector(v(), vhat()) ' angular momentum vectors vcross(r(), v(), hv()) uvector(hv(), hhat()) ' eccentricity vector for i% = 1 to 3 vtmp(i%) = v(i%) / mu next i% vcross(vtmp(), hv(), ecc()) for i% = 1 to 3 ecc(i%) = ecc(i%) - rhat(i%) next i% ' semimajor axis sma = 1.0 / (2.0 / rmag - vmag^2 / mu) ' keplerian period if (sma > 0.0) then period = pi2 / sqr(mu / sma^3) else period = 0.0 end if p = hhat(1) / (1.0 + hhat(3)) q = -hhat(2) / (1.0 + hhat(3)) const1 = 1.0 / (1.0 + p^2 + q^2) fhat(1) = const1 * (1.0 - p^2 + q^2) fhat(2) = const1 * 2.0 * p * q fhat(3) = -const1 * 2.0 * p ghat(1) = const1 * 2.0 * p * q ghat(2) = const1 * (1.0 + p^2 - q^2) ghat(3) = const1 * 2.0 * q h = vdot(ecc(), ghat()) xk = vdot(ecc(), fhat()) x1 = vdot(r(), fhat()) y1 = vdot(r(), ghat()) ' orbital eccentricity eccm = sqr(h * h + xk * xk) ' orbital inclination inc = 2.0 * atn(sqr(p * p + q * q)) ' true longitude xlambdat = atan3(y1, x1) ' check for equatorial orbit if (inc > 1.0e-10) then raan = atan3(p, q) else raan = 0.0 end if ' check for circular orbit if (eccm > 1.0e-10) then argper = modulo(atan3(h, xk) - raan) else argper = 0.0 end if ' true anomaly tanom = modulo(xlambdat - raan - argper) oev(1) = sma oev(2) = eccm oev(3) = inc oev(4) = argper oev(5) = raan oev(6) = tanom end sub '''''''''''''''''''''''''' '''''''''''''''''''''''''' sub orb2eci(oev(), r(), v()) ' convert orbital elements to state vector subroutine ' input ' oev(1) = semimajor axis (kilometers) ' oev(2) = orbital eccentricity (non-dimensional) ' oev(3) = orbital inclination (radians) ' (0 <= oev(3) <= pi) ' oev(4) = argument of perigee (radians) ' (0 <= oev(4) <= 2 pi) ' oev(5) = right ascension of ascending node or ' east longitude of ascending node (radians) ' (0 <= oev(5) <= 2 pi) ' oev(6) = true anomaly (radians) ' (0 <= oev(6) <= 2 pi) ' output ' r() = geocentric, equatorial position vector (kilometers) ' v() = geocentric, equatorial velocity vector (kilometers/second) ' special notes ' (1) eci coordinates if oev(5) = raan ' (2) ecf coordinates if oev(5) = east longitude ' (3) valid for elliptic and hyperbolic orbits ''''''''''''''''''''''''''''''''''''''''''''''' local sma, ecc, inc, argper, nangle, tanom local slr, rm, arglat, sarglat, carglat local sinc, cinc, snangle, cnangle, c1, c2, c3 ' "unload" orbital elements array sma = oev(1) ecc = oev(2) inc = oev(3) argper = oev(4) nangle = oev(5) tanom = oev(6) ' semiparameter slr = sma * (1.0 - ecc * ecc) ' radius magnitude rm = slr / (1.0 + ecc * cos(tanom)) ' argument of latitude arglat = argper + tanom sarglat = sin(arglat) carglat = cos(arglat) sinc = sin(inc) cinc = cos(inc) snangle = sin(nangle) cnangle = cos(nangle) c1 = sqr(mu / slr) c2 = ecc * cos(argper) + carglat c3 = ecc * sin(argper) + sarglat ' x, y and z components of the position vector r(1) = rm * (cnangle * carglat - snangle * cinc * sarglat) r(2) = rm * (snangle * carglat + cinc * sarglat * cnangle) r(3) = rm * sinc * sarglat ' x, y and z components of the velocity vector v(1) = -c1 * (cnangle * c3 + snangle * cinc * c2) v(2) = -c1 * (snangle * c3 - cnangle * cinc * c2) v(3) = c1 * c2 * sinc end sub ''''''''''''''''''' ''''''''''''''''''' sub svprint(r(), v()) ' print position and velocity vectors subroutine '''''''''''''''''''''''''''''''''''''''''''''''' local rmag, vmag print "position vector x-component ", str$(r(1), 0, -12) + " kilometers" print "position vector y-component ", str$(r(2), 0, -12) + " kilometers" print "position vector z-component ", str$(r(3), 0, -12) + " kilometers" rmag = vecmag(r()) print " " print "position vector magnitude ", str$(rmag, 0, -12) + " kilometers" print " " print "velocity vector x-component ", str$(v(1), 0, -12) + " kilometers/second" print "velocity vector y-component ", str$(v(2), 0, -12) + " kilometers/second" print "velocity vector z-component ", str$(v(3), 0, -12) + " kilometers/second" vmag = vecmag(v()) print " " print "velocity vector magnitude ", str$(vmag, 0, -12) + " kilometers/second" end sub '''''''''''''''' '''''''''''''''' sub oeprint(oev()) ' print classical orbital elements subroutine ''''''''''''''''''''''''''''''''''''''''''''' print "semimajor axis ", str$(oev(1), 0, 4) + " kilometers" print "eccentricity ", str$(oev(2), 0, 10) print "orbital inclination ", str$(rtd * oev(3), 0, 4) + " degrees" print "argument of perigee ", str$(rtd * oev(4), 0, 4) + " degrees" print "RAAN ", str$(rtd * oev(5), 0, 4) + " degrees" print "true anomaly ", str$(rtd * oev(6), 0, 4) + " degrees" print " " end sub ''''''''''''''''''''''''' ''''''''''''''''''''''''' 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 ''''''''''''''''''''''''''' ''''''''''''''''''''''''''' 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 vecmag(a()) as float ' vector magnitude function ' input ' { a } = column vector ( 3 rows by 1 column ) ' output ' vecmag = scalar magnitude of vector { a } vecmag = sqr(a(1) * a(1) + a(2) * a(2) + a(3) * a(3)) end function '''''''''''''''''''' '''''''''''''''''''' sub uvector (a(), b()) ' unit vector subroutine ' input ' a = column vector (3 rows by 1 column) ' output ' b = unit vector (3 rows by 1 column) ''''''''''''''''''''''''''''''''''''''' local i%, amag amag = vecmag(a()) for i% = 1 to 3 if (amag <> 0.0) then b(i%) = a(i%) / amag else b(i%) = 0.0 end if next i% end sub ''''''''''''''''''''''' ''''''''''''''''''''''' sub vcross(a(), b(), c()) ' vector cross product subroutine ' { c } = { a } x { b } ' input ' { a } = vector a ( 3 rows by 1 column ) ' { b } = vector b ( 3 rows by 1 column ) ' output ' { c } = { a } x { b } ( 3 rows by 1 column ) c(1) = a(2) * b(3) - a(3) * b(2) c(2) = a(3) * b(1) - a(1) * b(3) c(3) = a(1) * b(2) - a(2) * b(1) end sub '''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''' function vdot(a(), b()) as float ' vector dot product function ' c = { a } dot { b } ' input ' n% = number of rows ' { a } = column vector with n rows ' { b } = column vector with n rows ' output ' vdot = dot product of { a } and { b } '''''''''''''''''''''''''''''''''''''''' local i%, c = 0.0 for i% = 1 to 3 c = c + a(i%) * b(i%) next i% vdot = c end function '''''''''''''''''''''''' '''''''''''''''''''''''' function cotan(x) as float ' cotangent function ' x is input value '''''''''''''''''' cotan(x) = 1.0 / tan(x) end function |
||||
vegipete Guru Joined: 29/01/2013 Location: CanadaPosts: 1082 |
Ah, the remarkable march of technology. A few dollars worth of computer can calculate this sort of thing now. How many tons of computer at what cost at what power consumption did NASA need to do this 50 years ago? Visit Vegipete's *Mite Library for cool programs. |
||||
Print this page |