太阳系(重力)仿真

instemast · 发表于 2010-12-20 15:56:00

- 请按[N]适当加速(减速则按[M])
- 我放置了一个探测器，自地球附近发射，并且中途受月球引力swing by而加速
- 如果观察探测器，那么按[4]切换'主角'为地球，然后按[A]放大视图(缩小则按[Z])，并且别加速太大
- 详细的按键说明在窗口标题栏上有

//Test02.cpp
#include "stdafx.h"

using namespace std;
using namespace inst;
using namespace inst::i2d;

typedef double2 Vec;

bool g_bExit = false;
CGameWnd *g_pWnd;
CBuffer *g_pBuf;
CFont *g_pFont;

void Wnd_WndProc(UINT msg,WPARAM wparam,LPARAM lparam);
void InitLogic();
void Update1();
void Flush1();
void Update2();
void Flush2();
void OnReset(){}
void Init();
void Destroy(){}
bool MsgLoop();

void Init()
{
g_pWnd = new CGameWnd(L"[keys] A/Z/Q: zoom; N/M/B: time zoom; 0-9: follow with a star(1-9, 0:reset)");
g_pWnd->SetCallback(&Wnd_WndProc);
g_pWnd->Show();
g_pWnd->Activate();

CGraphics::CreateGraphics(800, 600, g_pWnd->GetHWnd());
g_pBuf = CGraphics::GetGraphics()->GetBackBuffer();
g_pFont = new CFont(L"Arial");

InitLogic();
}

double zoom = 1.0;
double t_zoom = 1.0;
int follow = -1;

void Wnd_WndProc(UINT msg,WPARAM wparam,LPARAM lparam)
{
Switch (msg)
{
Case WM_KEYDOWN:
      if (wparam=='A')
         zoom *= 2;
      else if (wparam=='Z')
         zoom /= 2;
      else if (wparam=='Q')
         zoom = 1.0;
      else if (wparam=='N')
         t_zoom *= 2;
      else if (wparam=='M')
         t_zoom /= 2;
      else if (wparam=='B')
         t_zoom = 1.0;
      else if (wparam >= '0' && wparam <= '9')
         follow = FromWString<int>(ToWString(wchar(wparam))) - 1;
      break;
Case WM_LBUTTONDOWN:
      break;
Case WM_DESTROY:
      g_bExit=true;
      break;
Case WM_DISPLAYCHANGE:
      CGraphics::GetGraphics()->Reset();
      OnReset();
      break;
}
}

int APIENTRY WinMain(HINSTANCE hInstance,HINSTANCE hPrevInstance,LPSTR lpCmdLine,int nCmdShow)
{
try
{
      Init();

      While (1)
         if(MsgLoop()) break;

      Destroy();

      return 0;
}
catch(...)
{
}
}

double g_dTime1 = 0;
double g_dTime2 = 0;

bool MsgLoop()
{
if(g_bExit) return true;

::MSG msg = {};
while(PeekMessage(&msg, NULL, 0, 0, PM_REMOVE))
{
      TranslateMessage(&msg);
      DispatchMessage(&msg);
}

const int N = 5000;

static int t1 = timeGetTime();
int t2 = timeGetTime();
if (t2 - t1 < 1)
{
      While (t2 - t1 < 1)
      {
         ::Sleep(1);
         t2 = timeGetTime();
      }
}
g_dTime2 = MIN(t2-t1, 100) / 1000.0;
g_dTime1 = g_dTime2 / N;
t1 = t2;

for (int i=0; i<N; i++)
{
      Update1();
      Flush1();
}

Update2();
Flush2();

return false;
}

/////////////////////////////////////////////////////////////////

const double _GameSecPerRealSec = 1;
double GameSPerRealS() { return _GameSecPerRealSec * t_zoom; }
double t = 0;
double dt() { return g_dTime1 * GameSPerRealS(); }

__forceinline byte3 MyGetRgb(dword color){ return byte3(GetR(color), GetG(color), GetB(color)); }
__forceinline dword MyRgb(byte3 rgb){ return Rgb(rgb(0), rgb(1), rgb(2)); }

const double KmPerAU = 149597870.691; // km
const double G = 6.67300e-20; // km3*kg-1*s-2

const double m_sun = 1.98892e30; // kg

const double m_mercury = 3.3022e23; // kg
const double d_mercury = 57910000; // avg, km
const double v_mercury = 47.8725; // avg, km/s

const double m_venus = 4.8685e24; // kg
const double d_venus = 108208930; // avg, km
const double v_venus = 35.0214; // avg, km/s

const double m_earth = 5.9742e24; // kg
const double d_earth = 149597870.7; // avg, km
const double v_earth = 29.780; // avg, km/s

