#rad
"Simulation of elastic damped pendulum by Runge-Kutta 4 method
'<table><tr><td>
'<h4>Input parameters</h4>
'Gravitational acceleration -'g = 9.80665m/s^2
'Pendulum unstretched length -'l = 1m
'Pendulum mass -'m = 1kg
'Strut stiffness -'k_s = 40*(N/.m)
'Initial angle -'θ_0 = -179.9°|rad
'Initial angular velocity -'ω_0 = 0rad/s
'Initial length -'r_0 = 1.1*l
'Initial radial velocity -'v_0 = 0m/s
'Damping factor (normalized, radial) -'β_r = 0.05
'Damping factor (normalized, angular) -'β_θ = 0.01
'</td><td>
#hide
PlotWidth = 220','PlotHeight = 140','PlotSVG = 1
w = 1.2*l', 'h_ = 2/3*w
W = 150', 'H = 150
k = W/w
L = r_0*k
R = L/12
x = L*sin(θ_0)
y = L*cos(θ_0)
δ = 16
δ_x = 0.25*δ
δ_y = -2*δ
y_g = y - δ_y + R
#def svg$ = '<svg viewbox="'-100 - W/2' '-H' 'W + 200' '2*H  + 100'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Georgia Pro; font-size:14px; width:'W + 100'pt; height:'H + 100'pt;margin-left:-10pt;">
#def thin_style$ = style="stroke:goldenrod; stroke-width:1; fill:none;"
#def thick_style$ = style="stroke:steelblue; stroke-width:1.5; fill:none;"
#def ball_style$ = style="stroke:steelblue; stroke-width:0.5;  fill:url(#ball);"
#def bal0_style$ = style="stroke:steelblue; stroke-width:0.25;  fill:url(#ball); fill-opacity:0.5;"
#def solid_style$ = style="stroke:black; stroke-width:1; fill:#eee;"
#def load_style$ = style="stroke:steelblue; stroke-width:1; fill:steelblue;"
#include svg_drawing.cpd
#show
#val
'<br/><br/>
svg$
'<defs><radialGradient id="ball" cx="35%" cy="35%"><stop offset="0%" stop-color="lightcyan"/><stop offset="100%" stop-color="lightblue"/></radialGradient><pattern id="concrete" x="4" y="4" width="8" height="8" patternUnits="userSpaceOnUse"><rect x="0" y="0" width="8" height="8" fill="#eae9e8" /><circle cx="1" cy="1" r="1.2" fill="#ccb" /><circle cx="5" cy="2" r="1.6" fill="#dadad0" /><circle cx="4" cy="6" r="0.8" fill="#aa9" /><circle cx="3" cy="4" r="0.4" fill="#884" /><circle cx="7" cy="5" r="1.2" fill="#cacaba" /><circle cx="5" cy="3" r="0.9" fill="#fffded" /></pattern></defs>
'<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
line$(0; 0; 0; 0.4*L; thin_style$)
'<path d="M 0 '0.2*L' A '0.2*L' '0.2*L' 0 0 'θ_0 < 0°' '0.2*x' '0.2*y'" thin_style$/>
texth$(-1.5*δ*sign(θ_0);0.25*L;'θ_0*180/π'°)
line$(0; 0; x; y; thick_style$)
text$(x/2-δ_x*y/L; y/2+δ_x*x/L;450°+θ_0;'r_0' ('l') m)
circle$(0; 0; R/2; solid_style$)
circle$(0; 0; R/8; solid_style$)
circle$(l*k*sin(θ_0); l*k*cos(θ_0); R; bal0_style$)
circle$(x; y; R; ball_style$)
'<polygon points="'x','y_g + δ_y' 'x','y_g' 'x - δ_x','y_g + δ_y/3' 'x + δ_x','y_g + δ_y/3' 'x','y_g'" load_style$/>
texth$(x+3*δ_x; y_g-δ;g)
'</svg>
#equ
'</td></tr></table>
'Critical damping (radial) -'b_cr_r = 2*m*sqrt(k_s/m)
'Critical damping (angular) -'b_cr_θ = 2*m*l*sqrt(g/l)
'Radial damping -'b_r = β_r*b_cr_r
'Angular damping -'b_θ = β_θ*b_cr_θ
'Maximum simulation time -'t_max = 20s
'<h4>Solution by Runge-Kutta RK4 method (explicit)</h4>
'Differential equations
'<div style="border-left:solid 1pt black; margin-left:24pt; padding-left:3pt;">
#noc
m*r″ - m*r*θ′^2 + k_s*(r - l) + b_r*r′ - m*g*cos(θ) = 0
m*r*θ″ + 2*m*r′*θ′ + b_θ*θ′ + m*g*sin(θ) = 0
#equ
'</div>
'The II order system of 2 equations is first reduced to the following I order system of 4 equations
'<div style="border-left:solid 1pt black; margin-left:24pt; padding-left:3pt;">
#noc
r′ = v
v′ = r*ω^2 - k/m*(r - l) - b_r/m*v + g*cos(θ)
θ′ = ω
ω′ = -((2*m*v*ω + b_θ*ω + m*g*sin(θ))/(m*r))
#equ
'</div>
'The solution is performed iteratively
'Step size -'h = 0.01s
'Number of steps -'n = t_max/h
#hide
N = min(n; 250)
k_N = floor(n/N)
N = n\k_N
#show
'Allocate vectors
r = vector(n)
v = vector(n)
θ = vector(n)
ω = vector(n)
'Set initial conditions
r.1 = r_0','θ.1 = θ_0/1rad
v.1 = v_0','ω.1 = ω_0/1rad
'Define the system of ODEs
f_r(r; v; θ; ω) = v
f_v(r; v; θ; ω) = r*ω^2 + g*cos(θ) - k_s/m*(r - l) - b_r/m*v
f_θ(r; v; θ; ω) = ω
f_ω(r; v; θ; ω) = -((2*m*v*ω + b_θ*ω + m*g*sin(θ))/(m*r))
'Perform Runge-Kutta 4 steps
#for i = 1 : n - 1
	r_i = r.i', 'v_i = v.i
	θ_i = θ.i', 'ω_i = ω.i
	k_1r = h*v_i', 'k_1v = h*f_v(r_i; v_i; θ_i; ω_i)
	k_1θ = h*ω_i', 'k_1ω = h*f_ω(r_i; v_i; θ_i; ω_i)
	r_i = r.i + k_1r/2', 'v_i = v.i + k_1v/2
	θ_i = θ.i + k_1θ/2', 'ω_i = ω.i + k_1ω/2
	k_2r = h*v_i', 'k_2v = h*f_v(r_i; v_i; θ_i; ω_i)
	k_2θ = h*ω_i', 'k_2ω = h*f_ω(r_i; v_i; θ_i; ω_i)
	r_i = r.i + k_2r/2', 'v_i = v.i + k_2v/2
	θ_i = θ.i + k_2θ/2', 'ω_i = ω.i + k_2ω/2
	k_3r = h*v_i', 'k_3v = h*f_v(r_i; v_i; θ_i; ω_i)
	k_3θ = h*ω_i', 'k_3ω = h*f_ω(r_i; v_i; θ_i; ω_i)
	r_i = r.i + k_3r', 'v_i = v.i + k_3v
	θ_i = θ.i + k_3θ', 'ω_i = ω.i + k_3ω
	k_4r = h*v_i', 'k_4v = h*f_v(r_i; v_i; θ_i; ω_i)
	k_4θ = h*ω_i', 'k_4ω = h*f_ω(r_i; v_i; θ_i; ω_i)
	'Update values for the next step
	r.(i + 1) = r.i + (k_1r + 2*k_2r + 2*k_3r + k_4r)/6
	v.(i + 1) = v.i + (k_1v + 2*k_2v + 2*k_3v + k_4v)/6
	θ.(i + 1) = θ.i + (k_1θ + 2*k_2θ + 2*k_3θ + k_4θ)/6
	ω.(i + 1) = ω.i + (k_1ω + 2*k_2ω + 2*k_3ω + k_4ω)/6
	#hide
