TI-Nspire版Pythonによる振り子のシミュレーション / ばね振り子

(https://qiita.com/cyq04000/items/a006b2d216f0fdec5afd の続き )    
参考: http://www.epii.jp/articles/note/physics/analytic/spring_pendulum      

#概要
今度はばね振り子をシミュレートしてみる。

#サンプルファイル

```example_pend_spring.py
from PEND_SPRING_Class import *
from CREATE_SIM_SPACE_Class import *
import math
import ti_draw as td

# 画面の大きさを決める。
WIDTH, HIGHT = td.get_screen_dim()
ASPECT = WIDTH / HIGHT
XMIN = -15
XMAX = -XMIN
YMIN = XMIN / ASPECT
YMAX = XMAX / ASPECT

# シミュレートすべきオブジェクトをいくつか作る。
#(自然長,ばね定数,錘質量,錘密度,初期角度,初期角速度,初期長さ,ステップ,ステップ分割数) 
h = 1/2**4
num_of_div = 2**3
p0 = PEND_SPRING_Class(2, 20, 2, 10, math.pi*0.9, 0, 0.2, h, num_of_div)
p1 = PEND_SPRING_Class(2, 20, 6, 10, math.pi*0.9, 0, 0.2, h, num_of_div)

# オブジェクトをシミュレーション空間に投入し、シミュレーションを始める。
CREATE_SIM_SPACE_Class([p0, p1], XMIN, XMAX, YMIN, YMAX).start_sim()
```

##実行結果
[![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193936/7df8b850-ac29-0896-632f-171b20e0fcbf.png)](https://youtu.be/wdXJKTwFJ-o)



#ばね振り子の質点のxy座標を計算するためのクラス

<details><summary>`ソースコード (PEND_SPRING_Class.py) を見る/隠す`</summary><div>

```PEND_SPRING_Class.py
from RK4_Class import *
import time
import ti_draw as td
import math

PI = math.pi

def draw_spring(x0, y0, x1, y1, N=5, r=0.2):
  num_of_points = N * 2 + 2
  last_index = num_of_points - 1

  length = math.sqrt((x1 - x0)**2 + (y1 - y0)**2)
  theta = math.atan2(y1 - y0, x1 - x0)

  N4 = N * 4
  x = [x0 if (i==-1) else (x0 + length * i / N4) for i in range(-1, N4+2, 2)]
  x[last_index] = x0 + length

  y0_plus_r = y0 + r
  y0_minus_r = y0 - r
  y = [y0_plus_r if (i%2) else y0_minus_r for i in range(num_of_points)]
  y[0] = y0
  y[last_index] = y0

  cos_theta = math.cos(theta)
  sin_theta = math.sin(theta)
  offset_x = cos_theta * x0 + sin_theta * y0 + x0
  offset_y = cos_theta * y0 - sin_theta * x0 + y0
  xlist = [cos_theta * x[i] - sin_theta * y[i] - offset_x for i in range(num_of_points)]
  ylist = [sin_theta * x[i] + cos_theta * y[i] - offset_y for i in range(num_of_points)]  

  td.draw_poly(xlist, ylist)

class PEND_SPRING_Class:
  # (自然長,ばね定数,錘質量,錘密度,初期角度,初期角速度,初期長さ,ステップ,ステップ分割数) 
  def __init__(self, len0, k, m, dens,   init_angle, init_vel, init_len, h=1/2**5, numOfDiv=1):
    dens = 2 # 重りの密度〔kg/m^3〕。
    vol = m/dens # 重りの体積〔m^3〕。
    self.MASS_R = ((3 * vol)/(4 * math.pi))**(1/3) # 重りが球体の場合の半径〔m〕。

    # リセット用に引数を取っておく。
    self.ARGs = [len0, k, m, dens,   init_angle, init_vel, init_len, h, numOfDiv]

    self.LENGTH = init_len
    self.K = k

    #更新すべき初期値を設定する。
    self.t0 = 0
    self.radInit    = init_angle
    self.radDotInit = init_vel
    self.lenInit    = init_len
    self.lenDotInit = 0
    self.inits = [self.radInit, self.radDotInit, self.lenInit, self.lenDotInit]

    # 解くべき微分方程式を定義する。
    def radDot(t, rad, radDot, len, lenDot):
      return radDot
    def radDotDot(t, rad, radDot, len, lenDot):
      return (-2 * lenDot * radDot - 9.80665 * math.sin(rad)) / len
    def lenDot(t, rad, radDot, len, lenDot):
      return lenDot
    def lenDotDot(t, rad, radDot, len, lenDot):
      return len * radDot**2 + 9.80665 * math.cos(rad) - k * (len - len0) / m
    self.FUNCS = [radDot, radDotDot, lenDot, lenDotDot]

    # 上記の微分方程式を上記の初期値で解くための函数を実体化する。
    self.sol = RK4_Class(self.FUNCS, self.inits, h, numOfDiv)

  def update(self):
    try:
      self.sol.update() # 初期値を更新して、
      self.t0 = self.sol.t0 # 更新された時間を取り出して、
      self.radInit = self.sol.inits[0] # 更新された角度を取り出して、
      self.lenInit = self.sol.inits[2] # 更新された長さを取り出す。
      return self
    except:
      td.draw_text(0, 0, "0 DIV OCCURRED!!")

  def draw(self):
    x = self.lenInit * math.sin(self.radInit) # 角度から質点のx座標を計算する。
    y = self.lenInit * -math.cos(self.radInit) # 角度から質点のy座標を計算する。
    #td.draw_line(0, 0, x, y) # 原点からその質点まで線を引く。
    draw_spring(0, 0, x, y, r=self.K/50) # 原点からその質点までギザギザ線を引く。
    td.draw_circle(x, y, self.MASS_R) # 質点の位置に円を描く。
    return self

  def _print(self):
    print(self.t0, self.radInit, self.lenInit) # 計算結果の確認用。

########
# test #
########
#if __name__ == "__main__":
def main():
   pend = PEND_SPRING_Class(2, 20, 2, 20,  math.pi*0.9, 0, 0.5, h=1/2**5, numOfDiv=1)
   for i in range(10):
     pend.update()._print()
    
```

</div></details>