const double m_moon = 7.347673e22; // kg
const double d_moon = 384400; // avg, km
const double v_moon = 1.022; // avg, km/s

const double m_mars = 6.4185e23; // kg
const double d_mars = 227936640; // avg, km
const double v_mars = 24.1309; // avg, km/s

class CObject;

vector<SmartPtr<CObject>> AllObjects;

class CObject : public virtual IDynamic
{
INST_DYNAMIC1(L"CObject");
public:
Vec s, v;
double m;

CObject(Vec _s, Vec _v, double _m) { s = _s; v = _v; m = _m; }

Vec f_g(const CObject *o) const
{
      INST_ASSERT(o != this);
      Vec D_s = o->s - s;
      double norm_D_s = Norm(D_s);
      return G * o->m * D_s / (norm_D_s * norm_D_s * norm_D_s);
}

Vec g_total() const
{
      Vec ret(0, 0);
      for(int i=0; i<AllObjects.size(); i++)
         if (this != AllObjects)
            ret += f_g(AllObjects);
      return ret;
}

virtual Vec F_others_total() const { return Vec(0, 0); }

Vec a_total() const { return g_total() + F_others_total() / m;}

Vec v_next() const { return v + (a_total()) * dt(); }
Vec s_next() const { return s + v * dt(); }

Vec _s_next;
Vec _v_next;

virtual void Update()
{
      //_s_next = s_next();
      //_v_next = v_next();
      s = s_next();
      v = v_next();
}

virtual void Flush()
{
      s = _s_next;
      v = _v_next;
}
};

class CObjectX : virtual public CObject
{
INST_DYNAMIC2(L"CObjectX", CObject);
public:
wstring name;
int radius;
dword color;

CObjectX(wstring _name, Vec _s, Vec _v, double _m, int _radius, dword _color)
      : CObject(_s, _v, _m)
{
      name = _name;
      radius = _radius;
      color = _color;
}
};

double Fps1() { return 1.0 / g_dTime1; }
double Fps2() { return 1.0 / g_dTime2; }

Vec Dir(double angle) { return Vec(cos(angle), sin(angle)); }

void InitLogic()
{
SmartPtr<CObjectX> pSun = new CObjectX(
      L"sun",
      Vec(0, 0),
      Vec(0, 0),
      m_sun,
      15,
      0xFFFFFF);
AllObjects.push_back(pSun.GetPtr());

SmartPtr<CObjectX> pMercury = new CObjectX(
      L"mercury",
      pSun->s + d_mercury * Dir(0),
      pSun->v + v_mercury * Dir(-PI/2),
      m_mercury,
      5,
      0x00FFFF);
AllObjects.push_back(pMercury.GetPtr());

SmartPtr<CObjectX> pVenus = new CObjectX(
      L"venus",
      pSun->s + d_venus * Dir(0),
      pSun->v + v_venus * Dir(-PI/2),
      m_venus,
      8,
      0xFFBB00);
AllObjects.push_back(pVenus.GetPtr());

SmartPtr<CObjectX> pEarth = new CObjectX(
      L"earth",
      pSun->s + d_earth * Dir(0),
      pSun->v + v_earth * Dir(-PI/2),
      m_earth,
      8,
      0x0000FF);
AllObjects.push_back(pEarth.GetPtr());

SmartPtr<CObjectX> pMoon = new CObjectX(
      L"moon",
      pEarth->s + d_moon * Dir(0),
      pEarth->v + v_moon * Dir(-PI/2),
      m_moon,
      2,
      0xFFFF00);
AllObjects.push_back(pMoon.GetPtr());

AllObjects.push_back(INST_NEW(CObjectX,
      L"tmp",
      pEarth->s + 8000 * Dir(0),
      pEarth->v + 9.91 * Dir(-PI/2),
      10.0,
      1,
      0x00FF00).GetPtr());

SmartPtr<CObjectX> pMars = new CObjectX(
      L"mars",
      pSun->s + d_mars * Dir(0),
      pSun->v + v_mars * Dir(-PI/2),
      m_mars,
      5,
      0xFF0000);
AllObjects.push_back(pMars.GetPtr());

}

void Update1()
{
for (int i=0; i<AllObjects.size(); i++)
      AllObjects->Update();
}

void Flush1()
{
t += dt();

return;

for (int i=0; i<AllObjects.size(); i++)
      AllObjects->Flush();
}

/////////////////////////////////////////////////////////////////

Vec s_view() { return (follow >= 0 && follow < AllObjects.size()) ? AllObjects[follow]->s : Vec(0, 0); }

Vec WorldToView(Vec world) { return world - s_view(); }
Vec ViewToWorld(Vec view) { return view + s_view(); }

