Let the distance between the centers of two identical ellipses with semi-axes $a, b$ be $a+b$. The first ellipse rotates around its center with a constant angular velocity $w$. The second ellipse rotates so that it is in contact with the first ellipse (i.e. there is a one point of contact between the two ellipses). Calculate the speed of the second ellipse depending on the rotation angle of the first ellipse. Find the locus of the points of contact of two ellipses.
To find the contact line, to solve the system of three equations $$x_1(t,p)=x_2(s,q)$$
$$y_1(t,p)=y_2(s,q)$$
$$\frac{x_1(t,p)_t}{x_2(s,q)_s}=\frac{y_1(t,p)_t}{y_2(s,q)_s}$$
Here $$p$$ and $$q$$ are the rotation angles of the ellipses. The system of equations is solved using standard Python tools.
The locus contact points are two lines - one in the form of an infinity symbol (1 line) and the second in the form of a line close to an ellipse (2 line), passes through points $(a,0)$, $(b,0)$. For each value $p$ of the rotation angle of the first ellipse, the corresponding value of the rotation angle $q$ for the second ellipse is sought.
I find set of $p,q$ - but this set mixture values for 1 line and values for 2 line. How to sort this set in two sets for values for 1 line and values for 2 line ? This sets give correct animation this problem by Python.
# two equals ellipses (a*cos(t),b*sin(t)), (a*cos(s)+a+b,b*sin(s))
# contact rotation about centers (0,0) and (a+b,0)
# front and rear contact - two variants
# the first ellipse has a constant angular velocity of rotation
import scipy.optimize
import numpy as np
import math as m
import matplotlib.pyplot as plt
import time
fig = plt.figure()
plt.xlabel('X')
plt.ylabel('Y')
ax = fig.gca()
plt.axis('equal')
P=[]
def xx(p,t,sd):
return a*np.cos(t)*np.cos(p)+b*np.sin(t)*np.sin(p)+sd
def yy(p,t,sd):
return -a*np.cos(t)*np.sin(p)+b*np.sin(t)*np.cos(p)+sd
def f(y):
t,q,s = y[:3]
eq1=yy(p,t,0)-yy(q,s,0)
eq2=xx(p,t,0)-xx(q,s,a+b)
eq3=(a*m.sin(p)*m.sin(t) + b*m.cos(p)*m.cos(t))*(-a*m.sin(s)*m.cos(q) + b*m.sin(q)*m.cos(s))-(
(a*m.sin(q)*m.sin(s) + b*m.cos(q)*m.cos(s))*(-a*m.sin(t)*m.cos(p) + b*m.sin(p)*m.cos(t)))
return(eq1,eq2,eq3)
w, a, b, h, N =0, 3, 2, 2*m.pi, 36
for k in range (N):
p=h/N*k
x0 = np.array([3.3,0.3,-p+m.pi/6])
sol = scipy.optimize.root(f, x0, method='lm')
t,q,s=sol.x
e1,e2,e3=f(sol.x)
if abs(e1)+abs(e2)+abs(e3)<10**(-7):
w=w+1
q=m.fmod(q,2*m.pi)
#print (p,q)
P.append([p,q])
x, y =xx(p,t,0), yy(p,t,0)
plt.plot(x,y,'ro',markersize=1.)
ksi = np.arange(0, 2*np.pi, 0.01)
#plt.plot(np.sin(ksi)/2+(a+b)/2, 1*np.cos(ksi), lw=1)
for pq in P:
p, q =pq[0],pq[1]
plt.plot(xx(q,ksi,a+b), yy(q,ksi,0), lw=1)
plt.plot(xx(p,ksi,0), yy(p,ksi,0), lw=1)
plt.plot(0, 0, 'ro',markersize=2.)
plt.plot(a+b, 0, 'ro',markersize=2.)
print(w)
plt.grid()
plt.show()
#time.sleep(1)
#plt.plot(xx(q,ksi,a+b), yy(q,ksi,0), lw=0)
#plt.plot(xx(p,ksi,0), yy(p,ksi,0), lw=0)
Example animation
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
N=20
N1=120
a=3
b=2
t = np.arange(0, 2*np.pi+4*np.pi/N1, 2*np.pi/N1)
x = a*np.cos(t)
y = b*np.sin(t)
fig, ax = plt.subplots()
plt.axis('equal')
line1, = ax.plot(x, y, color = "r")
line2, = ax.plot(y+2*a+b, x, color = "g")
def update(t, x, y, line1, line2):
line2.set_data(x*np.cos(t/N)+y*np.sin(t/N), -x*np.sin(t/N)+y*np.cos(t/N))
line1.set_data(y*np.cos(-t/N)+(x)*np.sin(-t/N)+a+b , -y*np.sin(-t/N)+x*np.cos(-t/N))
return [line1,line2]
ani = animation.FuncAnimation(fig, update, int(2*np.pi*N), fargs=[ x, y, line1, line2],
interval=295, blit=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()
