hr <- c(00,03,04,07,08,10,11,12,13,14,16,16,17,17,17,18,18,18,19,20,22)
min <- c(56,08,52,16,08,00,24,08,28,16,20,44,04,20,24,08,16,56,32,52,08)
eta <- 360*(hr+min/60)/24
eta <- pi*eta/180
n <- length(eta)
c <- sum(cos(eta))/n
s <- sum(sin(eta))/n
rho <- sqrt(c^2+s^2)
muhat <- 360+atan2(s,c)*180/pi
muhatdeg <- floor(muhat)
muhatmin <- round(60*(muhat-muhatdeg))
cat("Estimate mu as",muhatdeg,"degrees",muhatmin,"minutes",
    "and rho as",rho,"\n")