const double _AUPerBufH = 3.5;
double KmPerBufH() { return _AUPerBufH * KmPerAU / zoom; }
double KmPerBufPx() { return KmPerBufH() / g_pBuf->GetH(); }

Vec ViewToBuf(Vec world) { return Vec(double(world(0)/KmPerBufPx() + g_pBuf->GetW()*0.5), double(g_pBuf->GetH()*0.5 - world(1)/KmPerBufPx())); }
Vec BufToView(Vec buf) { return Vec(double(buf(0)*KmPerBufPx() - g_pBuf->GetW()*0.5), double(- buf(1)*KmPerBufPx() - g_pBuf->GetH()*0.5)); }

void RenderStars()
{
for (int i = 0; i < AllObjects.size(); i++)
{
      SmartPtr<CObjectX> pTmp = inst_cast<CObjectX>(AllObjects);
      if (pTmp)
      {
         Vec s_buf = ViewToBuf(WorldToView(AllObjects->s));
         g_pBuf->FillEllipse(s_buf(0), s_buf(1), pTmp->radius, pTmp->radius, pTmp->color);
         g_pBuf->DrawText(pTmp->name, s_buf(0) - g_pFont->GetStringWidth(pTmp->name) / 2, s_buf(1) - pTmp->radius - 20, g_pFont, 0x00FF00);
      }
}
}

void Update2()
{
g_pBuf->FillRect(0,0,g_pBuf->GetW(),g_pBuf->GetH(),Rgb(0,0,0));

RenderStars();

// coord axises
g_pBuf->DrawLine(g_pBuf->GetW()/2, 0, g_pBuf->GetW()/2, g_pBuf->GetH(), Rgb(0,255,0));
g_pBuf->DrawLine(0, g_pBuf->GetH()/2, g_pBuf->GetW(), g_pBuf->GetH()/2, Rgb(0,255,0));
g_pBuf->DrawLine(g_pBuf->GetW()/2, 0, g_pBuf->GetW()/2-8, 16, Rgb(0,255,0));
g_pBuf->DrawLine(g_pBuf->GetW()/2, 0, g_pBuf->GetW()/2+8, 16, Rgb(0,255,0));
g_pBuf->DrawLine(g_pBuf->GetW(), g_pBuf->GetH()/2, g_pBuf->GetW()-16, g_pBuf->GetH()/2-8, Rgb(0,255,0));
g_pBuf->DrawLine(g_pBuf->GetW(), g_pBuf->GetH()/2, g_pBuf->GetW()-16, g_pBuf->GetH()/2+8, Rgb(0,255,0));

g_pBuf->DrawText(wstring(L"Kernel FPS = ") + ToWString((int)(Fps1())), 2, 0, g_pFont, Rgb(0,255,0), 1);
g_pBuf->DrawText(wstring(L"Shell FPS = ") + ToWString((int)(Fps2())), 2, 16, g_pFont, Rgb(0,255,0), 1);
g_pBuf->DrawText(wstring(L"Space Scale = 100% : ") + ToWString(KmPerBufH()) + L" m (" + ToWString(KmPerBufH()/KmPerAU) + L" AU)", 2, 32, g_pFont, Rgb(0,255,0), 1);
g_pBuf->DrawText(wstring(L"Time Scale = 1s : ") + ToWString(GameSPerRealS()) + L"s (" + ToWString((double)GameSPerRealS()/3600/24) + L" day)", 2, 48, g_pFont, Rgb(0,255,0), 1);
g_pBuf->DrawText(wstring(L"Time = ") + ToWString(t) + L"s (" + ToWString(t/3600/24) + L" day)", 2, 64, g_pFont, Rgb(0,255,0), 1);

//--------------------------------------------------------------------

CGraphics::GetGraphics()->Present();
}

void Flush2()
{
}

thejinchao · 发表于 2010-12-20 20:37:00

难道是传说中的三体模拟算法？貌似最近很火的一个东西，:-)

-----------------
欢迎光临我的博客 http://www.thecodeway.com

instemast · 发表于 2010-12-20 20:46:00

thejinchao: Re:太阳系(重力)仿真

难道是传说中的三体模拟算法？貌似最近很火的一个东西，:-)

----------------...

多个物体的数值仿真并不涉及解三体或多体问题。

Miu.C · 发表于 2010-12-20 20:53:00

已经是完整的代码了？

instemast · 发表于 2010-12-20 23:28:00

Miu.C: Re:太阳系(重力)仿真

已经是完整的代码了？

除非你认为要把类库代码包含进来才算完整

小小C · 发表于 2010-12-21 09:09:00

instemast: Re: Re:太阳系(重力)仿真

