 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 - Program NPOE
| Author | Message | ||||
| cdeagle Senior Member  Joined: 22/06/2014 Location: United StatesPosts: 261 |
Program NPOE (Numerical Prediction of Orbital Events) was originally developed using FORTRAN. During my career I ported NPOE to other computer languages such as MATLAB and QuickBASIC. Now it has been ported to MMBASIC for both the MicroMite eXtreme and the Raspberry Pi. From the introduction in the user manual; NPOE is an interactive computer program which can model important orbital events and predict the long term evolution of satellites in Earth orbits. Program NPOE implements a special perturbation solution of orbital motion using a variable step size Runge-Kutta-Fehlberg (RKF78) integration method to numerically integrate Cowell’s form of the system of differential equations. Orbital events are predicted using Brent’s method for finding the root of a nonlinear equation. The user can control both the integration and root-finding criteria. The computer program can also simply propagate a user-defined orbit with spacecraft mass properties and perturbations provided by the user. Let's talk about the MMBASIC MicroMite eXtreme/Raspberry Pi version of NPOE. This is a subset of the full up FORTRAN version. This version consists of 4419 lines of comments, source code, data statements and white space. It's quite a number cruncher and I recommend using a Pi whenever possible. Here's an example of an event output computed by NPOE. This simulation used an 8 by 8 Earth gravity model during the propagation while finding and displaying flight conditions whenever the geodetic sub-point latitude of the spacecraft was 20 degrees north. program npoe - numerical prediction of orbital events ----------------------------------------------------- degree of gravity model 8 order of gravity model 8 integration tolerance 1.00000e-08 root-finding tolerance 0.0001 initial conditions ------------------ calendar date April 18 2017 UTC 0 hours 0 minutes 0.00 seconds position vector x-component 5.506706064204e+03 kilometers position vector y-component 5.506706064204e+03 kilometers position vector z-component 0.000000000000e+00 kilometers position vector magnitude 7.787658400000e+03 kilometers velocity vector x-component -4.504405239130e+00 kilometers/second velocity vector y-component 4.504405239130e+00 kilometers/second velocity vector z-component 3.458731500160e+00 kilometers/second velocity vector magnitude 7.248596878486e+00 kilometers/second semimajor axis 8000.0000 kilometers eccentricity 0.0265427000 orbital inclination 28.5000 degrees argument of perigee 0.0000 degrees argument of latitude 0.0000 degrees RAAN 45.0000 degrees true anomaly 360.0000 degrees Keplerian period 118.6847 minutes nodal period 118.3855 minutes anomalistic period 118.5519 minutes geodetic latitude 0.0000 degrees geodetic altitude 1409.5221 kilometers flight path angle 0.0000 degrees right ascension 45.0000 degrees declination 0.0000 degrees east longitude 198.7002 degrees Here's the link to the zipped archive of the MMBASIC source code and a few other text files. 2017-04-24_175154_npoe.zip Both the FORTRAN and MATLAB user manuals can be downloaded from my Google Drive collection at https://drive.google.com/drive/folders/0Bx2YPzS0EoweUFR6ZHZ0cVpFWjA?usp=sharing The initial orbit and epoch, and other important information must be defined by the user in an MMBASIC subroutine named "sim_definition". Here are the contents of a typical definition file. Hopefully the annotation in the source code will help. sub sim_definition(ri(), vi(), target, fxnew, tfinal) ' define simulation inputs subroutine ''''''''''''''''''''''''''''''''''''' local sma, ecc, xinc, argper, raan, tanom local cmonth, cday, cyear, uthr, utmin, utsec local jdate, oev(32) '''''''''''''''''''''''''''''''''''' ' initial classical orbital elements '''''''''''''''''''''''''''''''''''' ' semimajor axis (kilometers; semimajor axis > 0) sma = 8000.0 ' orbital eccentricity (non-dimensional); 0 <= eccentricity < 1) ecc = 0.0265427 ' orbital inclination (degrees); 0 <= inclination <= 180) xinc = 28.5 ' argument of perigee (degrees); 0 <= argument of perigee <= 360) argper = 0.0 ' right ascension of the ascending node (degrees); 0 <= RAAN <= 360) raan = 45.0 ' true anomaly (degrees); 0 <= true anomaly <= 360) tanom = 0.0 oev(1) = sma oev(2) = ecc oev(3) = xinc * dtr oev(4) = argper * dtr oev(5) = raan * dtr oev(6) = tanom * dtr '''''''''''''''''''''''''''''''''''''''''''' ' request initial calendar date and UTC time '''''''''''''''''''''''''''''''''''''''''''' cmonth = 4 cday = 18 cyear = 2017 uthr = 0.0 utmin = 0.0 utsec = 0.0 ' simulation period (days) ndays = 0.25 ' simulation period (seconds) tfinal = 86400.0 * ndays '''''''''''''''''''''''''''''' ' algorithm control parameters '''''''''''''''''''''''''''''' ' integration error tolerance; a value of 1.0e-8 is recommended) tetol = 1.0e-8 ' root-finding error tolerance; a value between 1.0e-4 and 1.0e-6 is recommended) rtol = 1.0e-4 ''''''''''''''''''''''''''''''''''''''''''' ' request degree and order of gravity model ''''''''''''''''''''''''''''''''''''''''''' ' degree of the gravity model (zonals); 0 <= zonals <= 8) lgrav% = 8 ' order of the gravity model (tesserals); 0 <= tesserals <= 8) mgrav% = 8 ''''''''''''''''''''''' ' orbital perturbations ''''''''''''''''''''''' ' include solar perturbation (1 = yes, 0 = no) isun% = 1 ' include lunar perturbation (1 = yes, 0 = no) imoon% = 1 ' include drag perturbation (1 = yes, 0 = no) idrag% = 1 ' include srp perturbation (1 = yes, 0 = no) isrp% = 1 ' drag coefficient (non-dimensional) cd = 2.2 ' cross-sectional area (square meters) areadrag = 40.0 ' spacecraft mass (kilograms) scmass = 6078.0 ' compute ballistic coefficient bcoeff = 1.0e-6 * areadrag * cd /scmass ' reflectivity constant (non-dimensional) reflect = 1.5 ' cross-sectional area (square meters) areasrp = 10.0 ' spacecraft mass (kilograms) scmass = 6078.0 ' compute srp coefficient csrp0 = 0.000001 * reflect * pss * aunit^2 * areasrp / scmass ' determine initial eci state vector orb2eci(oev(), ri(), vi()) ' compute julian date at initial time julian(cmonth, cday, cyear, jdate) jdate0 = jdate + uthr / 24.0 + utmin / 1440.0 + utsec / 86400.0 ' compute greenwich sidereal time at initial utc (radians) gast2(jdate0, gst0) ' determine complete set of initial orbital elements eci2orb2(0, ri(), vi(), oev()) ''''''''''''''''''''''''''''''''''''''' ' select the orbital element to predict ''''''''''''''''''''''''''''''''''''''' ' 1 = geodetic altitude ' 2 = geodetic latitude ' 3 = east longitude ' 4 = geocentric declination ' 5 = true anomaly ' 6 = argument of latitude ' 7 = flight path angle ' 8 = orbital speed ' 9 = right ascension ' NOTE: to propagate without searching ' for events, set ittype% = 0 ittype% = 2 '''''''''''''''''''''''''''''''''' ' define target value and compute ' initial value of objective function ''''''''''''''''''''''''''''''''''''' select case ittype% case 0 ' propagate without searching for events case 1 ' geodetic altitude (kilometers) talt = 200.0 fxnew = oev(16) - talt case 2 ' geodetic latitude (degrees); -90 <= latitude <= +90) xlat = 20.0 tlat = xlat * dtr fxnew = oev(14) - tlat case 3 ' east longitude (degrees); 0 <= longitude <= 360) elong = 120.0 tlong = elong * dtr fxnew = oev(15) - tlong target = tlong case 4 ' geocentric declination (degrees); -90 <= declination <= +90) dec = 20.0 tdec = dec * dtr fxnew = oev(13) - tdec case 5 ' true anomaly (degrees); 0 <= true anomaly <= 360) tanom = 260.0 ttanom = tanom * dtr fxnew = oev(6) - ttanom target = ttanom case 6 ' argument of latitude (degrees); 0 <= arglat <= 360) arglat = 200.0 targlat = arglat * dtr fxnew = oev(8) - targlat target = targlat case 7 ' flight path angle (degrees); -90 <= fpa <= +90) fpa = 1.0 tfpa = fpa * dtr fxnew = oev(11) - tfpa case 8 ' orbital speed (kilometers/second) tspeed = 0.0 fxnew = oev(28) - tspeed case 9 ' right ascension (degrees); 0 <= right ascension <= 360) rasc = 100.0 trasc = rasc * dtr fxnew = oev(12) - trasc target = trasc end select end sub |
||||
 Print this page |
