<正经代码> 三体0.5

· · 个人记录

思维导读:

三体问题 (天体力学中的基本力学模型)

三体问题是天体力学中的基本力学模型。

它是指三个质量、初始位置和初始速度都是任意的可视为质点的天体,在相互之间万有引力的作用下的运动规律问题。 现在已知,三体问题不能精确求解,即无法预测所有三体问题的数学情景,只有几种特殊情况已研究。

三体问题(three-body problem)

天体力学中的基本力学模型。研究三个可视为质点的天体在相互之间万有引力作用下的 运动规律 问题。这三个天体的质量、初始位置和初始速度都是任意的。最简单的一个例子就是太阳系中太阳、地球和月球的运动。在浩瀚的宇宙中,星球的大小可以忽略不计,所以我们可以把它们看成质点。如果不计太阳系其他星球的影响,那么它们的运动就只是在引力的作用下产生的,所以我们就可以把它们的运动看成一个三体问题。

研究方法

由于三体问题不能严格求解,在研究天体运动时,都只能根据实际情况采用各种近似的解法,研究三体问题的方法大致可分为3类:

第一类是分析方法,其基本原理是把天体的坐标和速度展开为时间或其他小参数的级数形式的近似分析表达式,从而讨论天体的坐标或轨道要素随时间的变化;

第二类是定性方法,采用微分方程的定性理论来研究长时间内三体运动的宏观规律和全局性质;

第三类是数值方法,这是直接根据微分方程的计算方法得出天体在某些时刻的具体位置和速度。

这三类方法各有利弊,对新积分的探索和各类方法的改进是研究三体问题中很重要的课题。

对此,我想说。。。

还有暴力模拟法啊!

程序解释:

自己改速度与放大倍数(出界会速度清零),白的是三体行星,还有三个太阳。

#include<bits/stdc++.h>
#include<windows.h>
#include<conio.h>
using namespace std;
int toint(float a){return ((int)(a*10+5))/10;}
void Setpos(float x,float y){COORD pos;pos.X=toint(y*2),pos.Y=toint(x);SetConsoleCursorPosition(GetStdHandle(STD_OUTPUT_HANDLE),pos);}
void Color(int a)
{
    if(a==4) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED|FOREGROUND_GREEN|FOREGROUND_BLUE);
    if(a==1) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_GREEN|FOREGROUND_BLUE);
    if(a==2) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_GREEN);
    if(a==3) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED);
//    if(a==4) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED|FOREGROUND_BLUE);
//    if(a==5) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED|FOREGROUND_GREEN);
//    if(a==6) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_BLUE);
//    if(a==7) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED|FOREGROUND_GREEN|FOREGROUND_BLUE);
//    if(a==7) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_GREEN|FOREGROUND_BLUE);
//    if(a==8) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY);
//    if(a==9) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_GREEN);
//    if(a==10) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_RED);
//    if(a==11) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_RED|FOREGROUND_GREEN|FOREGROUND_BLUE);
}
struct node
{
    float x,y,z;
    float vx,vy,vz;
    float r,m;
    bool life;
}Sun[1000001];
void Push(int a,int b)
{
    float Ax=Sun[a].x-Sun[b].x,Ay=Sun[a].y-Sun[b].y,Dis=sqrt(Ax*Ax+Ay*Ay)*1.0;
    if(Dis<=1) return;//未撞: Dis>=Sun[a].r||Dis>=Sun[b].r||

    float ac=Sun[a].m*Sun[b].m/(Dis*Dis)*1.0,afx,afy;
    if(abs(Ax)<=1) return;//afy=ac;
    else if(abs(Ay)<=1) return;//afx=ac;
    else
    {
        float d=abs(Ax/Ay*1.0);
        afy=sqrt(ac/(1+d*d))*1.0,afx=sqrt(ac/(1+d*d))*d*1.0;
    }
    if(Ax>0) afx*=-1;if(Ay>0) afy*=-1;
    Sun[a].vx+=afx/Sun[a].m*5.0;Sun[a].vy+=afy/Sun[a].m*5.0;
    Sun[b].vx-=afx/Sun[b].m*5.0;Sun[b].vy-=afy/Sun[b].m*5.0;
}
int b,T,More=10,Speed=100;
void Move()
{
    for(int i=1;i<=b;i++)
    {
        Setpos(Sun[i].x/More*1.0,Sun[i].y/More*1.0);cout<<"  ";
        Sun[i].x+=Sun[i].vx*More/10.0;Sun[i].y+=Sun[i].vy*More/10.0;
        if(Sun[i].x>=40*More) Sun[i].x=40*More,Sun[i].vx=0;
        if(Sun[i].x<=1) Sun[i].x=1,Sun[i].vx=-0;
        if(Sun[i].y>=40*More) Sun[i].y=40*More,Sun[i].vy=0;
        if(Sun[i].y<=1) Sun[i].y=1,Sun[i].vy=0;
        Setpos(Sun[i].x/More*1.0,Sun[i].y/More*1.0);Color(i);cout<<"●";
    }
    for(int i=1;i<=b;i++)
    for(int j=i+1;j<=b;j++) Push(i,j);
}
int main()
{
    system("mode con cols=82 lines=42");
    CONSOLE_CURSOR_INFO cursor_info={1,0};
    SetConsoleCursorInfo(GetStdHandle(STD_OUTPUT_HANDLE),&cursor_info);
    srand((unsigned)time(NULL));

//  自己改速度与放大倍数(出界会速度清零),白的是三体行星,还有三个太阳。 

    b=4;
    Sun[1].x=20,Sun[1].y=20,Sun[1].vx=5,Sun[1].vy=10,Sun[1].m=1000;
    Sun[2].x=5,Sun[2].y=5,Sun[2].vx=1,Sun[2].vy=-5,Sun[2].m=1000;
    Sun[3].x=10,Sun[3].y=15,Sun[3].vx=10,Sun[3].vy=-9,Sun[3].m=1000;
    Sun[4].x=13,Sun[4].y=30,Sun[4].vx=1,Sun[4].vy=1,Sun[4].m=1;//改这里! 
    while(1)
    {
        T++;
        Move();
        Sleep(Speed);
    }
}