除非你认为要把类库代码包含进来才算完整

inst::i2d是啥?

instemast · 发表于 2010-12-21 19:05:00

小小C: Re: Re: Re:太阳系(重力)仿真

inst::i2d是啥?

封装的DDraw:

class CGraphics;
class C2DErr;
class CSurfBase;
class ISurfCreator;
class CSurf;
class CBuffer;
class CFont;
class CSimpleSurfCreator;

instemast · 发表于 2011-4-2 01:56:00

话说，这两句展示了真正的functional

virtual Vec F_others_total() const { return Vec(0, 0); }
Vec a_total() const { return g_total() + F_others_total() / m;}

instemast · 发表于 2016-12-17 17:30:14

又整了个3D版的：

// Test03.cpp : メインプロジェクトファイルです。

#include "stdafx.h"

namespace Test03 {

      using namespace std;
      using namespace takaku;

      R dTimeLogic = 0;
      R dTimeShell = 0;

      R TimeZoom = 10000.0;

      R t = 0;
      R dt;

      vector<SmartPtr<CObject>> AllObjects;

      R3 DirXY(R angle) { return R3(cos(angle), sin(angle), 0.0); }

      void InitLogic()
      {
            SmartPtr<CObjectX> pSun = new CObjectX(
                     L"sun",
                     R3(0, 0, 0),
                     R3(0, 0, 0),
                     m_sun,
                     15,
                     0xFFFFFF);
            AllObjects.push_back(pSun.GetPtr());

            SmartPtr<CObjectX> pMercury = new CObjectX(
                     L"mercury",
                     pSun->s + d_mercury * DirXY(0),
                     pSun->v + v_mercury * DirXY(-PI/2),
                     m_mercury,
                     5,
                     0x00FFFF);
            AllObjects.push_back(pMercury.GetPtr());

            SmartPtr<CObjectX> pVenus = new CObjectX(
                     L"venus",
                     pSun->s + d_venus * DirXY(0),
                     pSun->v + v_venus * DirXY(-PI/2),
                     m_venus,
                     8,
                     0xFFBB00);
            AllObjects.push_back(pVenus.GetPtr());

            SmartPtr<CObjectX> pEarth = new CObjectX(
                     L"earth",
                     pSun->s + d_earth * DirXY(0),
                     pSun->v + v_earth * DirXY(-PI/2),
                     m_earth,
                     8,
                     0x0000FF);
            AllObjects.push_back(pEarth.GetPtr());

            SmartPtr<CObjectX> pMoon = new CObjectX(
                     L"moon",
                     pEarth->s + d_moon * DirXY(0),
                     pEarth->v + v_moon * DirXY(-PI/2),
                     m_moon,
                     2,
                     0xFFFF00);
            AllObjects.push_back(pMoon.GetPtr());

            AllObjects.push_back(TAKAKU_NEW(CObjectX,
                     L"tmp",
                     pEarth->s + 8000 * DirXY(0),
                     pEarth->v + 9.91 * DirXY(-PI/2),
                     10.0,
                     1,
                     0x00FF00).GetPtr());

            SmartPtr<CObjectX> pMars = new CObjectX(
                     L"mars",
                     pSun->s + d_mars * DirXY(0),
                     pSun->v + v_mars * DirXY(-PI/2),
                     m_mars,
                     5,
                     0xFF0000);
            AllObjects.push_back(pMars.GetPtr());

return;
            const int N = 4;

            for (int i=-N; i<=N; i++)
                     for (int j=-N; j<=N; j++)
                              for (int k=-N; k<N; k++)
                              AllObjects.push_back(new /*CElecObjectX*/CObjectX(
                                    L"",
                                    R3(j * KmPerAU / 1, i * KmPerAU / 1, k * KmPerAU / 1),
                                    R3(),
                                    m_sun,
                                    /*1e28 * (rand()%2 == 0 ? -1 : 1),*/
                                    1,
                                    0xFFFFFF));

      /*
            R tmp = atan(2.0) - PI / 4;

            AllObjects.push_back(new CObjectX(
                     L"A",
                     R3(0.0, d_earth),
                     sqrt(2.0 / sqrt(5.0/2.0)) * v_earth * R3(cos(PI+tmp), sin(PI+tmp)),
                     m_sun,
                     5,
                     0xFFFFFF));

            AllObjects.push_back(new CObjectX(
                     L"B",
                     R3(d_earth, 0.0),
                     R3(0.0, sqrt(2.0) * v_earth),
                     m_sun,
                     5,
                     0xFFFFFF));

            AllObjects.push_back(new CObjectX(
                     L"C",
                     R3(0.0, - d_earth),
                     sqrt(2.0 / sqrt(5.0/2.0)) * v_earth * R3(cos(-tmp), sin(-tmp)),
                     m_sun,
                     5,
                     0xFFFFFF));
      */
            //::GetPrivateProfileIntW();
            //::GetPrivateProfileStringW();
            //::GetPrivateProfileSectionW();
            //::GetPrivateProfileSectionNamesW();
      }

