Python之美数值模拟抛体运动

发表于 更新于

2023年05月18日星期四多云北京市北京师范大学,顺利完成Python之美作业，由我完成了Linux版本，程国栋负责完成Windows版本。

设计思路

  1. 据Ip获取本地城市名
  2. 根据城市由高德地图获取纬度
  3. 由纬度计算本地的g值
  4. 输入初速度大小
  5. 输入阻尼大小
  6. 绘制最大水平位移与夹角的关系
  7. 绘制海平面存在风阻时的模拟图像
  8. 绘制海平面上100米处存在风阻时的模拟图像
  9. 绘制海平面上200米处存在风阻时的模拟图像

代码实现

nspp_linux.py 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
 
#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
# Author:        冯振华 & 程国栋 
# StudentNumber: 202221140015 & 202221140011
# Version:       v2.2  for linux
# Linux:         冯振华
# Windows:       程国栋 
# Date:          2023年05月18日星期四多云北京市北京师范大学
# Copyright © 2023 feng <feng@archlinux>
#
# !!注意：本脚本在linux编写，vim编辑器下
#       :set ff      #查看文件格式
#       :set ff=dos  #设置为dos格式
#       :set ff=unix #设置为unix格式
# 如果不设置正确的格式地，运行程序会报错
#
import math
import matplotlib.pyplot as plt
import os
import requests
import json
# 输入初速度大小和空气阻力系数
v_0 = input("请输入初速度大小: ")
v_0 = float(v_0)
k = input("请输入空气阻力系数: ")
k = float(k)
# 获取本机Ip对应的地理位置
Ip_data=os.popen('curl -s  cip.cc').readlines()
Ip_local=Ip_data[1].split(' ')
city=Ip_local[-1].strip()
# 从高德地图获取纬度
GaoDeKey = '79e3576b80e0b544d0d17dd643e4654e'
url = 'https://restapi.amap.com/v3/geocode/geo'
params = {'key': GaoDeKey,
    'address': city}
res = requests.get(url,params)
jd = json.loads(res.text)
coords = jd['geocodes'][0]['location']
jd_value = coords.split(',')
latitude=float(jd_value[1])
latitude = latitude/180*math.pi
# 根据纬度计算本地的重力加速度g
g=9.7803*(1+0.005324*((math.sin(latitude))**2)-0.000005*((math.sin(2*latitude))**2))
# 计算最大水平位移和初速度与水平方向夹角的关系
def x_max(theta):
    rad = math.radians(theta)
    return (v_0**2* math.sin(2*rad)/ g)
    angles = range(0, 91, 1)
    x = [x_max(angle) for angle in angles]
#with the drag force
def simulate(v_1, theta,y_height):
    theta = math.radians(theta)
    v_x = v_1 * math.cos(theta)
    v_y = v_1 * math.sin(theta)
    x = 0
    y = float(y_height)
    t = 0
    delta_t = 0.01
    traj = [(x, y)]
    while y >= 0:
        a_x = -k * v_x ** 2 / (2 * v_1)
        a_y = -g - k * v_y ** 2 / (2 * v_1)
        v_x += a_x * delta_t
        v_y += a_y * delta_t
        x += v_x * delta_t
        y += v_y * delta_t
        t += delta_t
        traj.append((x, y))
    return traj
# 将各种情况集中于一张图上
plt.figure(figsize=(200,200),dpi=100)
plt.subplots_adjust(left=None,bottom=None,right=None,top=None,wspace=None,hspace=0.5)
plt.figure(1)
ax1 = plt.subplot(221)
ax1.margins(0.15)
ax1.plot(angles, x)
ax1.set_title("$\\theta$ with $x_{max}$ Free")
ax1.set_xlabel("$\\theta$(degrees)")
ax1.set_ylabel("$x_{max}$(meters)")
ax2 = plt.subplot(222)
ax2.margins(0.15)
angles = range(0, 91, 10)
for theta in angles:
    traj = simulate(v_0, theta,0) # 地表
    x_s = [t[0] for t in traj]
    y_s = [t[1] for t in traj]
    ax2.plot(x_s, y_s, label=f"${theta}^\circ$")
ax2.legend(loc="upper left")
ax2.set_title(f"Projectile Motion with Air Resistance ($v_0={v_0} m/s$)")
ax2.set_xlabel("Distance (meters)")
ax2.set_ylabel("Height (meters)")
ax3 = plt.subplot(223)  
ax3.margins(0.15)
angles = range(0, 91, 10)
for theta in angles:
    traj = simulate(v_0, theta,1)  # 100m高塔
    x_s = [t[0] for t in traj]
    y_s = [t[1] for t in traj]
    ax3.plot(x_s, y_s, label=f"${theta}^\circ$")
ax3.legend(loc="upper left")
ax3.set_title(f"Projectile Motion with Air Resistance ($v_0={v_0} m/s$)")
ax3.set_xlabel("Distance (meters)")
ax3.set_ylabel("Height (meters)")
ax4 = plt.subplot(224)
ax4.margins(0.15)
angles = range(0, 91, 10)
for theta in angles:
    traj = simulate(v_0, theta,2)  # 200m高塔
    x_s = [t[0] for t in traj]
    y_s = [t[1] for t in traj]
    ax4.plot(x_s, y_s, label=f"${theta}^\circ$")
