# duck.ode
# the infamous drinking duck
# modeled as a pendulum with some extra forces
# M is the total mass of the two parts of the duck
# m1 is the lower mass and M-m1 the upper
# mmin is the mass of the two parts when empty
# thus, M-2*mmin is the mass of the fluid
# we will call mu the mass of th is fluid in the lower vessel
# without loss of generality we assume this is 1 so M=1+2*mmin
mudot=-eps*abs(xdot)*min(beta*mu,1)
m1=mmin+mu
m2=mmin+mutot -mu
mu'=mudot
x'=xdot
xdot'=(-f*xdot-(m1*l1-m2*l2)*sin(x)+k/(1-cos(x-xlo))-mudot*(l1^2-l2^2))/(m1*l1^2+m2*l2^2)
par l1=2.6,l2=2.5,mmin=2,f=.5,k=.002,xlo=-.2,eps=.015,xdip=1.6,mutot=.4
par beta=50
global 1 x-xdip {xdot=0;mu=mutot}
init mu=.25,x=1.5
# some animation constants
x1=.5+.8*(l1/(l1+l2))*sin(x)
y1=.5-.8*(l1/(l1+l2))*cos(x)
x2=.5-.8*(l2/(l1+l2))*sin(x)
y2=.5+.8*(l2/(l1+l2))*cos(x)
xb=x2-.15*cos(x)
yb=y2-.15*sin(x)
@ meth=qualrk,dt=1,total=1000
done