另:娱乐代码:

#include<bits/stdc++.h>
#include<windows.h>
#include<conio.h>
using namespace std;
int toint(float a){return ((int)(a*10+5))/10;}
void Setpos(float x,float y){COORD pos;pos.X=toint(y*2),pos.Y=toint(x);SetConsoleCursorPosition(GetStdHandle(STD_OUTPUT_HANDLE),pos);}
void Color(int a)
{
    if(a%12==0) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED|FOREGROUND_GREEN|FOREGROUND_BLUE);
    if(a%12==1) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_GREEN|FOREGROUND_BLUE);
    if(a%12==2) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED|FOREGROUND_GREEN);
    if(a%12==3) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED);
    if(a%12==4) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED|FOREGROUND_BLUE);
    if(a%12==5) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_GREEN);
    if(a%12==6) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_BLUE);
    if(a%12==7) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_INTENSITY|FOREGROUND_RED|FOREGROUND_GREEN|FOREGROUND_BLUE);
    if(a%12==8) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_GREEN);
    if(a%12==9) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_RED);
    if(a%12==10) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_RED|FOREGROUND_GREEN|FOREGROUND_BLUE);
    if(a%12==11) SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),FOREGROUND_GREEN|FOREGROUND_BLUE);
}
struct node
{
    float x,y,z;
    float vx,vy,vz;
    float r,m;
    bool life;
}Sun[1000001];
void Push(int a,int b)
{
    float Ax=Sun[a].x-Sun[b].x,Ay=Sun[a].y-Sun[b].y,Dis=sqrt(Ax*Ax+Ay*Ay)*1.0;
    if(Dis<=20) return;//未撞: Dis>=Sun[a].r||Dis>=Sun[b].r||

    float ac=Sun[a].m*Sun[b].m/(Dis*Dis)*1.0,afx=0,afy=0;
    if(abs(Ay)<=1) Ay=1;
    float d=abs(Ax/Ay*1.0);
    afy=sqrt(ac/(1+d*d))*1.0,afx=sqrt(ac/(1+d*d))*d*1.0;
    if(Ax>0) afx*=-1;if(Ay>0) afy*=-1;
    if(abs(Ay)>1) Sun[a].vx+=afx/Sun[a].m*10.0;
    Sun[a].vy+=afy/Sun[a].m*10.0;
    if(abs(Ay)>1) Sun[b].vx-=afx/Sun[b].m*10.0;
    Sun[b].vy-=afy/Sun[b].m*10.0;
}
int b,T,More=10,Speed=15;
void Move()
{
    for(int i=1;i<=b;i++)
    {
        Setpos(Sun[i].x/More*1.0,Sun[i].y/More*1.0);cout<<"  ";
        Sun[i].x+=Sun[i].vx*More/10.0;Sun[i].y+=Sun[i].vy*More/10.0;
        if(Sun[i].x>=40*More) Sun[i].x=40*More,Sun[i].vx=0;
        if(Sun[i].x<=1) Sun[i].x=1,Sun[i].vx=-0;
        if(Sun[i].y>=40*More) Sun[i].y=40*More,Sun[i].vy=0;
        if(Sun[i].y<=1) Sun[i].y=1,Sun[i].vy=0;
        Setpos(Sun[i].x/More*1.0,Sun[i].y/More*1.0);Color(i);cout<<"●";
    }
    for(int i=1;i<=b;i++)
    for(int j=i+1;j<=b;j++) Push(i,j);
}
int rand(int a){return rand()%a;}
void Start(int a) {b++;Sun[a].x=rand(40)*More,Sun[a].y=rand(40)*More,Sun[a].vx=(rand(41)-20)*More/100.0,Sun[a].vy=(rand(41)-20)*More/100.0,Sun[a].m=1000;}//(rand(20)+1)*100;}
int main()
{
    system("mode con cols=82 lines=42");
    CONSOLE_CURSOR_INFO cursor_info={1,0};
    SetConsoleCursorInfo(GetStdHandle(STD_OUTPUT_HANDLE),&cursor_info);
    srand((unsigned)time(NULL));

//  自己改速度与放大倍数(出界会速度清零),白的是三体行星,还有三个太阳。 

    for(int i=1;i<=100;i++) Start(i);
    while(1)
    {
        Setpos(3,5);Color(0);cout<<"三 体 模 型 0.5     (  n*太阳 (质点)  )"; 
        T++;
        Move();
        Sleep(Speed);
    }
}