python引力模型网络分析引力模型建模

转载

数据探索家 2024-08-22 21:35:52

文章标签 python引力模型网络分析 c++ OpenGL 弹性碰撞 2d 文章分类 Python 后端开发

OpenGL（C++，freeglut库）实现多星引力运动模型

简介
物理基础

万有引力
三维空间的球体完全非弹性碰撞

模型的代码描述

向量
小球
运动系统

交互的代码描述

视角控制
模型选择

渲染过程
程序逻辑结构图

简介

学校的专业选修课让我接触到了计算机图形学这门课，让我终于有机会从冷冰冰的控制台应用程序中解放一回。于是我用OpenGL做了自己的第一个小项目，由于水平原因，效果可能没有想象中的惊艳，程序的重点部分都是C++代码。这里我们建立一个多星运动模型，利用OpenGL进行光照渲染，引力模型共有1~7星可供选择，待我接下来慢慢展开讲述。

物理基础

万有引力

牛顿提出的万有引力定律:将两个物体视为质点，两者的质量分别为 $python引力模型网络分析引力模型建模_弹性碰撞$ , $python引力模型网络分析引力模型建模_python引力模型网络分析_02$ ，两者的距离为 $python引力模型网络分析引力模型建模_OpenGL_03$ ，那么两者受到对方施加的引力为：
$python引力模型网络分析引力模型建模_python引力模型网络分析_04$
其中 $python引力模型网络分析引力模型建模_2d_05$ ， $python引力模型网络分析引力模型建模_c++_06$ 是受力方指向施力方的单位向量。
也可以表示为以下形式：
$python引力模型网络分析引力模型建模_OpenGL_07$
其中 $python引力模型网络分析引力模型建模_python引力模型网络分析_08$ 是以受力方质点位置为起点，指向施力方质点位置为终点的向量， $python引力模型网络分析引力模型建模_python引力模型网络分析_09$ 是其模。

三维空间的球体完全非弹性碰撞

假设两个小球碰撞瞬间参数如下图：

python引力模型网络分析引力模型建模_弹性碰撞_10

其中 $python引力模型网络分析引力模型建模_OpenGL_11$ 是小球1撞前速度， $python引力模型网络分析引力模型建模_c++_12$ 是小球1的撞后速度，小球2同理。

求解模型：

由动量守恒知：

$python引力模型网络分析引力模型建模_弹性碰撞_13$ 完全非弹性碰撞，能量守恒（只有动能，向量形式）：

$python引力模型网络分析引力模型建模_c++_14$ 将

$python引力模型网络分析引力模型建模_2d_15$ $python引力模型网络分析引力模型建模_2d_16$ 代入（2）式有：

$python引力模型网络分析引力模型建模_python引力模型网络分析_17$ 小球在碰撞的瞬间受力分析，可知两个球所受力的方向一定与 $python引力模型网络分析引力模型建模_python引力模型网络分析_08$ 共线，于是设

$python引力模型网络分析引力模型建模_python引力模型网络分析_19$ 结合（1）式代入（3）式中可解的：

$python引力模型网络分析引力模型建模_弹性碰撞_20$ 可得：

$python引力模型网络分析引力模型建模_OpenGL_21$ $python引力模型网络分析引力模型建模_弹性碰撞_22$ 求解完成

模型的代码描述

基于上面介绍的两个物理模型，我们在物理上已经得到了我们模型的所有约束条件，接下来就可以用代码来描述我们的模型了。vs2019中freegult库的安装就不介绍了。

向量

在我们的模型中，以及OpenGL渲染过程中描述对象的位置都是用的向量，这里我们就实现自己的一个向量类型。

class Vector {
private:
	float x, y, z;//三维空间向量坐标
	friend class RMatrix;
public:
	Vector();
	Vector(float _x, float _y, float _z);
	Vector(const Vector& V);//三个构造函数
	float X() const;
	float Y() const;
	float Z() const;//分别返回三个坐标的值

	void clear();//表示力的时候，每次计算新的力都得先将力变为0

	Vector operator^ (const Vector& V); //叉乘，定义成const貌似会更好

	float operator * (const Vector& V) const;

	Vector operator-(const Vector& V) const;

	Vector operator + (const Vector& V) const;
	Vector& operator += (const Vector& V);

	Vector operator * (float power) const;
	Vector& operator *= (float power);

