行星运动轨迹的程序实现

这里将以万有引力和势能动能守恒定律为基础,实现行星运动轨迹.然后再假设有两个固定的恒星,让行星在这两个恒星的力场中运动(这是三体问题的一种噢).前面我写过关于混沌曲线的文章:混沌数学及其软件模拟.这类混沌曲线的本质是一个导数方程,即我不知道这条曲线是什么样子,也不知道这条曲线最终通往何处去,我只知道,曲线上的任意一点的切线方向,从而得到它下一点的位置.从而得到整条曲线上的顶点位置.学过物理的人都知道,太阳系中行星的运动轨迹大致是一个椭圆,我们可以知道每个行星的椭圆方程.记得我刚学图形学时,就看过一个程序是太阳行星运动.但它不是以万有引力为基础的,只是让球体绕着固定的椭圆轨迹旋转.

本来我只打算使用万有引力定律,没有考虑势能动能守恒定律.但程序运行后发现,由于没有使用微积分,即使采样间隔时间再小,误差也很大.其实目前的算法误差也不小,呵呵,模拟一下吧,不要太计较.

先帖下基类代码,这代码与混沌数学及其软件模拟中的很相似,即是一个导数方程,用于由当前点计算下一点.

class DifferentiationFunction
{
public:
    virtual void Differentiate(float x, float y, float z, float t,
        float& outX, float& outY, float& outZ) = NULL;

    virtual float GetExtendT() const = NULL;

    virtual float GetStartX() const
    {
        return 0.0f;
    }

    virtual float GetStartY() const
    {
        return 0.0f;
    }

    virtual float GetStartZ() const
    {
        return 0.0f;
    }
};

行星方程:

// 行星的导数曲线
class Planet : public DifferentiationFunction
{
public:
    Planet()
    {
        m_star_weight = 1.0f;

        m_planet_x = 5.0f;
        m_planet_y = 8.0f;
        m_planet_z = 1.0f;
        m_planet_weight = 0.1f;
        m_planet_speed_x = 4.0f;
        m_planet_speed_y = 0.0f;
        m_planet_speed_z = 0.0f;

        m_g = 100.0f;

        m_ek = 0.5f*m_planet_weight*(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z); // 1/2*m*v*v
        float r = sqrt(m_planet_x*m_planet_x + m_planet_y*m_planet_y + m_planet_z*m_planet_z);
        m_ep = -/*0.5f**/m_g*m_star_weight*m_planet_weight/r;
        m_e = m_ek + m_ep;
    }

    void Differentiate(float x, float y, float z, float t,
        float& outX, float& outY, float& outZ)
    {
        t = t*10.0f;

        float sqd = x*x + y*y + z*z;
        float d = sqrt(sqd);

        float a = m_g*m_star_weight/sqd;
        float ax = -a*x/d;
        float ay = -a*y/d;
        float az = -a*z/d;

        outX = x + m_planet_speed_x*t + 0.5f*ax*t*t;
        outY = y + m_planet_speed_y*t + 0.5f*ay*t*t;
        outZ = z + m_planet_speed_z*t + 0.5f*az*t*t;

        m_planet_speed_x += ax*t;
        m_planet_speed_y += ay*t;
        m_planet_speed_z += az*t;

        float r = sqrt(outX*outX + outY*outY + outZ*outZ);
        m_ep = -/*0.5f**/m_g*m_star_weight*m_planet_weight/r;
        m_ek = m_e - m_ep;
        if (m_ek < 0.0f)
        {
            m_ek = 0.0f;
        }
        float v = sqrt(2*m_ek/m_planet_weight);

        float w = sqrt(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z);

        m_planet_speed_x *= v/w;
        m_planet_speed_y *= v/w;
        m_planet_speed_z *= v/w;
    }

    float GetExtendT() const
    {
        return 20.0f;
    }

    float GetStartX() const
    {
        return m_planet_x;
    }

    float GetStartY() const
    {
        return m_planet_y;
    }

    float GetStartZ() const
    {
        return m_planet_z;
    }

public:
    float m_star_weight;

    float m_planet_x;
    float m_planet_y;
    float m_planet_z;
    float m_planet_weight;
    float m_planet_speed_x;
    float m_planet_speed_y;
    float m_planet_speed_z;

    float m_g;      // 万有引力系数

    float m_e;
    float m_ek;     // 动能
    float m_ep;     // 引力势能
};

行星在这两个恒星的力场中运动(三体问题)

