# Solve a set of second-order differential equations

#dy/dt = F(x,y,t) y=ω
#dx/dt =G(x,y,t)　x=θ

# 重力加速度 g [m/s^2]
g=9.8 

# 竿の長さ　L [m]
L=1.000

# 錘の質量 m [kg]
m=0.15

# 粘性減衰係数 c [Ns]
c=0.000464

# 微分方程式の右辺
F(x,y,t)= -(g/L)*sin(x)-c/(m*L)*y
G(x,y,t)=y

# Runge-Kuttaステップの時間幅とステップ数
dt = 0.001
N  = 10000

# xの初期値 x0 [rad]
x0=0.03

# y=dx/dtの初期値 y0 [rad/s]
y0=0 

# Runge-Kutta法の変数の初期化
k1=0.0; k2=0.0; k3=0.0; k4=0.0
m1=0.0; m2=0.0; m3=0.0; m4=0.0

outputfile="differential-eq-output.dat"

set samples N
set xrange [0:dt*(N-1)]

# いったん別ファイルにデータを書き出す
set table outputfile
plot xold=x0,yold=y0,\
     "+" using ($1+dt):\
     (k1=dt*F(xold,yold,$1), \
     m1 =dt*G(xold,yold,$1), \
     k2 =dt*F(xold+m1*0.5,yold+k1*0.5, $1+dt*0.5), \
     m2 =dt*G(xold+m1*0.5,yold+k1*0.5, $1+dt*0.5), \
     k3 =dt*F(xold+m2*0.5,yold+k2*0.5, $1+dt*0.5), \
     m3 =dt*G(xold+m2*0.5,yold+k2*0.5, $1+dt*0.5), \
     k4 =dt*F(xold+m3,yold+k3, $1+dt), \
     m4 =dt*G(xold+m3,yold+k3, $1+dt), \
     ynew=yold+(1/6.0)*(k1+2.0*k2+2.0*k3+k4),\
     xnew=xold+(1/6.0)*(m1+2.0*m2+2.0*m3+m4),\
     xold=xnew, \
     yold=ynew, \
     xnew)

unset table

# x軸y軸の範囲指定
set xrange [0:10]
set yrange [-0.05:0.05]

# x軸y軸のラベル付け
set xlabel "t  [sec]"
set ylabel "{/=30 {/Symbol-Oblique q}}  [rad]"

# グラフ描画
plot outputfile u 1:2 w l linecolor -1 linewidth 2 \
     notitle,0 w l linecolor -1 notitle

#出力フォーマットとオプションの指定
set terminal jpeg

#出力ファイル名の指定
set output "pendulum-sample.jpeg"

#グラフの再描画
replot

#フォーマットと出力のリセット
set output
set terminal win
[EOF]