	Status Rotate(const Vector& n , float angle);//旋转绕轴旋转

	Status Translate(float dx, float dy, float dz);
	Status Normalize();//单位化
	float Angle(Vector _r) const;

	float R() const;

	friend std::ostream& operator<< (std::ostream& out, const Vector& V);
};

class RMatrix {
private:
	int row, col;
	float base[9];
	friend class Vector;
public:
	RMatrix();
	RMatrix(char axis, float sida);//构造绕X或Y或Z轴的旋转矩阵
//	RMatrix(const RMatrix& obj);
	RMatrix operator* (const RMatrix& other);
	RMatrix& operator*= (const RMatrix& other);
	Vector operator* (const Vector& v);
};

上面我们声明了一个向量类Vector以及一个辅助向量旋转的旋转矩阵类RMatix。Vector中包几乎含了描述空间位置以及空间运动的所有方法：内积，外积，数乘，加法，减法，绕任意轴旋转，作为点的时候的平移，以及求与另一向量的夹角，RMatrix辅助向量旋转，并没有太多的方法。基于这些方法，我们就可以描述一个球在空间中的位置，受力，速度，以及表面材质等，并且可以直接利用定义的方法完成受力计算和碰撞计算。由于这些声明需要的定义都比较简单，就不再这里展示，感兴趣的话详细定义代码见附录链接。¹

小球

小球的声明如下:

class Ball {
private:
	Vector o,//球心
		f,//受力
		v,//速度
		material;//材质
	float m, //质量
		r;//半径
	bool light;//自身是否发光
public :
	Ball();
	Ball(const Vector& _o, const Vector& _f, const Vector& _v , const Vector& _material, float _m, float _r , bool _light); 
	void Draw();//将小球的数据传入OpenGL渲染管线
	float distance(const Ball& b);//两球之间的距离
	Vector force(const Ball& b);//该球受到另一个球的引力
	void translateto(const Vector& _do);//移动
	void clearforce();//清除力
	void addforce(const Vector& _df);//添加力
	void addvelocity(const Vector& _dv);//添加速度
	void nextPosition(Vector& temppoistion, float dt);//下一个位置

	Vector O();
	Vector F();
	Vector V();
	float M();
	float R();
	bool L();
};

这个声明给出了小球类Ball包含的属性和方法，属性：球的球心，受力，速度，材质，质量，半径，是否发光，以及方法：求距离，求受力，移动，速度变化，受力变化，以及计算当前受力的下一个位置。用这些方法，通过循环模拟，我们可以描述没有碰撞的运动模型。这里补充一下模型假设：在极短的一个时间段内，认为小球受到的力不会改变。给出两个函数体的定义：

void Ball::Draw()
{
	if (light == true) {
		GLfloat lightPosition[] = { o.X() , o.Y() , o.Z() , 1 };
		glEnable(GL_LIGHTING);
		glEnable(GL_LIGHT0);

		glPushMatrix();
		glLightfv(GL_LIGHT0, GL_AMBIENT, lightAmbient);
		glLightfv(GL_LIGHT0, GL_DIFFUSE, lightDiffuse);
		glLightfv(GL_LIGHT0, GL_SPECULAR, lightSpecular);
		glLightfv(GL_LIGHT0, GL_POSITION, lightPosition);//这几个数组在其他部分是给出了定义的
		
		glLightf(GL_LIGHT0, GL_CONSTANT_ATTENUATION, 0.5);
		glLightf(GL_LIGHT0, GL_LINEAR_ATTENUATION, 0);
		glLightf(GL_LIGHT0, GL_QUADRATIC_ATTENUATION, 5);
		glPopMatrix();

		glMaterialfv(GL_FRONT, GL_EMISSION, EMISSION);
	}
	Material[0] = material.X();
	Material[1] = material.Y();
	Material[2] = material.Z();//全局的数组

	glPushMatrix();
	glMaterialfv(GL_FRONT, GL_SPECULAR, Material);
	glTranslatef(o.X(), o.Y(), o.Z());
	glutSolidSphere(r, 40, 32);
	glPopMatrix();
}
void Ball::nextPosition(Vector& temppoistion, float dt)
{
	float ax = f.X() / m,
		ay = f.Y() / m,
		az = f.Z() / m;
	temppoistion = o + Vector(v.X() * dt + 0.5 * ax * dt * dt, v.Y() * dt
		+ 0.5 * ay * dt * dt, v.Z() * dt + 0.5 * az * dt * dt);
}//按模型假设来的