// 三体问题的导数曲线
class ThreeBody : public DifferentiationFunction
{
public:
    ThreeBody()
    {
        m_star1_x = -10.0f;
        m_star1_y = 0.0f;
        m_star1_z = 0.0f;
        m_star1_weight = 1.0f;

        m_star2_x = 10.0f;
        m_star2_y = 0.0f;
        m_star2_z = 0.0f;
        m_star2_weight = 1.0f;

        m_planet_x = 5.0f;
        m_planet_y = 5.0f;
        m_planet_z = 0.1f;
        m_planet_weight = 0.1f;
        m_planet_speed_x = 0.0f;
        m_planet_speed_y = 2.0f;
        m_planet_speed_z = 0.0f;

        m_g = 50.0f;

        m_ek = 0.5f*m_planet_weight*(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z); // 1/2*m*v*v

        float d1x = m_star1_x - m_planet_x;
        float d1y = m_star1_y - m_planet_y;
        float d1z = m_star1_z - m_planet_z;
        float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z;
        float d1 = sqrt(sqd1);
        m_ep1 = -m_g*m_star1_weight*m_planet_weight/d1;

        float d2x = m_star2_x - m_planet_x;
        float d2y = m_star2_y - m_planet_y;
        float d2z = m_star2_z - m_planet_z;
        float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z;
        float d2 = sqrt(sqd2);
        m_ep2 = -m_g*m_star2_weight*m_planet_weight/d2;

        m_e = m_ek + m_ep1 + m_ep2;
    }

    void Differentiate(float x, float y, float z, float t,
        float& outX, float& outY, float& outZ)
    {
        t = t*20.0f;

        float d1x = m_star1_x - x;
        float d1y = m_star1_y - y;
        float d1z = m_star1_z - z;
        float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z;
        float d1 = sqrt(sqd1);

        float d2x = m_star2_x - x;
        float d2y = m_star2_y - y;
        float d2z = m_star2_z - z;
        float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z;
        float d2 = sqrt(sqd2);

        float a1 = m_g*m_star1_weight/sqd1;
        float a1x = a1*d1x/d1;
        float a1y = a1*d1y/d1;
        float a1z = a1*d1z/d1;

        float a2 = m_g*m_star2_weight/sqd2;
        float a2x = a2*d2x/d2;
        float a2y = a2*d2y/d2;
        float a2z = a2*d2z/d2;

        outX = x + m_planet_speed_x*t + 0.5f*(a1x + a2x)*t*t;
        outY = y + m_planet_speed_y*t + 0.5f*(a1y + a2y)*t*t;
        outZ = z + m_planet_speed_z*t + 0.5f*(a1z + a2z)*t*t;

        m_planet_speed_x += (a1x + a2x)*t;
        m_planet_speed_y += (a1y + a2y)*t;
        m_planet_speed_z += (a1z + a2z)*t;

        {
            float d1x = m_star1_x - outX;
            float d1y = m_star1_y - outY;
            float d1z = m_star1_z - outZ;
            float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z;
            float d1 = sqrt(sqd1);
            m_ep1 = -m_g*m_star1_weight*m_planet_weight/d1;

            float d2x = m_star2_x - outX;
            float d2y = m_star2_y - outY;
            float d2z = m_star2_z - outZ;
            float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z;
            float d2 = sqrt(sqd2);
            m_ep2 = -m_g*m_star2_weight*m_planet_weight/d2;

            m_ek = m_e - m_ep1 - m_ep2;
            if (m_ek < 0.0f)
            {
                m_ek = 0.0f;
            }
            float v = sqrt(2*m_ek/m_planet_weight);

            float w = sqrt(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z);

            m_planet_speed_x *= v/w;
            m_planet_speed_y *= v/w;
            m_planet_speed_z *= v/w;
        }
    }

    float GetExtendT() const
    {
        return 20.0f;
    }

    float GetStartX() const
    {
        return m_planet_x;
    }

    float GetStartY() const
    {
        return m_planet_y;
    }

    float GetStartZ() const
    {
        return m_planet_z;
    }

public:
    float m_star1_x;
    float m_star1_y;
    float m_star1_z;
    float m_star1_weight;

    float m_star2_x;
    float m_star2_y;
    float m_star2_z;
    float m_star2_weight;

    float m_planet_x;
    float m_planet_y;
    float m_planet_z;
    float m_planet_weight;
    float m_planet_speed_x;
    float m_planet_speed_y;
    float m_planet_speed_z;

    float m_g;      // 万有引力系数

    float m_e;
    float m_ek;     // 动能
    float m_ep1;    // 引力势能
    float m_ep2;    // 引力势能
};

这是我写的测试程序,写它是为下一步的三体模拟软件做准备.

时间: 09-17

行星运动轨迹的程序实现的相关文章

三体运动的程序模拟

前几天看了<三体>,很不错的科幻小说.说到三体,我想到我大学的一个舍友叫王晶,和香港那个导演同名同姓同性别.记得有一次几个同学在一块聊天,有个女生问他:父母为什么给他取名叫晶.他说叫晶是父母希望能有三个太阳守护着他.那时我还很单纯,不会用五行缺什么的话来讽刺他,只是说,如果给他起名叫王晶晶的话,那就有6个太阳守护他了.现在对三体有了一些了解,才意识到被三个太阳罩着,那不叫守护,应该是被蹂躏.三体内的行星,感觉是被三个恒星玩弄于股掌之间,如同球一样踢来踢去. 空间中三个星体,受万有引力作用下的运

混沌图像---三体的纠结

