PyOpenGL三体模拟

给定多星系统的初始状态,以一定的时间步,计算在引力作用下的星体运动,并使用openGL实时可视化。

实验环境

python37+OpenGL https://www.cnblogs.com/GraceSkyer/p/9235582.html

numpy

PIL

初始条件

使用一个数组p表示多星系统的初始条件

p = [[1.0, 0, 0, 0, 0, 0, 0],[1.0, 0, 10, 0, 0.2, 0, 0],[1.0, 10, 0, 0, 0, 0, 0.2],]

数组每行表示一个星体的:质量,x坐标,y坐标,z坐标,x速度,y速度,z速度

单位分别是:kg,m,m,m,m/s,m/s,m/s

新建openGL运行类NBody,传入初始条件

class NBody:def __init__(self, node, k, time_step, width=640, height=480, title='NBody'.encode()):self.ptcN = len(node)self.path = [deque() for _ in range(self.ptcN)]self.ptcM = node[:, 0]self.ptcP = node[:, 1:4]self.ptcV = node[:, 4:]self.step = time_stepself.K = k

状态更新

计算每个node受到其他node的引力的合力,并以此计算node在时间步内的位置变化。

for item in range(self.ptcN):force = np.array([0.0, 0, 0])for i in range(self.ptcN):if i == item:continuedistance = self.ptcP[item] - self.ptcP[i]ss = 1/(sum(np.square(distance))**0.5+0.01)force -= ss**3*self.ptcM[i]*self.K*distanceself.ptcV[item] += force
self.ptcP += self.ptcV*self.step

可视化

初始化

self.color = []
for i in range(self.ptcN):self.color.append(np.random.randint(min(2 * i % 6, 2 * (i + 1) % 6), max(2 * i % 6, 2 * (i + 1) % 6)))
self.run = 1
glutInit(sys.argv)
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH)
glutInitWindowSize(width, height)
self.window = glutCreateWindow(title)
glutDisplayFunc(self.Draw)
glutMouseFunc(self.myMouse)
glutMotionFunc(self.onMouseMove)
glutSpecialFunc(self.mySpecial)
glutIdleFunc(self.Draw)
self.InitGL(width, height)
self.xR = 0.0
self.yR = 0.0
self.Oldx = 0
self.Oldy = 0
self.scale = 1.0
self.tz = 50

绘制画布,先计算新的状态,再根据进行视角旋转,再绘制每个node,记录node位置的历史数据,根据历史数据在绘制运动轨迹。

def Draw(self):if self.run:self.refresh()glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)glLoadIdentity()glEnable(GL_DEPTH_TEST)glDepthFunc(GL_LEQUAL)glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST)# 沿z轴平移gluLookAt(0, 0, -self.tz,  # 相机在世界坐标系中的位置坐标0, 0, 0,  # 相机正对着的世界坐标系中的点的位置坐标,成像后这一点会位于画板的中心位置0, 1, 0)  # 定义相机本身的朝向# 分别绕x,y,z轴旋转glRotatef(self.xR, 1.0, 0.0, 0.0)glRotatef(self.yR, 0.0, 1.0, 0.0)glScalef(self.scale, self.scale, self.scale)glPushMatrix()for i, (x, y, z) in enumerate(self.ptcP):glBindTexture(GL_TEXTURE_2D, self.color[i])glTranslatef(x, y, z)glutSolidSphere((self.ptcM[i] + 1) ** 0.3, 10, 10)glTranslatef(-x, -y, -z)if np.random.randint(10) == 3:self.path[i].append([x, y, z])if len(self.path[i]) > 40:self.path[i].popleft()glBegin(GL_LINES)for t in range(len(self.path[i]) - 1):glVertex3f(self.path[i][t][0], self.path[i][t][1], self.path[i][t][2])glVertex3f(self.path[i][t + 1][0], self.path[i][t + 1][1], self.path[i][t + 1][2])glEnd()glPopMatrix()glFlush()# 刷新屏幕,产生动画效果glutSwapBuffers()# time.sleep(0.001)