      void UpdateLogic()
      {
            dt = dTimeLogic * TimeZoom;

            for (int i=0; i<AllObjects.size(); i++)
                     AllObjects->Update();
      }

      void FlushLogic()
      {
            t += dt;

            //for (int i=0; i<AllObjects.size(); i++)
            //       AllObjects->Flush();
      }

      R FovDeg = 75.0;

      int Stands = -1, Looks = -1;

      R3 ChildCamPos = ChildCamPos0, ChildCamFront = ChildCamFront0, ChildCamUp = ChildCamUp0;

      void Init()
      {
            G::Form1 = gcnew Form1();
            G::Form1->Show();

            PresentParameters ^pp = gcnew PresentParameters();
      pp->Windowed = true;
            pp->SwapEffect = SwapEffect::Copy;
      pp->BackBufferCount = 1;
      pp->BackBufferWidth = 800;
      pp->BackBufferHeight = 600;
            pp->AutoDepthStencilFormat = DepthFormat::D32;
            pp->DeviceWindow = G::Form1->pictureBox1;

            G::Dev = gcnew Device(0, DeviceType::Hardware, G::Form1->pictureBox1, CreateFlags::HardwareVertexProcessing,
                     gcnew array<PresentParameters^> { pp });

            G::RT = gcnew Texture(G::Dev, 800, 600, 0, Usage::RenderTarget, Format::X8R8G8B8, Pool::Default);

            G::VBRT = gcnew VertexBuffer(VertRT::typeid, 4, G::Dev, Usage::None,
                     VertexFormats::PositionW | VertexFormats::Texture1,
                     Pool::Managed);

            G::VBRT->SetData(
                     gcnew array<VertRT>(4) {
                              VertRT(Vector2(0, 0), 0, Vector2(0, 0)),
                              VertRT(Vector2(800.0f, 0), 0, Vector2(1.0f, 0)),
                              VertRT(Vector2(0, 600.0f), 0, Vector2(0, 1.0f)),
                              VertRT(Vector2(800.0f, 600.0f), 0, Vector2(1.0f, 1.0f)), },
                     0,
                     LockFlags::Discard);

            G::VBLocus = gcnew VertexBuffer(VertLocus::typeid, 1, G::Dev, Usage::None,
                     VertexFormats::Position | VertexFormats::Diffuse | VertexFormats::PointSize,
                     Pool::Managed);

            G::VBLocus->SetData(gcnew VertObject(Vector3()/*, 1.0f, 0xFF000000*/), 0, LockFlags::Discard);

            G::VBAxises = gcnew VertexBuffer(VertAxis::typeid, 6, G::Dev, Usage::None,
                     VertexFormats::Position | VertexFormats::Diffuse,
                     Pool::Managed);

            G::VBAxises->SetData(
                     gcnew array<VertAxis>(6) {
                              VertAxis(Vector3(-100 * KmPerAU, 0, 0), 0xFFFF0000),
                              VertAxis(Vector3(100 * KmPerAU, 0, 0), 0xFFFF0000),
                              VertAxis(Vector3(0, -100 * KmPerAU, 0), 0xFF00FF00),
                              VertAxis(Vector3(0, 100 * KmPerAU, 0), 0xFF00FF00),
                              VertAxis(Vector3(0, 0, -100 * KmPerAU), 0xFF0000FF),
                              VertAxis(Vector3(0, 0, 100 * KmPerAU), 0xFF0000FF), },
                     0,
                     LockFlags::Discard);

            InitLogic();

            G::VBObject = gcnew VertexBuffer(VertObject::typeid, 1, G::Dev, Usage::None, VertexFormats::Position, Pool::Managed);

            G::VBObject->SetData(
                     gcnew array<VertObject>(6) {
                              VertObject(Vector3(0, 1, 0)),
                              VertObject(Vector3(1, 0, 0)),
                              VertObject(Vector3(0, 0, 1)),
                              VertObject(Vector3(-1, 0, 0)),
                              VertObject(Vector3(0, 0, -1)),
                              VertObject(Vector3(0, -1, 0)),},
                     0,
                     LockFlags::Discard);

            G::IBObject = gcnew IndexBuffer(Int16::typeid, 24, G::Dev, Usage::None, Pool::Managed);

            /*G::IBRT->SetData(
                     gcnew array<Int16>(24){
                              0,1,2,       5,2,1,
                              0,2,3,       5,3,2,
                              0,3,4,       5,4,3,
                              0,4,1,       5,1,4
                     },
                     0, LockFlags::Discard);*/
      }