其余部分的定义代码感兴趣的小伙伴们就转附录链接吧（其实所有代码都在这里面）。¹

运动系统

上面给出了小球的描述，但是我们最终要描述的是一个运动系统。于是我们又定义了一个系统类System用于将多个小球联系起来。声明如下:

struct System {
	int n;
	Ball* balls;
	System();
	System(int _n, const Ball* _balls);
	System(const System& s);
	void simulate();
	System& operator= (const System& s);

	~System();
};

void recall(Ball& bone, Vector& tpone, Ball& btwo, Vector& tptwo, float low, float high, float top);

其中System属性包含构成系统的小球的数量，以及指向这些小球的指针，方法最主要的是模拟一次短时间（约定好的）系统的变化。后面定义的recall函数是用来解决碰撞问题的一个递归回溯。之所以需要递归回溯，是因为可能就在模拟的这一段时间内，小球在前一部分的假设下，可能早就出现了碰撞并且甚至可能出现了重叠，这显然是不符合物理中非场物质的空间排他性，所以需要搜索碰撞发生的时间点并用碰撞模型改变两个小球的速度方向。下面给出这两个函数的函数体：

void System::simulate()
{
	float dt = 0.05;
	Vector* tempposition = new Vector[n]; //记录新位置
	int i, j;
	for (i = 0; i < n; i++)
		balls[i].clearforce(); //清除受力
	for(i = 0 ; i < n ; i++)
		for (j = i + 1; j < n; j++) {
			Vector tempf = balls[i].force(balls[j]);
			balls[i].addforce(tempf);
			balls[j].addforce(tempf * -1.0);//对每一个小球计算受力
		}
	for (i = 0; i < n; i++)
		balls[i].nextPosition(tempposition[i], dt);//将小球的新位置存在数组中

	for (i = 0; i < n; i++)
		for (j = i + 1 ; j < n; j++) {
			Vector sub = tempposition[j] - tempposition[i];
			if (sub.R() < (balls[i].R() + balls[j].R()) * 0.99)
				recall(balls[i], tempposition[i], balls[j], tempposition[j], 0, dt, dt);
		}//检测是否有小球相撞，如果有递归调用函数recall
	for (i = 0; i < n; i++) {
		float ax = balls[i].F().X() / balls[i].M();
		float ay = balls[i].F().Y() / balls[i].M();
		float az = balls[i].F().Z() / balls[i].M();
		balls[i].addvelocity(Vector(ax * dt, ay * dt, az * dt));
		balls[i].translateto(tempposition[i]);//改变小球位置及速度
	}
	delete[] tempposition;//释放数组
}
void recall(Ball& bone, Vector& tpone, Ball& btwo, Vector& tptwo, float low, float high, float top)
{
	double mid = (low + high) / 2;
	double R;
	bone.nextPosition(tpone, mid);
	btwo.nextPosition(tptwo, mid);
	R = (tpone - tptwo).R();
	if (high - low < accuracy) {//如果搜索时间段已经足够小，则认为小球在此时发生碰撞
		bone.translateto(tpone);
		btwo.translateto(tptwo);
		Vector temp = tptwo - tpone;
		float k = (temp * (bone.V() - btwo.V())) * 2 / (R * R * (bone.M() + btwo.M()));
		btwo.addvelocity(temp * (bone.M() * k));
		bone.addvelocity(temp * (btwo.M() * k * -1.0));
		bone.nextPosition(tpone, top - mid);
		btwo.nextPosition(tptwo, top - mid);
		return;
	}
	if (R > (bone.R() + btwo.R()) / 0.99) {//未碰撞，搜索更晚时间
		recall(bone, tpone, btwo, tptwo, mid, high, top);
		return;
	}
	else if (R < (bone.R() + btwo.R()) * 0.99) {//早已发生碰撞，搜索更早时间
		recall(bone, tpone, btwo, tptwo, low, mid, top);
		return;
	}
	else {//碰撞时，使用碰撞模型
		bone.translateto(tpone);
		btwo.translateto(tptwo);
		Vector temp = tptwo - tpone;
		double k = (temp * (bone.V() - btwo.V())) * 2 / (R * R * (bone.M() + btwo.M()));
		btwo.addvelocity(temp * (bone.M() * k));
		bone.addvelocity(temp * (btwo.M() * k * -1.0));
		bone.nextPosition(tpone, top - mid);
		btwo.nextPosition(tptwo, top - mid);
		return;
	}
}