鼠标控制,实现鼠标拖动时视角旋转。

def myMouse(self, btn, state, x, y):if btn == GLUT_LEFT_BUTTON and state == GLUT_DOWN:self.Oldx, self.Oldy = x, yif state == GLUT_UP and btn == 3:self.scale += 1if state == GLUT_UP and btn == 4:self.scale -= 1glutPostRedisplay()def onMouseMove(self, x, y):self.yR += x - self.Oldxself.yR = min(self.yR, 180.0)self.yR = max(self.yR, -180.0)glutPostRedisplay()self.Oldx = xself.xR -= y - self.Oldyself.xR = min(self.xR, 180.0)self.xR = max(self.xR, -180.0)glutPostRedisplay()self.Oldy = y

键盘控制,实现视角缩放和暂停

def mySpecial(self, key, x, y):if key == GLUT_KEY_UP:self.tz -= 1if key == GLUT_KEY_DOWN:self.tz += 1if key == GLUT_KEY_F4:self.run = 1 - self.runglutPostRedisplay()def myIdle(self):pass

加载纹理,每个node的纹理使用jpg纹理文件加载

def LoadTexture(self):imgFiles = ['./color/' + str(i) + '.jpg' for i in range(6)]for i in range(6):img = Image.open(imgFiles[i])width, height = img.sizeimg = img.tobytes('raw', 'RGBX', 0, -1)glGenTextures(2)glBindTexture(GL_TEXTURE_2D, i)glTexImage2D(GL_TEXTURE_2D, 0, 4,width, height, 0, GL_RGBA,GL_UNSIGNED_BYTE, img)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_S, GL_CLAMP)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_T, GL_CLAMP)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_S, GL_REPEAT)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_T, GL_REPEAT)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER, GL_NEAREST)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER, GL_NEAREST)glTexEnvf(GL_TEXTURE_ENV,GL_TEXTURE_ENV_MODE, GL_DECAL)

设置InitGL

def InitGL(self, width, height):self.LoadTexture()glEnable(GL_TEXTURE_2D)glClearColor(1.0, 1.0, 1.0, 0.0)glClearDepth(1.0)glDepthFunc(GL_LESS)glShadeModel(GL_SMOOTH)# 背面剔除,消隐glEnable(GL_CULL_FACE)glCullFace(GL_BACK)glEnable(GL_POINT_SMOOTH)glEnable(GL_LINE_SMOOTH)glEnable(GL_POLYGON_SMOOTH)glMatrixMode(GL_PROJECTION)glLoadIdentity()gluPerspective(45.0, float(width) / float(height), 0.1, 100.0)glMatrixMode(GL_MODELVIEW)

运行效果

在这里插入图片描述

可实现功能:

  • 对N体问题进行模拟
  • 随机选择预设的纹理,进行可视化展示
  • 使用鼠标控制视角变化
  • 对历史轨迹可视化

总结

基于openGL实现的N体问题模拟。在低速,较少节点下有很好的运行效果。如果节点较多可以考虑使用八叉树进行优化。基本实现了openGL的基本图形显示,视角控制和简单纹理功能。

三体问题(three-body problem)是天体力学中的基本力学模型。它是指三个质量、初始位置和初始速度都是任意的可视为质点的天体,在相互之间万有引力的作用下的运动规律问题。现已知三体问题不能精确求解,即无法预测所有三体问题的数学情景,只有几种特殊情况已研究。三体问题最简单的一个例子就是太阳系中太阳、地球和月球的运动。在浩瀚的宇宙中,星球的大小可以忽略不记,所以我们可以把它们看成质点。如果不计太阳系其他星球的影响,那么它们的运动就只是在引力的作用下产生的,所以我们就可以把它们的运动看成一个三体问题。研究三体问题的方法大致可分为分析方法、定性方法、数值方法三类。

