# Solve a set of second-order differential equations

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

# moment of inertia of balance wheel [Kgm^2]
I=1.4e-9

# spring constant of hairspring k [Nm/rad]
k=4.97428e-7

# solid friction force R [Nm]
R=3e-9
r=R/k

# viscous damping coefficient　c [Nms]
c=1e-10

# viscous damping factor j
j=c/(2*sqrt(I*k))

# frequency f
f=6

# natural frequency wn
wn=f*pi

# time width and step numbers of Runge-Kutta step
dt = 0.0001
N = 20000

# default value of x  x0 [rad]
x0=5.5

# default value of y=dx/dt y0 [rad/s]
y0=0

# default value of z (theoretical value) z0 [rad]
z0=5.5

# default value of the time shift s0 [s]
s0=0

# right side of differential equation
F(x,y,t)=-c*y/I-k/I*x-sgn(yold)*R/I
G(x,y,t)=y

#  initializing the variable of Runge-Kutta method
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)]

# writing down the data to separate file
set table outputfile
plot xold=x0,yold=y0,zold=z0,Aold=z0,sold=s0,\
"+" 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),\
ynew*yold <= 0 ? Anew=xnew : Anew=Aold,\
ynew*yold <= 0 ? snew=$1 : snew=sold,\
znew=(Anew+sgn(yold)*r)*exp(-j*wn*($1-snew)) \
*cos((sqrt(1-j**2))*wn*($1-snew) \
-atan(j/(sqrt(1-j**2))))/(sqrt(1-j**2))-sgn(yold)*r,\
xold=xnew,\
yold=ynew,\
zold=znew,\
Aold=Anew,\
sold=snew,\
xnew):(ynew):(znew):(Anew):(snew):(0) w xyerrorbars

unset table

# setting the range of x and y axes
set xrange [0:2]
set yrange [-6:6]

# setting the label of x and y axes
set xlabel "t [sec]"
set ylabel "θ [rad]"

# setting the kinds of lines of graph
set style line 1 linetype 1
set style line 1 linecolor -1
set style line 1 linewidth 2

set style line 2 linetype 1
set style line 2 linecolor 1
set style line 2 linewidth 1

# graphing
plot outputfile u 1:2 w l linestyle 1 title "rk-4",\
0 w l linecolor -1 notitle
replot outputfile u 1:4 w l linestyle 2 title "theoretical"

# designating the output format and option
set terminal jpeg

# designating the output file name
set output "temp.jpeg"

# regraphing
replot

# resetting the format and output
set output
set terminal win
[EOF]