最简单而引人注目的混沌,莫过于三体运动.仅仅三颗星体的运动,就能变得复杂而眩目.这种复杂曾令数学家们在百年间困惑不已.如果只有两个天体,那么一切是多么简单,18世纪的伯努利就已解出了运动的所有可能轨迹,用合适的坐标,就能用简单的曲线描述.但仅仅是多了一个天体,就要等到19世纪的庞加莱,才给出了差强人意的答案:没有漂亮的解(正式术语是三体系统是不可积的).这并非因为人类的智慧所限,而是从本质上来说,三个天体之间的运动轨迹不可能用简单的式子表达.自然并不像原来期盼的那么简单,它的复杂性令人绝望.小说

fsd

三体运动的程序模拟 前几天看了<三体>,很不错的科幻小说.说到三体,我想到我大学的一个舍友叫王晶,和香港那个导演同名同姓同性别.记得有一次几个同学在一块聊天,有个女生问他:父母为什么给他取名叫晶.他说叫晶是父母希望能有三个太阳守护着他.那时我还很单纯,不会用五行缺什么的话来讽刺他,只是说,如果给他起名叫王晶晶的话,那就有6个太阳守护他了.现在对三体有了一些了解,才意识到被三个太阳罩着,那不叫守护,应该是被蹂躏.三体内的行星,感觉是被三个恒星玩弄于股掌之间,如同球一样踢来踢去. 空间中三个星体,

【CSS3】变换transform---小案例,行星运动效果

transform变换 rotate旋转 rotate(angle) 2D旋转:rotateX(angle) 沿着x轴旋转:rotateY(angle) 沿着y轴旋转:rotate(angle) 沿着z轴旋转: deg角度.rad弧度.turn圈.grad梯度 90deg = 100grad = 0.25turn ≈  1.570796326794897rad scale缩放 0<取值<1 缩小:取值>1 放大. translate 位移 translate(0px,0px);trans

第二章 希腊化——罗马科学

第二章 希腊化--罗马科学2.1希腊数理天文学的起源希腊天文学的基本动机来自于拯救现象(将纷乱的现象规整)欧多克索斯--同心球模型,一个行星做多个匀速圆周运动,叠加成为复杂的运动,同心球模型不仅解决了柏拉图的拯救现象,同时奠定了希腊天文学的基本模式2.1希腊数理天文学的起源(习题)2.2希腊数理天文学的代表阿基米德(约前287-前212年)反射杠杆原理阿基米德重要的数学成就阿基米德的物理成就阿基米德的科学,把抽象的.演绎的科学,与实用联系起来托勒密托勒密的希腊数理天文学的基本模型托勒密体系代表了

巧妙利用CSS3实现太阳系行星公转运动轨迹

前段时间在博客园看到一篇很不错的文章,就是介绍如何用CSS3来实现太阳系中行星的运动.本小姐处于好奇心从头到完整的看完了.突然发现好奇不仅仅害死猫,也可以增长知识嘛(你这是什么谬论!!!).也许有的朋友看过有关的文章那么今天我们就一起来看看怎么巧妙的利用CSS3 来实现这么漂亮的动画的吧! 我不清楚大家对于CSS3的动画了解多少,但是我对于 CSS3的东西还是限制与用什么查什么的阶段.css3的animation动画,结合transform的大小.旋转.位移.斜切,通过两三行代码,便可做出很多有

谬论之程序员的眼光看世界

曾经阅读霍金先生的时间简史,书中说世界的第四维是时间.当时感觉高深莫测,现在作为程序员,我来说说以我程序员的眼光是怎么来看宇宙的吧. 如今的程序员大部分在使用面向对象的语言来编程.所以,万物皆对象对于程序员来说,还是比较容易理解的.我重新衡量世界观也是从此开始的. 三维世界中,给定原点后,所有的物体都可以使用长宽高来表示.即如今我们的世界只有长宽高三个维度.用数组来表示,可以看做一个三维数组.即用三维数组中的各个点,就可以表示出所有物体的形状. 大家都知道,物质的最小单位是质子,电子,中子等  

程序员职业规划:让自己变得重要

昨天在Better Software Magazine上看到一篇文章Make Yourself Essential(杂志不在手边,文章名和杂志期号都记不清了),谈到了IT市场变幻,程序员风光不再,身为程序员应该怎么应对.作者开篇明义:工作外包,公司裁员,是公司的错么?当然不是!一切都是程序员的错.当我们技术沦为大宗商品(所谓大宗商品,就是说商品和商品间的区别主要就是价格了)的时候,我们被收费更低的竞争对手代替的日子也就不远了.而随着技术的进步,很多编程技术也就不可避免地成为大宗商品.而要想增加自

Java_太阳系_行星模型_小游戏练习_详细注释

1 //实现MyFrame--实现绘制窗口,和实现重写 重画窗口线程类 2 3 package cn.xiaocangtian.Test; 4 5 import java.awt.Frame; 6 import java.awt.event.WindowAdapter; 7 import java.awt.event.WindowEvent; 8 9 10 public class MyFrame extends Frame { 11 12 //加载窗口 13 public void launc