#loop
t(i) = i*h - h
r2d = 180/π
#show
'Kinetic energy -'E_k(i) = m/2*(v.i^2 + (r.i*ω.i)^2)|J
'Potential energy -'E_p(i) = k_s/2*(r.i - l)^2 - m*g*(r.i*cos(θ.i) - l)|J
'...
'<h4>Plot results</h4><table><tr><td>Angle - <var>θ</var>(<var>t</var>), [°]<br/>
$Plot{t(i)|θ.i*r2d @ i = 1 : n}
'</td><td>Length - <var>r</var>(<var>t</var>), [m]<br/>
$Plot{t(i)|r.i @ i = 1 : n}
'</td></tr><tr><td>Trajectory - <var>x</var>(<var>t</var>), <var>y</var>(<var>t</var>), [m]<br/>
$Plot{r.i*sin(θ.i)|-r.i*cos(θ.i) @ i = 1 : n}
'</td><td>Energy - <var>E</var>(<var>t</var>), [J]<br/>
$Plot{t(i)|E_p(i) + E_k(i) & 0s|0J & 0s|0J & 0s|0J & t(i)|E_p(i) & t(i)|E_k(i) @ i = 1 : n}
'<br/><b style="color:MediumVioletRed">━━</b> Potential - <var>E</var><sub>p</sub>, 
'<b style="color:MediumSpringGreen">━━</b> Kinetic - <var>E</var><sub>k</sub>, 
'<b style="color:Tomato">━━</b> Total - <var>E</var><sub>tot</sub>
'</td></tr></table>
'<style>.Series3, .Series4{stroke:none!Important; fill:none!Important;} #start{height:1.8em; padding-top:0;} .fr{display:none;} .fr svg{align:center;} .trace{stroke:red; stroke-width:1; fill:red;}</style>
'<h4>Animate results &emsp;<button id="start">Start</button></h4>
#format 0.00
#hide
k = W/w
x = r*sin(θ)
y = r*cos(θ)
n_tr = 40
#for i_N = 1 : N
	#hide
	i = i_N*k_N
	xi = x.i*k
	yi = y.i*k
	t = t(i)
	i_0 = max((i_N - n_tr)*k_N; 2)
	n_0 = i - i_0
	Δr = abs(r_0 - l) + 0.2*l
	#show
	#val
	'<div class="fr" id="fr-'i_N'">
	#equ
	#nosub
	'Simulation time -'t'<br/>
	'<table><tr><td>Angle - <var>θ</var>(<var>t</var>), [°]<br/>
	$Plot{t(i)|θ.i*r2d & t(j)|θ.j*r2d & t_max|θ_0/1rad*r2d & t_max|-θ_0/1rad*r2d @ j = 1 : i}
	'<br/>Length - <var>r</var>(<var>t</var>), [m]<br/>
	$Plot{t(i)|r.i & t(j)|r.j & t_max|l + Δr & t_max|l - Δr @ j = 1 : i}
	'</td><td>
	#val
	svg$
	#for j = 1 : n_0/2
		#hide
		j2 = i_0 + 2*j
		j1 = j2 - 2
		o = round(j/n_0*100)/100
		#show
		'<line x1="'x.j1*k'" y1="'y.j1*k'" x2="'x.j2*k'" y2="'y.j2*k'" class="trace" opacity="'o'"/>
	#loop
	'<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
	line$(0; 0; xi; yi; thick_style$)
	circle$(0; 0; R/2; solid_style$)
	circle$(0; 0; R/8; solid_style$)
	circle$(xi; yi; R; ball_style$)
	'</svg></td></tr></table></div>
#loop
'<script>$("#fr-1").show();$("#start").click(function(){$(this).hide(); var i='N';var fr=$("#fr-"+i);setInterval(function(){fr.hide();if(++i>'N')i=1;fr=$("#fr-"+i);fr.show();}, '1000*h*k_N');});</script>
#equ