# 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 I [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
h=c/k
# viscous damping factor j
j=c/(2*sqrt(I*k))

# frequency f
f=6

# natural frequency wn
wn=f*pi

# impulse torque of balance M [Nm]
M=2.18e-7
m=M/k

# angle of balance at the start point of unlocking a [degree]
a=-25 #degree

# angle of balance at the end point of unlocking b [degree]
b=-15 #degree

# angle of balance at the start point of impulse g [degree]
g=-13 #degree

# angle of balance at the end point of impulse q [degree]
q=22 #degree

# ratio of unlock torque and impulse torque of balance e
e=0.3

# time width and step numbers of Runge-Kutta step
dt = 0.000001
N = 340000

# 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 time shift  s0 [s]
s0=0

# right side of differential equation 
# (changing the equation of motion here)
F(x,y,t)= y > 0 ? (x < a*pi/180 ? -wn**2*(h*y \
+x+sgn(yold)*r)\
:(x < b*pi/180 ? -wn**2*(h*y+x+sgn(yold)*r \
+sgn(yold)*e*m)\
:(x < g*pi/180 ? -wn**2*(h*y+x+sgn(yold)*r) \
:(x < q*pi/180 ? -wn**2*(h*y+x+sgn(yold)*r \
-sgn(yold)*m)\
: -wn**2*(h*y+x+sgn(yold)*r)))))\
:(x > -a*pi/180 ? -wn**2*(h*y \
+x+sgn(yold)*r)\
:(x > -b*pi/180 ? -wn**2*(h*y+x+sgn(yold)*r \
+sgn(yold)*e*m)\
:(x > -g*pi/180 ? -wn**2*(h*y+x+sgn(yold)*r) \
:(x > -q*pi/180 ? -wn**2*(h*y+x+sgn(yold)*r \
-sgn(yold)*m)\
:-wn**2*(h*y+x+sgn(yold)*r)))))
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(ynew)*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(ynew)*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:0.34]
set yrange [-6:6]

# setting the aspect ratio of graph
set size ratio 1

# setting the grid of graph
set grid

# 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
#plot 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 "escapement-error.jpeg"

# regraphing
replot

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