druzice.cc


Simple 3D model demo (not valid)

//////////////////////////////////////////////////////////////////////////// // druzice3D.cc -- příklad použití knihovny SIMLIB pro spojitou simulaci 3D // // model družice pohybující se v soustavě Země-Měsíc // [upozornění - bez záruky] // // (c) 1997-2001 Petr Peringer #include "simlib.h" #include "simlib3D.h" // pojmenování typů: typedef Value3D Position; typedef Value3D Speed; typedef Value3D Force; // konstanty: const double gravity_constant = 6.67e-11; // gravitační konstanta const double m0 = 1000; // hmotnost družice [kg] const double m1 = 5.983e24; // hmotnost Země [kg] const double m2 = 7.374e22; // hmotnost Měsíce [kg] // inicializace: const Position p0(36.0e6, 0, 0); // poloha družice const Position p2(384.405e6, 0, 0); // poloha Měsíce const Speed v0(0, 4200, 1600); // počáteční rychlost družice const Speed v2(0, 1022, 0); // oběžná rychlost Měsíce const double YEAR = 3.6e3 * 24 * 365; // 1 rok v sekundách const double MAXTime = 1*YEAR; // doba simulace Constant3D Zero(0, 0, 0); // pomocný objekt //////////////////////////////////////////////////////////////////////////// struct MassPoint { double m; // hmotnost Expression3D inforce; // vstupní síla Integrator3D v; // rychlost Integrator3D p; // poloha MassPoint(const double mass, Position p0, Speed v0 = Speed(0, 0, 0)): m(mass), inforce(Zero), v(inforce / m, v0), p(v, p0) {} void SetInput(Input3D i) { inforce.SetInput(i); } }; //////////////////////////////////////////////////////////////////////////// struct MyWorld { enum { MAX = 10 }; // pozor: nekontroluje se překročení MassPoint *m[MAX]; // pole ukazatelů na objekty v systému int n; // počet objektů -- hmotných bodů MyWorld(); Position Mcenter() { // výpočet těžiště systému Position p = Value3D(0, 0, 0); double sum = 0; for(int i = 0; i < n; i++) { p = p + m[i]->m * m[i]->p.Value(); sum += m[i]->m; } return p / sum; } }; MyWorld *w; // svět vznikne až později :-) //////////////////////////////////////////////////////////////////////////// // gravitační síla pusobící na hmotný bod p // class GravityForce:public aContiBlock3D { MassPoint *p; MyWorld *w; public: GravityForce(MassPoint * _p, MyWorld * _w):p(_p), w(_w) {} Force Value() { Force f(0, 0, 0); // gravitační síla for (int i = 0; i < w->n; i++) { MassPoint *m = w->m[i]; if(m == p) continue; Value3D distance = m->p.Value() - p->p.Value(); double d = abs(distance); // vzdálenost f = f + (distance * gravity_constant * p->m * m->m / (d * d * d)); } return f; } }; //////////////////////////////////////////////////////////////////////////// // pohyblivé planety // typedef MassPoint Planet; typedef MassPoint Moon; typedef MassPoint Satellite; MyWorld::MyWorld(): n(0) { // vytvoříme digitální svět m[n++] = new Planet(m1, Position(0, 0, 0), Speed(0, 0, 0)); m[n++] = new Moon(m2, p2, v2); m[n++] = new Satellite(m0, p0, v0); for (int i = 0; i < n; i++) // zapneme silové působení m[i]->SetInput(new GravityForce(m[i], this)); } //////////////////////////////////////////////////////////////////////////// // sledování stavu modelu ... // extern void Sample(); Sampler s(Sample, 3600); // vzorkování stavu modelu void Sample() { Position t = w->Mcenter(); // těžiště systému Position p = (w->m[2])->p.Value() - t; Position p1 = (w->m[1])->p.Value() - t; Print("%g ", Time); Print("%g %g %g ", p1.x() / 1e6, p1.y() / 1e6, p1.z() / 1e6); Print("%g %g %g ", p.x() / 1e6, p.y() / 1e6, p.z() / 1e6); Print("\n"); } //////////////////////////////////////////////////////////////////////////// // popis experimentu ... // int main() { SetOutput("druzice3D.dat"); w = new MyWorld; Print("# Model oběhu družice kolem soustavy Země-Měsíc \n"); Init(0, 3.6e3*24*100); // inicializace experimentu SetMethod("abm4"); // vícekroková integrační metoda SetAccuracy(1e-7); // je nutná vysoká přesnost Run(); // start simulace } ////////////////////////////////////////////////////////////////////////////

Simple 3D model results (visualized using Gnuplot):

[results] [results] [results]


Last modification:

If you have some suggestions to improve this page, please mail to peringer AT fit.vutbr.cz