      void RenderObjects(bool bLocus)
      {
            for (int i = 0; i < AllObjects.size(); i++)
            {
                     G::Dev->Transform->World = DXMatrix::Translation(Vector3(DE_R3(AllObjects->s)));

                     if (bLocus)
                     {
                              G::Dev->Material.Diffuse = Color::FromArgb(0xFF | dynamic_cast<CObjectX *>(AllObjects.GetPtr())->color);

                              G::Dev->RenderState->DiffuseMaterialSource = ColorSource::Material;
                     }

                     if (bLocus)
                     {
                              G::Dev->SetStreamSource(0, G::VBLocus, 0);
                              G::Dev->VertexFormat = G::VBLocus->Description.VertexFormat;
                              G::Dev->DrawPrimitives(PrimitiveType::PointList, 0, 1);
                     }
                     else
                     {
                              G::Dev->SetStreamSource(0, G::VBObject, 0);
                              G::Dev->Indices = G::IBObject;
                              G::Dev->VertexFormat = G::VBObject->Description.VertexFormat;
                              G::Dev->DrawIndexedPrimitives(PrimitiveType::TriangleList, 0, 0, 24, 0, 8);
                     }

                     if (bLocus)
                     {
                              G::Dev->RenderState->DiffuseMaterialSource = ColorSource::Color1;
                     }
            }
      }

      void UpdateShell()
      {
            G::Dev->RenderState->CullMode = Cull::None;
            G::Dev->RenderState->Lighting = false;
            G::Dev->RenderState->DiffuseMaterialSource = ColorSource::Color1;

            G::Dev->Clear(ClearFlags::Target, 0, 0, 0);

            G::Dev->BeginScene();

            G::Dev->Transform->View = DXMatrix::LookAtLH(ToDXVec(CamPos()), ToDXVec(CamPos() + CamFront()), ToDXVec(CamUp()));
            G::Dev->Transform->Projection = DXMatrix::PerspectiveFovLH(Fov(), 4.0f / 3, KmPerAU / 50, KmPerAU * 200);

            G::Dev->SetRenderTarget(0, G::RT->GetSurfaceLevel(0));

            RenderObjects(true);

            G::Dev->SetRenderTarget(0, G::Dev->GetBackBuffer(0, 0, BackBufferType::Mono));

            //G::Dev->TextureState[0]
            G::Dev->SetStreamSource(0, G::VBRT, 0);
            G::Dev->VertexFormat = G::VBRT->Description.VertexFormat;
            G::Dev->DrawPrimitives(PrimitiveType::TriangleFan, 0, 2);

            RenderObjects(false);

            G::Dev->Transform->World = DXMatrix::Identity;

            G::Dev->SetStreamSource(0, G::VBAxises, 0);
            G::Dev->VertexFormat = G::VBAxises->Description.VertexFormat;
            G::Dev->DrawPrimitives(PrimitiveType::LineList, 0, 3);

            G::Dev->EndScene();

            G::Dev->Present();

            G::Form1->uiLooks->Text = Looks.ToString();
            G::Form1->uiStands->Text = Stands.ToString();
            G::Form1->uiFov->Text = FovDeg.ToString();
            G::Form1->uiTimeZoom->Text = TimeZoom.ToString();

            G::Form1->uiInfo->Text = "Kernel FPS: " + FpsLogic() + "\r\n"
                     + "Shell FPS: " + FpsShell();
      }

      void Destroy()
      {
      }

}

using namespace Test03;

[STAThreadAttribute]
int main(array<System::String ^> ^args)
{
      // コントロールが作成される前に、Windows XP ビジュアル効果を有効にします
      Application::EnableVisualStyles();
      Application::SetCompatibleTextRenderingDefault(false);

      // メインウィンドウを作成して、実行します
      //Application::Run(gcnew Form1());

      Init();

      while (G::Form1->Visible)
      {
            Application::DoEvents();

            const int N = 1;

            static int t1 = timeGetTime();
            int t2 = timeGetTime();
            if (t2 - t1 < 1)
            {
                     while (t2 - t1 < 1)
                     {
                              ::Sleep(1);
                              t2 = timeGetTime();
                     }
            }
            dTimeShell = MIN(t2-t1, 100) / 1000.0;
            dTimeLogic = dTimeShell / N;
            t1 = t2;
            for (int i=0; i<N; i++)
            {
                     UpdateLogic();
            }
            UpdateShell();
      }

      Destroy();

      return 0;
}

instemast · 发表于 2016-12-17 17:30:37