我看到了我的爱恋
我飞到她的身边
我捧出给她的礼物
那是一小块凝固的时间
时间上有美丽的条纹
摸起来像浅海的泥一样柔软

她把时间涂满全身
然后拉起我飞向存在的边缘
这是灵态的飞行
我们眼中的星星像幽灵
星星眼中的我们也像幽灵

完整代码

import sys
from OpenGL.GL import *
from OpenGL.GLUT import *
from OpenGL.GLU import *
from PIL import Image
import numpy as np
from collections import deque
import timeclass NBody:def __init__(self, node, k, time_step, width=640, height=480, title='NBody'.encode()):self.ptcN = len(node)self.path = [deque() for _ in range(self.ptcN)]self.ptcM = node[:, 0]self.ptcP = node[:, 1:4]self.ptcV = node[:, 4:]self.step = time_stepself.K = kself.color = []for i in range(self.ptcN):self.color.append(np.random.randint(min(2 * i % 6, 2 * (i + 1) % 6), max(2 * i % 6, 2 * (i + 1) % 6)))self.run = 1glutInit(sys.argv)glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH)glutInitWindowSize(width, height)self.window = glutCreateWindow(title)glutDisplayFunc(self.Draw)glutMouseFunc(self.myMouse)glutMotionFunc(self.onMouseMove)glutSpecialFunc(self.mySpecial)glutIdleFunc(self.Draw)self.InitGL(width, height)self.xR = 0.0self.yR = 0.0self.Oldx = 0self.Oldy = 0self.scale = 1.0self.tz = 50def refresh(self):for item in range(self.ptcN):force = np.array([0.0, 0, 0])for i in range(self.ptcN):if i == item:continuedistance = self.ptcP[item] - self.ptcP[i]ss = 1 / (sum(np.square(distance)) ** 0.5 + 0.001)force -= ss ** 3 * self.ptcM[i] * self.K * distanceself.ptcV[item] += forceself.ptcP += self.ptcV * self.stepdef setup(self):glClearColor(1.0, 1.0, 1.0, 1.0)glShadeModel(GL_SMOOTH)ambient = [0.4, 0.6, 0.2, 1.0]position = [1.0, 1.0, 3.0, 1.0]glEnable(GL_LIGHTING)glEnable(GL_LIGHT0)glLightfv(GL_LIGHT0, GL_AMBIENT, ambient)glLightfv(GL_LIGHT0, GL_POSITION, position)mat_ambient = [0.5, 0.2, 0.5, 1.0]mat_diffuse = [0.8, 0.6, 0.3, 1.0]mat_specular = [0.8, 0.6, 0.3, 1.0]mat_shininess = [45.0]glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient)glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse)glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular)glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess)glEnable(GL_COLOR_MATERIAL)glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE)glDepthFunc(GL_LESS)glEnable(GL_DEPTH_TEST)glEnable(GL_AUTO_NORMAL)glEnable(GL_NORMALIZE)# 绘制图形def Draw(self):if self.run:self.refresh()glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)glLoadIdentity()glEnable(GL_DEPTH_TEST)glDepthFunc(GL_LEQUAL)glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST)# 沿z轴平移gluLookAt(0, 0, -self.tz,  # 相机在世界坐标系中的位置坐标0, 0, 0,  # 相机正对着的世界坐标系中的点的位置坐标,成像后这一点会位于画板的中心位置0, 1, 0)  # 定义相机本身的朝向# 分别绕x,y,z轴旋转glRotatef(self.xR, 1.0, 0.0, 0.0)glRotatef(self.yR, 0.0, 1.0, 0.0)glScalef(self.scale, self.scale, self.scale)glPushMatrix()for i, (x, y, z) in enumerate(self.ptcP):glBindTexture(GL_TEXTURE_2D, self.color[i])glTranslatef(x, y, z)glutSolidSphere((self.ptcM[i] + 1) ** 0.3, 10, 10)glTranslatef(-x, -y, -z)if np.random.randint(10) == 3:self.path[i].append([x, y, z])if len(self.path[i]) > 40:self.path[i].popleft()glBegin(GL_LINES)for t in range(len(self.path[i]) - 1):glVertex3f(self.path[i][t][0], self.path[i][t][1], self.path[i][t][2])glVertex3f(self.path[i][t + 1][0], self.path[i][t + 1][1], self.path[i][t + 1][2])glEnd()glPopMatrix()glFlush()# 刷新屏幕,产生动画效果glutSwapBuffers()# time.sleep(0.001)def myMouse(self, btn, state, x, y):if btn == GLUT_LEFT_BUTTON and state == GLUT_DOWN:self.Oldx, self.Oldy = x, yif state == GLUT_UP and btn == 3:self.scale += 1if state == GLUT_UP and btn == 4:self.scale -= 1glutPostRedisplay()def onMouseMove(self, x, y):self.yR += x - self.Oldxself.yR = min(self.yR, 180.0)self.yR = max(self.yR, -180.0)glutPostRedisplay()self.Oldx = xself.xR -= y - self.Oldyself.xR = min(self.xR, 180.0)self.xR = max(self.xR, -180.0)glutPostRedisplay()self.Oldy = ydef mySpecial(self, key, x, y):if key == GLUT_KEY_UP:self.tz -= 1if key == GLUT_KEY_DOWN:self.tz += 1if key == GLUT_KEY_F4:self.run = 1 - self.runglutPostRedisplay()def myIdle(self):passdef LoadTexture(self):# 提前准备好的6个图片imgFiles = ['./color/' + str(i) + '.jpg' for i in range(6)]for i in range(6):img = Image.open(imgFiles[i])width, height = img.sizeimg = img.tobytes('raw', 'RGBX', 0, -1)glGenTextures(2)glBindTexture(GL_TEXTURE_2D, i)glTexImage2D(GL_TEXTURE_2D, 0, 4,width, height, 0, GL_RGBA,GL_UNSIGNED_BYTE, img)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_S, GL_CLAMP)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_T, GL_CLAMP)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_S, GL_REPEAT)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_T, GL_REPEAT)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER, GL_NEAREST)glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER, GL_NEAREST)glTexEnvf(GL_TEXTURE_ENV,GL_TEXTURE_ENV_MODE, GL_DECAL)def InitGL(self, width, height):self.LoadTexture()glEnable(GL_TEXTURE_2D)glClearColor(1.0, 1.0, 1.0, 0.0)glClearDepth(1.0)glDepthFunc(GL_LESS)glShadeModel(GL_SMOOTH)# 背面剔除,消隐glEnable(GL_CULL_FACE)glCullFace(GL_BACK)glEnable(GL_POINT_SMOOTH)glEnable(GL_LINE_SMOOTH)glEnable(GL_POLYGON_SMOOTH)glMatrixMode(GL_PROJECTION)glLoadIdentity()gluPerspective(45.0, float(width) / float(height), 0.1, 100.0)glMatrixMode(GL_MODELVIEW)def MainLoop(self):glutMainLoop()if __name__ == '__main__':# weight, position_x, y, z, velocity_x, y, zp = [[1.0, 0, 0, 0, 0, 0, 0],[0, 0, 10, 0, 0.2, 0, 0],[0, 10, 0, 0, 0, 0, 0.2],[0, 0, 0, 10, 0, 0.2, 0],[0, 0, -10, 0, -0.2, 0, 0],[0, -10, 0, 0, 0, 0, -0.2],[0, 0, 0, -10, 0, -0.2, 0],]p = np.array(p)w = NBody(p, 0.1, 0.1)w.MainLoop()

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/45969.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