ax4.legend(loc="upper left")
ax4.set_title(f"Projectile Motion with Air Resistance ($v_0={v_0} m/s$)")
ax4.set_xlabel("Distance (meters)")
ax4.set_ylabel("Height (meters)")
plt.show()
nspp_windows.py 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
 
#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
# Author:        冯振华 & 程国栋 
# StudentNumber: 202221140015 & 202221140011
# Version:       v2.2  for windows
# Linux:         冯振华
# Windows:       程国栋 
# Date:          2023年05月18日星期四多云北京市北京师范大学
# Copyright © 2023 feng <feng@archlinux>
#
# !!注意：本脚本在linux编写，vim编辑器下
#       :set ff      #查看文件格式
#       :set ff=dos  #设置为dos格式
#       :set ff=unix #设置为unix格式
# 如果不设置正确的格式地，运行程序会报错
#
import math
import numpy as np
import matplotlib.pyplot as plt
import requests
import json
from bs4 import BeautifulSoup
from scipy.integrate import odeint


# 输入初速度大小和空气阻力系数
v_0 = float(input("请输入初速度大小: "))
k = float(input("请输入空气阻力系数: "))


# 获取本机Ip对应的地理位置
url = 'http://www.cip.cc/'
header = {'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:109.0) Gecko/20100101 Firefox/112.0'}
html = requests.get(url, headers=header)  # 访问网址
bsObj = BeautifulSoup(html.content, 'lxml')  # 使用lxml解析器
msg = bsObj.find('div', class_="data kq-well")
msg = [i for i in msg]
city = msg[1]


# 从高德地图获取纬度
GaoDeKey = '79e3576b80e0b544d0d17dd643e4654e'
url = 'https://restapi.amap.com/v3/geocode/geo'
params = {'key': GaoDeKey,
          'address': city}
res = requests.get(url, params)
jd = json.loads(res.text)
coords = jd['geocodes'][0]['location']
jd_value = coords.split(',')
latitude = float(jd_value[1])
latitude = latitude/180*math.pi  # 转换成弧度


# 根据纬度计算本地的重力加速度g
g = 9.7803*(1+0.005324*((math.sin(latitude))**2)-0.000005*((math.sin(2*latitude))**2))


def simulate(v_1, theta, y_height):
    # 轨迹图
    theta = math.radians(theta)
    v_x0 = v_1 * math.cos(theta)
    v_y0 = v_1 * math.sin(theta)

    def motion1(w, t, k):
        v_x, v_y = w
        return[-k * v_x ** 2, -g - k * v_y ** 2]
    t = np.linspace(0, 20, 300000)
    track1 = odeint(motion1, (v_x0, v_y0), t, args=(k, ))
    v_x1 = [i for i in track1[:, 0]]
    v_y1 = [j for j in track1[:, 1] if j > 0]
    v_x1 = v_x1[:len(v_y1)]
    delta_t = (20 - 0) / 300000
    i = 1
    x = 0
    y = y_height
    y_list1 = [y]
    x_list1 = [0]
    while i < len(v_y1) - 1:
        x += ((v_x1[i] + v_x1[i-1]) / 2) * delta_t
        y += ((v_y1[i] + v_y1[i-1]) / 2) * delta_t
        if y < 0:
            break
        x_list1.append(x)
        y_list1.append(y)
        i += 1

    def motion2(w, t, k):
        v_x, v_y = w
        return[-k * v_x ** 2, -g + k * v_y ** 2]

    track2 = odeint(motion2, (v_x1[-1], 0), t, args=(k, ))
    v_x2 = [i for i in track2[:, 0]]
    v_y2 = [j for j in track2[:, 1] if j < 0]
    j = 1
    x = x_list1[-1]
    y = y_list1[-1]
    y_list2 = [y]
    x_list2 = [x]
    while j < len(v_y2) - 1:
        x += ((v_x2[j] + v_x2[j-1]) / 2) * delta_t
        y += ((v_y2[j] + v_y2[j-1]) / 2) * delta_t
        if y < 0:
            break
        x_list2.append(x)
        y_list2.append(y)
        j += 1
    x_list = x_list1 + x_list2
    y_list = y_list1 + y_list2
    return x_list, y_list


def x_MaxTheta(v, height):
    x_max = []
    for angle in range(1, 91, 1):
        x_list, y_list = simulate(v_1=v, theta=angle, y_height=height)
        x_max.append(x_list[-1])
    return x_max


heght = float(input('输入初始高度:'))
x_max = x_MaxTheta(v=v_0, height=heght)
plt.plot(range(1, 91, 1), x_max, 'g-')
plt.show()
x = 0
j = 1
for i in x_max:
    j += 1
    if i > x:
        x = i
        num = j - 1
print('当前初始条件下抛得最远的角度是{},最远距离是{}'.format(num, x))
print('输入要模拟的轨迹的抛出角度:')
a = float(input())
x, y = simulate(v_1=v_0, theta=a, y_height=heght)
plt.plot(x, y, 'r-')
plt.show()