// Test03.h

#pragma once

namespace Test03 {

      using namespace std;
      using namespace takaku;

      typedef double                      R;
      typedef Matrix<R,3,1>       R3, R31;
      typedef Matrix<R,1,3>       R13;
      typedef Matrix<R,3,3>       R33;
      #define DE_R3(X)       (X)(0), (X)(1), (X)(2)

      inline R3 Cross(const R3& a, const R3& b) { return R3(a(1,0) * b(2,0) - a(2,0) * b(1,0), a(2,0) * b(0,0) - a(0,0) * b(2,0), a(0,0) * b(1,0) - a(1,0) * b(0,0)); }

      extern R dTimeLogic, dTimeShell;
      static R FpsLogic() { return 1.0 / dTimeLogic; }
      static R FpsShell() { return 1.0 / dTimeShell; }

      extern R TimeZoom;

      extern R t;
      extern R dt; // = dTimeLogic * TimeZoom

      void Init();
      void UpdateShell();
      void Destroy();

      void InitLogic();
      void UpdateLogic();
      void FlushLogic();

      const R KmPerAU = 149597870.691; // km
      const R G = 6.67300e-20; // km3*kg-1*s-2
      const R k = 1.3806503e-17; // km2*kg*s-2*K-1

      const R m_sun = 1.98892e30; // kg

      const R m_mercury = 3.3022e23; // kg
      const R d_mercury = 57910000; // avg, km
      const R v_mercury = 47.8725; // avg, km/s

      const R m_venus = 4.8685e24; // kg
      const R d_venus = 108208930; // avg, km
      const R v_venus = 35.0214; // avg, km/s

      const R m_earth = 5.9742e24; // kg
      const R d_earth = 149597870.7; // avg, km
      const R v_earth = 29.780; // avg, km/s

      const R m_moon = 7.347673e22; // kg
      const R d_moon = 384400; // avg, km
      const R v_moon = 1.022; // avg, km/s

      const R m_mars = 6.4185e23; // kg
      const R d_mars = 227936640; // avg, km
      const R v_mars = 24.1309; // avg, km/s

      class CObject;

      extern vector<SmartPtr<CObject>> AllObjects;
      typedef vector<SmartPtr<CObject>>::iterator IterAllObjects;

      class CObject : public virtual IRefCount
      {
      public:
            R3 s, v;
            R m;

            CObject(R3 _s, R3 _v, R _m) { s = _s; v = _v; m = _m; }

            R3 f_g(const CObject *o) const
            {
                     TAKAKU_ASSERT(o != this);

                     R3 Δs = o->s - s;
                     R norm_Δs = Norm(Δs);
                     return G * o->m * Δs / (norm_Δs * norm_Δs * norm_Δs);
            }

            R3 g_total() const
            {
                     R3 ret;
                     for(int i=0; i<AllObjects.size(); i++)
                              if (this != AllObjects[i])
                                    ret += f_g(AllObjects[i]);
                     return ret;
            }

            virtual R3 f_F_others(const CObject *o) const { return R3(); }

            R3 F_others_total() const
            {
                     R3 ret;
                     for(int i=0; i<AllObjects.size(); i++)
                              if (this != AllObjects[i])
                                    ret += f_F_others(AllObjects[i]);
                     return ret;
            }

            R3 a_total() const { return g_total() + F_others_total() / m;}

            R3 v_next() const { return v + (a_total()) * dt; }
            R3 s_next() const { return s + v * dt; }

            R3 _s_next;
            R3 _v_next;

            virtual void Update()
            {
                     //_s_next = s_next();
                     //_v_next = v_next();
                     s = s_next();
                     v = v_next();
            }

            //virtual void Flush()
            //{
            //       s = _s_next;
            //       v = _v_next;
            //}
      };

      class CElecObject : virtual public CObject
      {
      public:
            R e;

            CElecObject(R3 _s, R3 _v, R _m, R _e)
                     : CObject(_s, _v, _m)
            {
                     e = _e;
            }

            virtual R3 f_F_others(const CObject *o) const
            {
                     TAKAKU_ASSERT(o != this);

                     SmartPtr<const CElecObject> _o = dynamic_cast<const CElecObject *>(o);

                     if (_o == null) return R3();

                     R3 D_s = _o->s - s;
                     R norm_D_s = Norm(D_s);
                     return - k * _o->e * e * D_s / (norm_D_s * norm_D_s * norm_D_s);
            }
      };

      class CObjectX : virtual public CObject
      {
      public:
            wstring name;
            int radius;
            dword color;

            CObjectX(wstring _name, R3 _s, R3 _v, R _m, int _radius, dword _color)
                     : CObject(_s, _v, _m)
            {
                     name = _name;
                     radius = _radius;
                     color = _color;
            }
      };

