Rotating disks

June 9, 2011
By

(This article was first published on Xi'an's Og » R, and kindly contributed to R-bloggers)

My neighbour is an half-retired entrepreneur who still runs his electric engine company. A few weekends ago, he came to me with the following physics question related with one of those engines: given a primary disk rotating at the angular speed of ω0 and a secondary disk located on the first one with a centre O1, a distance r0 between both centres, a radius of r1 and a relative angular speed of ω1, what is the absolute speed of a point M on the periphery of the secondary disk? Since I could not make sense of the solutions given in Wikipedia, I wrote a small (and crude) R code to show the location of M and to derive an approximation to the speed… Any suggestion for improvement is welcome!

#location of the second centre A
a=function(t,R=c(2,1),om=c(1,13)){
   sum(R)*c(cos(om[2]*t),sin(om[2]*t))}
#location of the peripheral point M
b=function(t,R=c(2,1),om=c(1,13)){
  a(t,R,om)+R[1]*c(cos(om[1]*t),sin(om[1]*t))}
#plot of the location of M as above
draw_position=function(R=c(1,2), om=c(3,1), phi=c(0,0) ){
position=function(t,index,R=c(1,2),om=c(3,1),phi=c(0,0)){
 relative_position=function(){
 return(list(R[index]*cos(om[index]*t+phi[index]),R[index]*sin(om[index]*t+phi[index])))
 }
 if(index == 1){
 return(relative_position() )
 } else {
 p1=position(t,index-1,R=R,om=om,phi=phi)
 p2=relative_position()
 return(list(p1[[1]]+p2[[1]],p1[[2]]+p2[[2]]))
 }}
 tes=seq(0,2*pi,length=10**3)
 xy_range=c(-1,1)*sum(abs(R))
 plot(0,0,pch="x",xlim=xy_range,ylim=xy_range,axes=F,xlab="",ylab="")
 pA=position(tes,1,R,om,phi)
 pB=position(tes,2,R,om,phi)
 lines(pB[[1]],pB[[2]],pch=19,cex=.4,col="tomato")
 lines(pA[[1]],pA[[2]],pch=19,cex=.2,col="blue")
 }
draw_position(om=c(1,30))
#angle at time t
the=function(t,R=c(2,1),om=c(1,13)){
  bb=b(t,R=R,om=om)
  bb=bb/sqrt(sum(bb^2))
  theta=acos(bb[1])
  if (bb[2]  theta}
#angular speed at time t
#by very crude differenciating
dthe=function(t,R=c(2,1),om=c(1,13)){
     dtes=mean(diff(tes))
    (the(t+dtes)-the(t-dtes))/(2*dtes)}
#new plot
plot(apply(as.matrix(tes),1,dthe,om=c(-55,2)),type="l",ylim=c(-2*pi,pi))

Antoine Dreyer actually contributed to improve the above code from an earlier version and he also derived the (Cartesian) components of the speed for me:

- \omega_0 r_0 \sin(\omega_0 t + \phi_0) - (\omega_0 + \omega_1)r_1 \sin(\omega_0 t + \phi_0 + \omega_1 t + \phi_1)

and

\omega_0 r_0 \cos(\omega_0 t + \phi_0) + (\omega_0 + \omega_1)r_1 \cos(\omega_0 t + \phi_0 + \omega_1 t + \phi_1)

if O1(0) has polar coordinates (r00) and M(0) has polar coordinates (r11) with respect to O1(0).

Filed under: R Tagged: angular speed, R

To leave a comment for the author, please follow the link and comment on their blog: Xi'an's Og » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: ,

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)