使用python进行字频统计和词频统计

问题描述 读取给定的语料库,根据制表符’\t’划分其文本与标签,将获得的文本仅保留汉字部分,并按字划分,保存在列表中,至少使用一种方法,统计所有汉字的出现次数,并按照从高到低的顺序排序&…

Imagenet VGG-19图片识别实例展示

资源: 1.相关的vgg模型下载网址 http://www.vlfeat.org/matconvnet/models/beta16/imagenet-vgg-verydeep-19.mat 2.ImageNet 1000种分类以及排列 https://github.com/sh1r0/caffe-Android-demo/blob/master/app/src/main/assets/synset_words.txt(如果…

外滩画报:揭秘全球电子垃圾坟墓

在西方发达国家,有这样一个不为人知的秘密:当你把电子垃圾送给回收商而不是扔进垃圾箱里后,很快,大约 80% 的电子垃圾就会被装上集装箱船,运往尼日利亚、印度、巴基斯坦和中国那些常年被毒烟笼罩的垃圾场。师从人道主义…

《荒野猎人》影评

如果没有小李子和奥斯卡数十年相爱相杀赚足眼球这件事,《荒野猎人》最大的看点应该是导演亚利桑德罗冈萨雷斯伊纳里图和摄影师艾曼努尔卢贝兹基的再度合作。 伊纳里图的电影履历如晴朗夏夜的星空一般漂亮璀璨,执导座《爱情是狗娘》一举拿下2000年东京国际…

