#rad
"Simulation of simple damped pendulum by Runge-Kutta 4 method
'<table><tr><td>
'<h4>Input parameters</h4>
'Gravitational acceleration -'g = 9.80665m/s^2
'Pendulum length -'l = 1m
'Pendulum mass -'m = 1kg
'Initial angle -'θ_0 = -179°|rad
'Initial angular velocity -'ω_0 = 1.834rad/s
'Maximum simulation time -'t_max = 20s
'Damping factor (normaized) -'β = 0.005
'Critical damping -'b_cr = 2*m*sqrt(g/l)
'Damping -'b = β*b_cr
'</td><td>
#hide
Precision = 10^-5
PlotWidth = 300','PlotHeight = 150','PlotSVG = 1
w = 1.2*l', 'h_ = 2/3*w
W = 180', 'H = 150
k = W/w
L = l*k
R = L/10
x = L*sin(θ_0)
y = L*cos(θ_0)
δ = 16
δ_x = 0.25*δ
δ_y = -2*δ
y_g = y - δ_y + R
#def svg$ = '<svg viewbox="'-20-W/2' '-H-20' 'W+40' '2*H+40'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Georgia Pro; font-size:14px; width:'W+50'pt; height:'H+50'pt; margin-left:20pt;">
#def thin_style$ = style="stroke:goldenrod; stroke-width:1; fill:none"
#def thick_style$ = style="stroke:steelblue; stroke-width:2; fill:none"
#def ball_style$ = style="stroke:steelblue; stroke-width:0.5;  fill:url(#ball)"
#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
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; l*k; 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.2*δ*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;if(θ_0<0°;450°+θ_0;450°-θ_0);'l'm)
circle$(0; 0; R/2; solid_style$)
circle$(0; 0; R/8; solid_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> 
'<h4>Solution by Runge-Kutta RK4 method (explicit)</h4>
#noc
'Differential equation -'θ″ + b/m*θ′ + g/l*sin(θ) = 0
#equ
'The II order equation is first reduced to the following system of I order equations
#noc
θ′ = ω'and'ω′ = (g/l)*sin(-θ) - b/m*ω
#equ
'The solution is performed iteratively
'Step size -'h = 0.02s
'Number of steps -'n = t_max/h
#hide
N = min(n; 250)
k_N = floor(n/N)
N = n\k_N
#show
'For each time step'n'the values for the next step will be obtained by using the following equstions
#noc
'<table class="eq" style="margin-left:24pt;">
'<tr><td>First step (k₁) -</td><td>'k_1,θ = ω_i', </td><td>' _
k_1,ω = (g/l)*sin(-θ_i) - b/m*k_1,θ'</td></tr>
'<tr><td>Second step (k₂) -</td><td>'k_2,θ = ω_i + 0.5*h*k_1,ω', </td><td>'k_2,ω = (g/l)*sin(-(θ_i + 0.5*h*k_1,θ)) - b/m*k_2,θ'</td></tr>
'<tr><td>Third step (k₃) -</td><td>'k_3,θ = ω_i + 0.5*h*k_2,ω', </td><td>'k_3,ω = (g/l)*sin(-(θ_i + 0.5*h*k_2,θ)) - b/m*k_3,θ'</td></tr>
'<tr><td>Fourth step (k₄) -</td><td>'k_4,θ = ω_i + h*k_3,ω', </td><td>'k_4,ω = (g/l)*sin(-(θ_i + h*k_3,θ)) - b/m*k_4,θ'</td></tr>
'</table>
'Update values using weighted averages
'<p class="eq" style="border-left:solid 1pt black; margin-left:24pt; padding-left:3pt;"><var>θ</var><sub>n+1</sub> ='θ_i + (h/6)*(k_1,θ + 2*k_2,θ + 2*k_3,θ + k_4,θ)'<br/> <var>ω</var><sub>n+1</sub> = 'ω_i + (h/6)*(k_1,ω + 2*k_2,ω + 2*k_3,ω + k_4,ω)'</p>
#equ
'Allocate vectors
θ = vector(n)'<br/>'ω = vector(n)
'Set initial conditions
θ.1 = θ_0/1rad','ω.1 = ω_0/1rad
'Perform Runge-Kutta 4 steps 
#for i = 1 : n - 1
	'RK4 factors
	k_1,θ = ω.i', ' _
	k_1,ω = (g/l)*sin(-θ.i) - b/m*k_1,θ
	k_2,θ = ω.i + 0.5*h*k_1,ω
	k_2,ω = (g/l)*sin(-(θ.i + 0.5*h*k_1,θ)) - b/m*k_2,θ
	k_3,θ = ω.i + 0.5*h*k_2,ω
	k_3,ω = (g/l)*sin(-(θ.i + 0.5*h*k_2,θ)) - b/m*k_3,θ
	k_4,θ = ω.i + h*k_3,ω
	k_4,ω = (g/l)*sin(-(θ.i + h*k_3,θ)) - b/m*k_4,θ
	'Update values using weighted averages
	'Rotation -'θ.(i + 1) = θ.i + (h/6)*(k_1,θ + 2*k_2,θ + 2*k_3,θ + k_4,θ)
	'Angular velocity -'ω.(i + 1) = ω.i + (h/6)*(k_1,ω + 2*k_2,ω + 2*k_3,ω + k_4,ω)' < /p > 
	#hide
#loop
t(i) = i*h - h
r2d = 180/π
#show
'Kinetic energy -'E_k(i) = m/2*l^2*ω.i^2|J
'Potential energy-'E_p(i) = m*l*g*(1 - cos(θ.i))|J
'...
'<h4>Plot results</h4>
'Rotation - <var>θ</var>(<var>t</var>), [°]
$Plot{t(i)|θ.i*r2d @ i = 1 : n}
'Energy - <var>E</var>(<var>t</var>), [J]
$Plot{t(i)|E_k(i) + E_p(i) & 0s|0J & 0s|0J & 0s|0J & t(i)|E_p(i) & t(i)|E_k(i) & 0s|0J @ 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>
'<style>#start{height:1.8em; padding-top:0;}.Series3, .Series4{stroke:none!Important; fill:none!Important;} .fr{display:none;} .fr svg{margin-right:64pt; align-center;}</style>
'<h4>Animate results &emsp;<button id="start">Start</button></h4>
#format 0.00
#hide
k = W/w
#val
#for i_N = 1 : N
	#hide
	i = i_N*k_N
	θ_i = θ.i
	x = L*sin(θ_i)
	y = L*cos(θ_i)
	t = t(i)
	#show
	#val
	'<div class="fr" id="fr-'i_N'">
	#equ
	#nosub
	'Simulation time -'t'<br/>
	$Plot{t(i)|θ.i*r2d & t(j)|θ.j*r2d & t_max|θ_0/1rad*r2d & t_max|-θ_0/1rad*r2d @ j = 1 : i}
	#val
	svg$
	'<polygon points="0,0 '-δ','-δ' 'δ','-δ'" solid_style$/>
	line$(0; 0; x; y; thick_style$)
	circle$(0; 0; R/2; solid_style$)
	circle$(0; 0; R/8; solid_style$)
	circle$(x; y; R; ball_style$)
	'</svg>
	'</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