      class CElecObjectX : virtual public CElecObject, virtual public CObjectX
      {
      public:
            CElecObjectX(wstring _name, R3 _s, R3 _v, R _m, R _e, int _radius, dword _color)
                     : CElecObject(_s, _v, _m, _e), CObjectX(_name, _s, _v, _m, _radius, _color), CObject(_s, _v, _m)
            { }
      };

      inline SmartPtr<CObject> GetObject(int i) { return (i >=0 && i < AllObjects.size()) ? AllObjects[i] : null; }

      extern R FovDeg;
      inline R Fov() { return FovDeg / 180; }

      extern int Stands, Looks;
      inline R3 ParentCamPos() { return GetObject(Stands) ? GetObject(Stands)->s : R3(0.0, 0.0, - KmPerAU * 5.0); }
      inline R3 ParentCamTarget() { return GetObject(Looks) ? GetObject(Looks)->s : R3(); }
      inline R3 ParentCamK() { return Normalize(ParentCamTarget() - ParentCamPos()); }
      inline R3 ParentCamI() { return Normalize(Cross(GetObject(Stands) ? R3(0,0,1) : R3(0,1,0), ParentCamK())); }
      inline R3 ParentCamJ() { return Cross(ParentCamK(), ParentCamI()); }
      inline R33 MatParentCamBasis() { return ParentCamI() || ParentCamJ() || ParentCamK(); }

      //const R3 ChildCamPos0(0.0, 0.0, 0.0), ChildCamFront0(0.0, 0.0, KmPerAU), ChildCamUp0(0.0, KmPerAU, 0.0);
      const R3 ChildCamPos0(0.0, -KmPerAU/30, -KmPerAU/10), ChildCamFront0(0.0, KmPerAU/5, KmPerAU), ChildCamUp0(0.0, KmPerAU, 0.0);
      extern R3 ChildCamPos, ChildCamFront, ChildCamUp;

      inline R3 CamPos() { return ParentCamPos() + MatParentCamBasis() * ChildCamPos; }
      inline R3 CamFront() { return MatParentCamBasis() * ChildCamFront; }
      inline R3 CamUp() { return MatParentCamBasis() * ChildCamUp; }

      using namespace System;
      using namespace System::ComponentModel;
      using namespace System::Collections;
      using namespace System::Collections::Generic;
      using namespace System::Windows::Forms;
      using namespace System::Data;
      using namespace System::Drawing;
      using namespace Microsoft::DirectX;
      using namespace Microsoft::DirectX::Direct3D;

      typedef Microsoft::DirectX::Matrix DXMatrix;

      template <class T> inline Vector3 ToDXVec(const T& x) { return Vector3(x(0), x(1), x(2)); }
      template <class T> inline T FromDXVec(Vector3 x) { return T(x.X, x.Y, x.Z); }

      inline String^ ToString(const R3& x) { return L"(" + x(0).ToString() + L", " + x(1).ToString() + L", " + x(2).ToString() + L")"; }

      ref class Form1;

      public ref class G
      {
      public:
            static Form1 ^Form1;

            static Device ^Dev;
            static Texture ^RT;
            static VertexBuffer ^VBRT;
            static VertexBuffer ^VBObject;
            static IndexBuffer ^IBObject;
            static VertexBuffer ^VBLocus;
            static VertexBuffer ^VBAxises;
      };

      value class VertRT
      {
      public:
            Vector4 P;
            Vector2 T;

            VertRT(Vector2 p, float z, Vector2 t)
            {
                     P = Vector4(p.X, p.Y, z, 1.0f);
                     T = t;
            }
      };

      value class VertObject
      {
      public:
            Vector3 P;

            VertObject(Vector3 p)
            {
                     P = p;
            }
      };

      value class VertLocus
      {
      public:
            Vector3 P;
            float Size;
            UInt32 Color;

            VertLocus(Vector3 p, float size, UInt32 color)
            {
                     P = p;
                     Size = size;
                     Color = color;
            }
      };

      value class VertAxis
      {
      public:
            Vector3 P;
            UInt32 Color;

            VertAxis(Vector3 p, UInt32 color)
            {
                     P = p;
                     Color = color;
            }
      };

}

账号		自动登录	找回密码
密码			立即注册

太阳系(重力)仿真

Re:太阳系(重力)仿真

Re: Re:太阳系(重力)仿真

Re:太阳系(重力)仿真

Re: Re:太阳系(重力)仿真

Re: Re: Re:太阳系(重力)仿真

Re: Re: Re: Re:太阳系(重力)仿真

Re:太阳系(重力)仿真