通过间隔一段时间，调用一次simulate()函数，我们就能模拟整个系统运动的过程了。剩余定义的代码就见附录链接。¹

交互的代码描述

在描述完我们的运动模型之后，接下来就剩下OpenGL渲染了。渲染过程就是freeglut库使用的基本过程：初始化，建立窗口，自定义的初始化，注册各种回调函数，进入大循环。
这部分后面讨论。
除了模型描述，还有很重要的交互功能要实现：选择14个模型中的一个，或者切换另一个模型：输入stable1，。。。stable
7，unstable1,…,unstable7选择或切换模型。视角控制：前进或后退，以及旋转。光线衰减模式选择：输入yes光线衰减，并且在此模式下输入+光想变亮，输入-光线变暗；输入no光线不衰减。

视角控制

我们采用glLookAt()函数来观察我们建立的模型，为了更好的使用这个函数，我们定义了一个类来对应其九个参数，声明如下：

class mycamera {
private :
	Vector o, //位置
		n, //视线朝向
		v; //头朝向
public :
	mycamera();
	void reset();//每次切换模型调用，使得照相机回到最初的状态
	Vector eye();//投影点，对应glLookAt()前三个参数
	Vector center();//视角中心，对应glLookAt()中间三个参数
	Vector up();//上方向，对应glLookAt()后三个参数
	void forward(float step);//前后平移
	void nod(float sida);//上下旋转
	void shade(float sida);//左右旋转
};

通过这个mycamera类配合glLookAt()函数，我们可以实现视角的平移与旋转（缺少了左右平移，有兴趣的可以小伙伴可以自己加上）。

模型选择

模型选择的功能我们采用一个有限状态机（有限自动机）来实现，声明如下：

class FSM {
public:
	static int u, n, s, t, a, b, l, e, y, o,
		one, two, three, four, five, six, other, up, down , seven;//对应输入字母表
	static int transition[36][20];//对应转移函数
	bool attenuation, stable;//是否选择光照衰减，是否是stable类
	int select;//选择的模式
	int state;//当前状态
	float factor;//衰减因子
	FSM() : attenuation(false), stable(false), select(1), state(0) , factor(0.0005){}
	int nextState(char ch);//下一状态
};

通过上面这个有限自动机，对每一个键盘字符输入进行下一状态转移，再配合一定的选择结构（回调函数中的结构）就可以实现模型选择功能了。

渲染过程

渲染过程以代码展开
main函数结构：

int main(int argc , char* argv[])
{
    glutInit(&argc , argv);//初始化glut库
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);//选择显示模式，这里选择的是双缓冲，rgba颜色表示

    glutInitWindowSize(W_WIDTH, W_HIGHT);//设置窗口大小
    glutInitWindowPosition(100, 100);//窗口位置
    glutCreateWindow("多星系统");//创建一个以“多星系统”为标题的窗口

    myinit();

    glutReshapeFunc(myReshape);
    glutDisplayFunc(display);
    glutTimerFunc(1, move, 0);
    glutKeyboardFunc(processNormalKeys);
    glutSpecialFunc(processSpecialKeys);
//    glutPassiveMotionFunc(passiveMotion);

    glutMainLoop();

}

先介绍回调函数，所谓回调函数就是被别人（函数）自动调用的函数。glutReshapeFunc()注册一个回调函数在窗口形状发生变化时自动调用。glutDisplayFunc()注册一个回调函数在窗口收到重绘消息的时候自动调用。glutTimerFunc()注册一个回调函数，过了设定的时间就自动执行一次该回调函数。glutKeyboardFunc()注册一个回调函数在输入普通字符时自动调用。glutSpecialFunc()注册一个回调函数在输入特殊字符时自动调用。glutMainLoop()让所有与事件相关的函数无限循环，也就是说会一直检测事件，并且自动调用回调函数。
下面给出这几个函数的函数体：

void myinit(void) {//初始化一些基本设置
    glEnable(GL_DEPTH_TEST);//开启深度测试

    glEnable(GL_LIGHTING);//开启光照模式

    glShadeModel(GL_SMOOTH);//平滑明暗处理
}

