# なだれ 仕様
大まかな処理の流れはデータの処理の後、演算と描画を繰り返すだけ。
この繰り返しの1回を1Tickとする。
データ処理では実験で得たデータから演算に必要な確率を出す。
演算では各粒子の位置を計算し、更新する。物理系の処理はここに入る。
描画では各粒子の位置をもとに画像を生成し、出力する。
<details>
<summary>ソースコード</summary>
```python=
import cv2
import json
import os
import random
import math
import csv
import numpy as np
path = os.getcwd()
with open(path + "/setting.json") as f:
print("Loading settings...")
settings = json.load(f)
print(settings)
def magnitude(array: list[float]):
"""与えられた配列を成分としたベクトルの大きさを返す。
:param array: ベクトルの成分
:typep array: list
:returns: ベクトルの大きさ
:rtype: float
"""
sum = 0
for a in array:
sum += a ** 2
sum = math.sqrt(sum)
return sum
def getComponent(component, magnitude):
return math.sqrt(magnitude ** 2 - component ** 2)
#def getAngle()
path = os.getcwd()
with open(path + "/data/standard_deviations1.csv") as f:
reader = csv.reader(f)
temp = [row for row in reader]
sd = temp[0]
sd = [float(v) for v in sd]
alpha = temp[1]
alpha = [float(v) for v in alpha]
beta = temp[2]
beta = [float(v) for v in beta]
#print(sd[-1])
#print(random.normalvariate(mu=0, sigma=sd[45]))
def horizontalAngle(angle):
"""入射角がangle度の時の水平方向の反射角を確率を加味して返す。
Parameters
----------
angle : int
垂直方向の入射角
Returns
-------
float
水平方向の角度の変化
"""
return random.normalvariate(mu=0, sigma=sd[round(angle)])
def verticalAngle(angle):
"""入射角がangle度のときの垂直方向の反射角を確率を加味して返す。
Parameters
----------
angle : int
垂直方向の入射角
Returns
-------
float
水平方向の角度の変化
"""
return random.gammavariate(alpha=alpha[round(angle)], beta=beta[round(angle)])
#for i in range(100):
# print(horizontalAngle(30))
class Particle:
tickrate = settings["physics"]["tickrate"]
#重力加速度
ga = settings["physics"]["gravity_acceleration"]
#反発係数
rc = settings["physics"]["repulsion_coefficient"]
def __init__(self, posX: float, posY: float, posZ: float):
self.posX = posX
self.posY = posY
self.posZ = posZ
self.motionX = 0
self.motionY = settings["physics"]["initial_motionY"]
self.motionZ = 0
self.horizontalAngle = 90
self.verticalAngle = settings["physics"]["initial_vertical_angle"]
def update(self):
"""自身の位置を更新する。
"""
self.motionZ -= self.ga * self.tickrate
if self.posZ <= 0:
a = horizontalAngle(self.verticalAngle)
b = verticalAngle(self.verticalAngle)
# 水平方向の衝突
self.horizontalAngle += a
mag = magnitude([self.motionX,self.motionY])
self.motionX = mag * math.cos(math.radians(self.horizontalAngle))
self.motionY = mag * math.sin(math.radians(self.horizontalAngle))
self.motionZ = self.motionZ * -1
# 垂直方向の衝突
self.verticalAngle = b
if b > 90:
b -= 90
self.motionX *= -1
self.motionY *= -1
mag = magnitude([self.motionX, self.motionY, self.motionZ])
self.motionX = mag * math.cos(math.radians(self.horizontalAngle)) * math.cos(math.radians(self.verticalAngle))
self.motionY = mag * math.sin(math.radians(self.horizontalAngle)) * math.cos(math.radians(self.verticalAngle))
self.motionZ = mag * math.sin(math.radians(self.verticalAngle))
# 非弾性衝突
self.motionX *= self.rc
self.motionY *= self.rc
self.motionZ *= self.rc
# print(a,self.verticalAngle,mag)
self.posX += self.motionX * self.tickrate
self.posY += self.motionY * self.tickrate
self.posZ += self.motionZ * self.tickrate
def setMotion(self,x,y,z):
self.motionX = x
self.motionY = y
self.motionZ = z
def render(particles: list[Particle]):
"""particleを描画する。
Parameters
----------
particles : list[Particle]
描画する粒子(Particle)の配列
"""
height = settings["rendering"]["height"]
width = settings["rendering"]["width"]
radius = settings["rendering"]["particle_radius"]
color = settings["rendering"]["particle_rgb"]
img = np.full((height, width, 3),255, dtype=np.uint8)
#circle(img,中心,半径,色,<アンチエイリアスをいれるならlineType=cv2.LINE_AA>)
# cv2.circle(img, (100,100), settings["particle_radius"],(settings["particle_rgb"][2], settings["particle_rgb"][1], settings["particle_rgb"][0]),ineType=cv2.LINE_AA)
for p in particles:
posX = round(p.posX)
posY = round(p.posY)
cv2.circle(img, (posX,posY), radius, (color[2],color[1],color[0]), lineType=cv2.LINE_AA )
cv2.imshow("test", img)
cv2.moveWindow("test",0,0)
supporting_img = np.full((height, 300,3),255, dtype=np.uint8)
for p in particles:
posX = round(p.posZ)
posY = round(p.posY)
cv2.circle(supporting_img, (posX,posY), radius, (color[2],color[1],color[0]), lineType=cv2.LINE_AA)
cv2.imshow("supporting_img", supporting_img)
cv2.moveWindow("supporting_img", width, 0)
# cv2.waitKey(0)
# return False
testparticle = Particle(500,0,10)
particles = []
particles.append(testparticle)
def update():
"""メインループ。"""
for p in particles:
if magnitude([p.motionX, p.motionY, p.motionZ]) < settings["physics"]["minimum_speed"]:
particles.remove(p)
break;
p.update()
if p.posZ < 0:
for i in range(settings["physics"]["bound_quantity"]):
a = Particle(p.posX, p.posY, p.posZ)
a.motionX = p.motionX
a.motionY = p.motionY
a.motionZ = p.motionZ
a.verticalAngle = p.verticalAngle
particles.append(a)
# print(magnitude([p.motionX, p.motionY, p.motionZ]))
if __name__ == "__main__" :
while cv2.waitKey(1) != 27:
update()
render(particles)
```
</details>
## 要検討
- [ ] 速度、距離の尺度、単位はどうするか 重力加速度も
- [x] ~~ベクトルの成分から角度を出すのはどうすればいいの?~~
成分と角度の変換をするのではなく、衝突前後のベクトルの大きさが変わらないことを利用して求める。
- [x] ~~XY面で衝突後の角度がn度(nは整数)になる確率はどう出すか
正規分布の特定の範囲の確率を出せることは確認した。でもこの場合複数の条件それぞれに正規分布があるから、ただの正規分布として確率を出せない。~~
n度の時において、実験のデータから平均、標準偏差の近似を取ることにする。
- [x] ガンマ分布？確率は出せる?
## TODO
- [x] データ処理実装
- [x] 実験1の標準偏差を集める
- [x] 正規分布から確率を出す
- [ ] 物理系実装
- [ ] 定数決定 場所は作った。
- [x] 衝突
- [x] XY面
- [x] Z軸
- [x] ベクトル計算
- [ ] この文章上でのいろいろな値に共通した文字を置く
- [x] すべての変数、関数、クラスのドキュメント(pydoc)を書く
## データ処理
**Pythonでは実装せず、Excelで確率を計算している。**
**確率を出すExcelファイルは`prob_calc.xlsx`。これから直接確率を読み込んでいるわけではないので注意。**
~~実験データはCSV形式で一つのファイルに1種の実験を保存しておく。~~
~~実験データの必要な値だけ取り出したものを保存する。確率まではExcelで出し、それをCSV形式で保存することにした。~~
正規分布があればpython標準のrandomが使えることがわかった。そのため標準偏差だけを`standard_deviations1.csv`に保存している。
一行目は正規分布を出すための標準偏差。
二行目はガンマ分布を出すためのパラメータしーた。
三行目はガンマ分布を出すためのパラメータけー。
間違っているかもしれないところ。
限られた条件でしか実験していないため、実験していない値はある程度の予測(近似?)で求める。
下のグラフのように一回目の実験の標準偏差は指数関数で近似を取るのが最も適していそう。
---
### 指数関数で近似
断片的で誤差を含む関数を指数関数( $y=a^{x}$ )とみなす。以下はその手順。
実験1の40cmのデータは以下のようになる。
| 角度 | 標準偏差 | 自然対数 |
| ------ | -------- | -------- |
| 0 | | |
| ...... | | |
| 15 | 7.51 | 2.02 |
| ...... | | |
| 30 | 8.94 | 2.19 |
| ...... | | |
| 45 | 16.94 | 2.83 |
| ...... | | |
| 60 | 16.53 | 2.80 |
| ...... | | |
| 75 | 32.15 | 3.47 |
| ...... | | |
| 90 | | |
最初に値の自然対数を取る。自然対数( $\ln{x}$ )とはネイピア数( $e$ )を底とする対数( $\log_e{x}$ )。
$e$ を使うと微分積分が楽になるらしい。[Wikipedia ネイピア数](https://ja.wikipedia.org/wiki/%E3%83%8D%E3%82%A4%E3%83%94%E3%82%A2%E6%95%B0)
値が指数関数で近似する場合、値の自然対数は一次関数に近似する(一次近似、線形近似)。これは以下のように表せる。
$$y=ae^{bx}$$
$$\begin{aligned}
\ln{y}
&=\ln{ae^{bx}}\\
&=\ln{a}+bx
\end{aligned}$$
次に自然対数に近似する一次関数を最小二乗法を用いて求める。データ $\left ( x,y \right )$ が $n$ 組与えられたとき、最小二乗法による直線の式は以下のようになる。
$$y=Ax+B$$
$$A=\frac{Cov(X,Y)}{\sigma_X^{2} }\\
B=\mu_Y-A\mu_X$$
ただし、 $Cov\left ( X,Y \right )$ は共分散、 $\mu_X$ 、 $\sigma_X$は $x$の平均と標準偏差、 $\mu_Y$ は $y$ の平均。
そうすれば $A=0.023488121$ $B=1.605480755$ という値が求められる(ExcelのLINEST関数を用いた)。
上2つの式から $A=b$ 、 $B=\ln{a}$より
$$a=e^{1.605480755}=4.980253310$$
$$b=0.023488121$$
最後に $y=ae^{bx}$ に代入して
$$y=4.980253310e^{0.023488121x}$$
が求めていた標準偏差に近似する指数関数になる!!
というのをexcelでやったらできた。
微妙にずれている。
---
### 衝突後の角度がmになるの確率Pを出す
~~以下は現状使っておらず確率はExcelを用いて求めている。Excelの`NORM.DIST`の仕様を理解する必要がある。~~
今はpythonの`random.normalvariate()`をこれに用いているため、以下は必要ない。この関数の仕様を理解する必要がある。[random.normalvariate](https://docs.python.org/ja/3/library/random.html#random.normalvariate)
1. 衝突前の角度がn度のときの標準偏差を取る。
1. 標準偏差と平均(実験1の場合は常に0とみなす)を用いて標準正規分布を取る。
標準正規分布は $\boldsymbol{Z}= \frac{\boldsymbol{X}-\boldsymbol{m}}{\sigma }$ と表せる。 $\boldsymbol{m}$ は平均、 $\sigma$ は標準偏差になる。 $\boldsymbol{X}$ は確率変数。
1. 数B P.80 例題3のようにPを求める。m-0.5からm+0.5の範囲にXがある確率をmとみなす。
$$\begin{aligned}
P\left ( m-0.5< X\leq m+0.5 \right )\\
&= P\left ( \frac{\left ( m-0.5 \right )-m}{\sigma }< Z\leq \frac{\left ( m-0.5 \right )-m}{\sigma }m+0.5 \right )\\
&=P\left ( \frac{-0.5}{\sigma }< Z\leq \frac{0.5}{\sigma } \right )\\
&=P\left ( 0\leq Z \leq \frac{0.5}{\sigma } \right ) - P\left ( 0\leq Z \leq \frac{-0.5}{\sigma } \right )\\
&=p\left ( \frac{0.5}{\sigma } \right )-p\left ( \frac{-0.5}{\sigma } \right )
\end{aligned}$$
この式書くのに30分かかった。けどあってる自信がまったくない。
---
pythonでは以下の計算をしないことにした。
~~実行のたびに計算し直すのは良くないので計算結果はどこかに保存しておく。~~
~~計算するコードを別に用意してメインにこの処理を含めないのがいいかも。~~
~~計算で出すのは~~
- ~~角度がn度(nは整数)になる確率~~
## 数学?系
### ベクトル(Vector)
現状使っていない。
よく考えたらこのクラスって必要か?
X、Y、Zの成分を持つ。
`changeComponents(x,y,z)`でベクトルの大きさを保ったまま成分を変える。引数をすべて指定すると大きさが変わってしまうので1つか2つ指定するように。この実装は微妙。
成分3つを変える関数は作らないつもり。新しくインスタンスしてください。
よく考えたらこれって必要か?
### mathutil
数学系の計算をする。
`magnitude()`は与えられた配列の各要素を2乗した和の平方根を返す。ベクトルの大きさを取るときとかに使う。
`getComponent()`第一引数との`magnitude`が第二引数になる数を返す。ベクトルの大きさを変えずに成分を変えるときに使う。名前が微妙。
## 物理系
Z-UPのXYZ座標上とする。
画面右向きがX軸正方向、画面下向きがY軸正方向。
画面左上が原点。
粒子は位置、各軸方向の速度、~~速度のベクトル~~を持つ。
速度のベクトルは都度出したほうがいいかも。
粒子はZ座標を持つが描画はXY面しかされない。
定数としてsettings.jsonに反発係数(`repulsion_coefficient`)を定義する。~~軸によって異なるようにするのがいいかもしれない。~~ おそらく必要ない。
定数として重力加速度(`gravity_acceleration`)を定義する。
定数として1回の更新に時間をどれくらい進めるか(`tickrate`)を定義する。
XY面とZ軸で分けて書いているが分けて処理できないと思われる。うまく実装する必要がある。
速度が`minimum_speed`を超えると消える。
### XY面上の動き
粒子は面に当たったとき(=Z座標が0のとき)XY方向の速度が変わる。
それ以外のときは速度を保って等速直線運動する。
---
面との衝突の処理
条件: Z<=0
1. 実験データから得た確率と乱数を用いてで衝突後のXY方向の速度の角度を決める。
1. 速度のベクトルが等しくなるようにX、Y軸の速度を変更する。
1. X、Y軸の速度に反発係数をかける。
1. 通常処理へ
---
通常処理
条件: なし
毎Tick実行する。
1. X、Y座標を各軸の速度とTickrateをかけた値にする。
### Z軸方向の動き
粒子は面に当たったとき(=Z座標が0のとき)Z方向の速度が変わる。
それ以外のときは等加速度直線運動する。
ここらへんあってるかどうか怪しい。
---
面との衝突の処理
条件: Z<=0
~~XY面での衝突より先に処理する必要があるかも。~~
ない。同時にしなければならない。
1. 実験データから得た確率と乱数を用いて衝突後の角度(?)を求める。 ここでは少数の値が出る。角度を求めるには整数しか使えないので四捨五入をしている。
1. 速度のベクトルが等しくなるようにX、Y、Zの速度を変更する。
1. Z軸の速度に反発係数と-1をかける。
1. 通常処理へ
---
通常処理
条件: なし
毎Tick実行する。
1. Z速度に重力加速度`tickrate`を足す。
1. Z座標をZ軸の速度と`tickrate`をかけた値にする。
↑どっちが先にするのがいいの?
## 描画系
描画するときは少数は使えない。描画のときだけ四捨五入する。
とりあえず直前に四捨五入した値を使う。まとめたほうが良さそうだったら一つにまとめる。
Sign in with Wallet
Connect another wallet