印度之行(一) 印度是个很大的国家

(baidu真渣。。。) 首先,我们弄清了,老王吃了一个月咖喱的地方,下面是老王的咖喱味思考: 在印度浦内市(又译浦那市)待了5周,基本是在公司的浦内office做项目,…

HTML标签

HTML标签友情链接 如果在 HTML 中需要文字或者图片垂直居中时&#xff0c;可以使用 align "center"&#xff0c;可以对文字或着图片进行垂直居中&#xff01; 这是一个很好的例子&#xff1a; <!doctype html> <html lang"en"> <head&…

苏州免费景点

文章目录 苏州免费景点姑苏区相门古城墙平江历史街区曲园畅园七里山塘朴园五峰园拥翠山庄天香小筑明轩实样北寺塔定慧寺城隍庙神仙庙苏州博物馆忠王府苏州民俗博物馆中国昆曲博物馆苏州公园桂花公园桐泾公园 吴中区沐春园(原园博会苏州园)瑞园道勤小筑乡畦小筑石湖风景区灵岩山…

又一个程序员被判刑了!运维违规操作被判5年半,IT从业需要懂法律!

点击上方 "程序员小乐"关注, 星标或置顶一起成长 每天凌晨00点00分, 第一时间与你相约 每日英文 Sleeping is nice. You forget about everything for a little while. 还是睡觉好&#xff0c;能暂时忘掉一切烦恼。 每日掏心话 人生如梦&#xff0c;岁月匆匆&#x…

又背锅了,一个“锁表” 损失 800万,程序员被判5年半

图片来自 Pexels 出处&#xff1a;云头条&#xff0c;知乎综合整理, 51CTO 链接&#xff1a;https://www.zhihu.com/question/389167387/answer/1170852426 近日&#xff0c;云头条发布的“一个违规操作、损失 800 万、被判五年半&#xff1a;运维夏某某致郑大一附院智慧医院系…

一个“锁表”损失800万,运维被判5年半

“ 近日&#xff0c;云头条发布的“一个违规操作、损失 800 万、被判五年半&#xff1a;运维夏某某致郑大一附院智慧医院系统瘫痪 2 个小时&#xff0c;判破坏计算机信息系统罪”一文引发了技术圈的热议。 图片来自 Pexels 事件经过 夏某某任职北京中科某某科技有限公司&#x…

大运河的一些总结