void myReshape(GLsizei w, GLsizei h) {
    glViewport(0, 0, w, h);//定义视窗口

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(90, (GLfloat)w / h, 1, 2000);//定义投影模式

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    gluLookAt(camera.eye().X(), camera.eye().Y(), camera.eye().Z(),
        camera.center().X(), camera.center().Y(), camera.center().Z(),
        camera.up().X(), camera.up().Y(), camera.up().Z()
    );//在哪看
}
void move(int value) {
    S.simulate();//模拟一次短时间变化
    
    glutPostRedisplay();//发送重绘消息
    glutTimerFunc(1, move, 0);//再次调用
}

void display(void) {//显示的重要函数
    glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
    glClearDepth(1.0f);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);//清空屏幕

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glEnable(GL_LIGHT1);
    glLightfv(GL_LIGHT1, GL_POSITION, lightpositon);
    glLightfv(GL_LIGHT1, GL_AMBIENT, lightambient);
    glLightfv(GL_LIGHT1, GL_DIFFUSE, lightdiffuse);
    glLightfv(GL_LIGHT1, GL_SPECULAR, lightspecular);//开灯

    if (myFsm.attenuation == true) {
        glLightf(GL_LIGHT1, GL_CONSTANT_ATTENUATION, 1);
        glLightf(GL_LIGHT1, GL_LINEAR_ATTENUATION, 0);
        glLightf(GL_LIGHT1, GL_QUADRATIC_ATTENUATION, myFsm.factor);
    }
    else{
        glLightf(GL_LIGHT1, GL_CONSTANT_ATTENUATION, 1);
        glLightf(GL_LIGHT1, GL_LINEAR_ATTENUATION, 0);
        glLightf(GL_LIGHT1, GL_QUADRATIC_ATTENUATION, 0);
    }//衰减设置

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    gluLookAt(camera.eye().X(), camera.eye().Y(), camera.eye().Z(),
        camera.center().X(), camera.center().Y(), camera.center().Z(),
        camera.up().X(), camera.up().Y(), camera.up().Z()
    );//观察模型

    
    int i;
    for (i = 0; i < S.n; i++) {
        S.balls[i].Draw();
    }//向渲染管线输入数据
    glutSwapBuffers();//交换缓冲区
}


void processSpecialKeys(int key, int x, int y) {
    switch (key) {//按键控制视角
    case GLUT_KEY_UP: {
        camera.nod(alpha);
        break;
    }
    case GLUT_KEY_DOWN: {
        camera.nod(-alpha);
        break;
    }
    case GLUT_KEY_LEFT: {
        camera.shade(-alpha);
        break;
    }
    case GLUT_KEY_RIGHT: {
        camera.shade(alpha);
        break;
    }
    }
    lightpositon[0] = camera.eye().X();
    lightpositon[2] = camera.eye().Y();
    lightpositon[3] = camera.eye().Z();
//    cout << camera << std::endl;
    glutPostRedisplay();
}


void processNormalKeys(unsigned char key, int x, int y) {
    switch (key) {//按键控制视角
    case 'w': {
        camera.forward(step);
        break;
    }
    case 's': {
        camera.forward(-step);
        break;
    }
    }
    lightpositon[0] = camera.eye().X();
    lightpositon[2] = camera.eye().Y();
    lightpositon[3] = camera.eye().Z();
    if (myFsm.nextState(key) == 1) {//按键选择模式
        camera.reset();
        switch (myFsm.stable){
        case true: {
            S = stable_S[myFsm.select - 1];
            break;
        }
        case false: {
            S = unstable_S[myFsm.select - 1];
            break;
        }
        }
    }
    glutPostRedisplay();
}

之所以这样设计，是为了让各个模块的功能更加独立，例如mycamera的移动与旋转并没有直接涉及到操作，只是在操作后得到它的输入，这就使得既便是换了一种操作方法，只要有对的转化方法变成这个类的输入，仍然可以实现视角控制，而不需要该一大堆代码来实现操作模式的变换。同理FSM类也是如此。Ball类和System类更是独立于操作之外。这样结构化的设计程序，对应功能做成一个模块，这样程序的逻辑结构会非常明确，会使得程序在调试的时候思路更加明确：每个结构模块的独立的排除错位。再次强调一边所有代码都在附录的链接里。¹