Orbiter-Forum Yet another propagation problem (atmospheric drag mostly)
 Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

 Math & Physics Mathematical and physical problems of space flight and astronomy.

 05-15-2018, 03:44 PM #1 Raccoon Orbinaut Yet another propagation problem (atmospheric drag mostly) Hello! I’ve write a program for orbit propagation and choose numerical integration method runge-kutta 4. I wrote only 2 perturbation effect: 2-body and atmospheric drag. Here is the code of main cycle: Code: ```Satellite sat(QVector3D(7000,0,500), QVector3D(0,9,0),0.0025, 3725.0); QVector< QVector3D > satellite; double h = 1.0; QDateTime dt = QDateTime(QDate(1983,9,23),QTime(0,0,0,0)); for (int j = 0.0; j < 100; j++) { double Mjd_UTC = MJD(dt); // 32.184 – UTC-TAI double Mjd_TT = Mjd_UTC + 32.184/86400.0; RK4step(h, sat, Mjd_TT); dt = dt.addDays(j); satellite.push_back(sat.loc); satellite.push_back(sat.vel); }``` My implementation of RK4: Code: ```void RK4step(const double& h, Satellite& y, const double &Mjd_TT) { QVector3D loc(y.loc); QVector3D vel(y.vel); QVector3D k1,k2,k3,k4, acceleration; // rotation matrix for satellite relative to Earth at time Mjd_TT // O.Montenbruck, E. Gill; Satellite Orbits - Models, Methods, and Applications QMatrix3x3 T = NutationMatrix(Mjd_TT) * PreccessionMatrix(MJD_J2000, Mjd_TT); k1 = partialStep(y, T, y.size, y.mass); k2 = partialStep(y + k1 * h/2, T, y.size, y.mass); k3 = partialStep(y + k2 * h/2, T, y.size, y.mass); k4 = partialStep(y + k3 * h, T, y.size, y.mass); acceleration += (k1 + (k2 + k3) * 2 + k4) / 6; y.vel += acceleration * h; y.loc += y.vel * h; }``` Partial calculation: Code: ```QVector3D partialStep(const Satellite &y, const QMatrix3x3 &T, const double& satSize, const double& satMass ) { QVector3D acceleration (0.0, 0.0, 0.0); double r = Magnitude(y.loc); // MU - Earth gravity = 398602 double temp = -MU / (r * r * r); QVector3D acc2Body = y.loc * temp; double Cd = 2.2; QVector3D accFromAir = AccelDrag(0.0, y.loc, y.vel, T, satSize, satMass, Cd); acceleration = acc2Body + accFromAir; return acceleration; }``` Implementation of atmospheric drag which I find in “O.Montenbruck, E. Gill; Satellite Orbits - Models, Methods, and Applications” : Code: ```QVector3D AccelDrag ( double Mjd_TT, const QVector3D &loc, const QVector3D &vel, const QMatrix3x3& T, double Area, double mass, double CD ) { // W3 – Earth rotation = 0.0000729211 rad/sec const QVector3D omega (0.0, 0.0, W3); double vel_abs, dens; QVector3D loc_tod, vel_tod; QVector3D vel_rel, acc_tod; QMatrix3x3 T_trp; // for converting to ICRF/EME2000 T_trp = T.transposed(); loc_tod = T * loc; vel_tod = T * vel; // velocity relative to Earth rotation vel_rel = vel_tod - QVector3D::crossProduct(omega, loc_tod); vel_abs = Magnitude2(vel_rel); // Harris-Priester atmospheric density model dens = Density_HP (Mjd_TT, loc_tod); // atmospheric acceleration acc_tod = -0.5*CD*(Area/mass)*dens*vel_abs*vel_rel; return T_trp * acc_tod; }``` With solo 2-body perturbation effect orbit change correctly to check this I find other several implementations of it, for example matlab simple propagator of this guy: https://github.com/komrad36/SimpleOrbitalPropagator I know that in matlab it’s dorman-prince method but difference in values only -+1-5, it’s okay for me.The problem is that with solo atmospheric drag satellite position change not quit correct, I guess. Following common sense satellite must go down sooner or later, but it fly away from Earth all the time position vector just increase. And I don’t find simple example to compare results with mine. May be I am not right and propagate with only atmospheric drag is incorrect? Or it's error in code?
 05-15-2018, 05:31 PM #2 martins Orbiter Founder What do you mean by incorrect? You are free to simulate whatever you like. And yes, adding drag should reduce the orbital energy and lead to a decaying orbit. If yours doesn't then there must be a bug. Since you say your 2-body propagation works, the problem is likely in your drag term. Did you try to debug if it produces reasonable results? For example, is the vector returned by AccelDrag approx. opposite your y.vel value, as expected? In the first instance, you can probably drop the second term in the calculation of your vel_rel value (i.e. non-rotating planet). It removes one potential source of error, and should still produce the expected result. Are you sure that your definition of matrix T is as intended (i.e. transforms from global to local planet system and not the other way round)?
 05-15-2018, 10:48 PM #3 cristiapi Orbinaut Quote: Originally Posted by Raccoon  The problem is that with solo atmospheric drag satellite position change not quit correct, I guess. Following common sense satellite must go down sooner or later, but it fly away from Earth all the time position vector just increase. It usually happens when the time step is too big. I see: for (int j = 0.0; j < 100; j++) { ... dt = dt.addDays(j); } but I think that you should do: for (int j = 0; j < 100; j++) { ... dt = dt.addDays(t_step); } where t_step is a small constant time step.
 05-16-2018, 03:47 PM #4 Raccoon Orbinaut martins, i understand that i can simulate whatever i want, but under incorrect i mean that result will be not "good" for understanding theme. I am programmer not physicist, i need to understand how propagation work on simple example. With 2-Body was simple but other perturbation effects are hard to calculate. I want understand them 1-by-1. Vector returned by AccelDrag is opposite to my y.vel value, but power is E-20 or E-17. Values are too low and with my example where velocity is (0, 9, 0) after 2500 iteration it will be (-1E-14, 7.996... , -2E-13). Functions for calculation matrix T i get from this book "O.Montenbruck, E. Gill; Satellite Orbits - Models, Methods, and Applications". May be they are not what i needed. In different sources all this coordinates systems call by different names, it's confusing for me. I even found atmospheric models (US standart and Russian GOST) which doesn't depend on time. And still same results! cristiapi, i try all variants even crazy ones to change years, but again values which i get too low.
 05-18-2018, 02:30 PM #5 Raccoon Orbinaut martins, vector returned by AccelDrag is opposite to mine y.vel value but values are too small from E-20 to E-17. And their influence is small too. about matrix T i am not sure, but i take it from book "O.Montenbruck, E. Gill; Satellite Orbits - Models, Methods, and Applications", may be i miss something. cristiapi, there is nothing change when i take another step, i even try years and get same results.
 05-24-2018, 10:25 PM #6 cristiapi Orbinaut Quote: Originally Posted by Raccoon  cristiapi, there is nothing change when i take another step, i even try years and get same results. Did you try 1 or 2 seconds steps? If you use a low order integrator (like RK4) the step must be very small for a LEO satellite.
 05-25-2018, 11:40 AM #7 martins Orbiter Founder Quote: Originally Posted by Raccoon  martins, vector returned by AccelDrag is opposite to mine y.vel value but values are too small from E-20 to E-17. And their influence is small too. about matrix T i am not sure, but i take it from book "O.Montenbruck, E. Gill; Satellite Orbits - Models, Methods, and Applications", may be i miss something. How do you know that the values are too small? Surely this depends on altitude. You didn't specify your orbit altitude. What are the correct values for the case you are debugging? Also, this seems to be at odds with your previous statement that your propagator works ok for the frictionless 2-body case, and only has problems when including the atmospheric drag term. If the drag deceleration is in the correct direction, and negligibly small anyway, I don't see how it can lead to an increase in orbital energy.
 05-25-2018, 12:18 PM #8 cristiapi Orbinaut Quote: Originally Posted by martins  You didn't specify your orbit altitude. What are the correct values for the case you are debugging? The position vector is QVector3D(7000,0,500), which means that the altitude is about 647 km. According to the NRLMSISE-00 model, the density is about 10-14 kg/m3.

 Tags atmospheric drag, c++, propagation

 Posting Rules BB code is On Smilies are On [IMG] code is On HTML code is Off You may not post new threads You may not post replies You may not post attachments You may not edit your posts
 Forum Jump User Control Panel Private Messages Subscriptions Who's Online Search Forums Forums Home Orbiter-Forum.com     Announcements     Meets & Greets Orbiter Space Flight Simulator     Orbiter Web Forum         OFMM         Orbiter Forum Space Station         Simpit Forum     General Questions & Help     MFD Questions & Help     Hardware & Software Help     Tutorials & Challenges     Orbiter SDK     Orbiter Visualization Project     Orbiter Beta » Orbiter Project Orbiter Addons     OrbitHangar Addons & Comments     Addons     Addon Development     Addon Requests     Addon Support & Bugs         Addon Developer Forums             Project Apollo - NASSP     Orbiter Lua Scripting Far Side of the Moon     Spaceflight News     Math & Physics     Astronomy & the Night Sky     Backyard Rocketry     Brighton Lounge     International Forum

All times are GMT. The time now is 02:20 AM.