京杭大运河&#xff0c;是一条承载着历史&#xff0c;流淌着过去和未来的河流&#xff0c;是我非常喜欢的一条人工河&#xff0c;但一直以来不太清楚其具体分段和水流流向&#xff0c;近期通过了解南水北调东线工程&#xff0c;并经过查阅资料终于弄清楚了京杭大运河的各个河段…

科技周刊第六期:接近本质的东西才会长远

这里记录每周值得分享的东西&#xff0c;每周五发布。 封面图 中国西南西藏自治区山南市扎南县的雅鲁藏布江&#xff08;出处&#xff09; 本周话题&#xff1a;接近本质的东西才会长远 我想说三个现象&#xff1a; 1、为什么很多明星能够一直红下去&#xff1f;而有的明星只…

python调用mysql并在前台做数据展示

今天是学习python的第二天。 根据自己的需要&#xff0c;将前段时间的扇形图稍微升华一下&#xff0c;从而可以从mysql数据库中查询数据&#xff0c;并作图形的展示。 以下为图形展示&#xff1a; #导入库--注意本段代码不适用于python2 import pymysql import matplotlib.py…

中国农民丰收节交易会暨“日照有礼”功能性特色产品展示

中国农民丰收节交易会暨“日照有礼”功能性特色产品展示 齐鲁网闪电新闻讯 新闻中国采编网 中国新闻采编网 谋定研究中国智库网 经信研究 国研智库 国情讲坛 万赢信采编&#xff1a;为庆祝农民丰收节&#xff0c;进一步推动日照乡村振兴建设步伐&#xff0c;倡导功能性农业农业…

C++ STL的简单运用——学习记录

本周学了一下STL&#xff0c;学的内容有string&#xff0c;stack&#xff0c;queue&#xff0c;vector&#xff0c;priority queue&#xff0c;map&#xff0c;set与其简单应用&#xff0c; 速度比较快&#xff0c;&#xff0c;还好学之前我预习了一些&#xff0c;&#xff0c;…

ArcGIS水文分析实战教程(12)河网分级流程

ArcGIS水文分析实战教程(12)河网分级流程 本章导读&#xff1a;如果说河流提取是面对没有数据后者数据匮乏的用户&#xff0c;那么河网分级就完全属于为水文研究而生的一个工具。河流具有干流和支流之分&#xff0c;河网分级能够将这些干支关系理顺&#xff0c;并从中找到其水文…

[论文总结] 深度学习在农业领域应用论文笔记2

文章目录 1. A comparative study of fruit detection and counting methods for yield mapping in apple orchards &#xff08;IF3.581, 2019&#xff09;1.1 介绍1.2 实验与结果1.3 定性结果1.4 失败案例1.4.1 检测1.4.2 计数 1.5 结论与未来工作 2.Wheat crop yield predic…

Java 反射慢?它到底慢在哪?

往期热门文章&#xff1a; 1、GitHub 被超火的 ChatGPT 霸榜&#xff01; 2、Java使用 try catch会影响性能&#xff1f; 3、原来count(*)是接口性能差的真凶&#xff01; 4、大公司病了&#xff0c;这也太形象了吧&#xff01;&#xff01;&#xff01; 5、全球最大资源站创始…

运营人必懂 | TikTok运营指南

“TikTok之前确实很火&#xff0c;现在呢&#xff1f;” 最新数据告诉你&#xff1a; Sensor Tower商店情报数据显示&#xff0c;2022年9月抖音及海外版TikTok在全球App Store和Google Play吸金超过3.15亿美元&#xff0c;是去年同期的1.7倍&#xff0c;蝉联全球移动应用&…

抖音号运营爆量爆单技巧

泛知识付费2.0时代&#xff0c;短视频、直播间成为了知识传播的重要阵地&#xff0c;只要有技能干货&#xff0c;不论是行业大咖&#xff0c;还是精通某领域的普通人&#xff0c;都有机会成为大众的“老师”&#xff0c;依靠输出视频、音频等内容课程来知识变现。而抖音生态&am…