druzice.cc
Simple 3D model demo (not valid)
<![CDATA[
////////////////////////////////////////////////////////////////////////////
// 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]](druzice3D-1.png)
![[results]](druzice3D-2.png)
![[results]](druzice3D-3.png)
Last modification:
If you have some suggestions to improve this page, please mail to
peringer AT fit.vutbr.cz