Elliptical Billiards: Triangular Orbits
Dan S. Reznik
April, 2019
Source files available here
Load packages
Shape of error function
Note: goes to zero as \(t\) approaches \(\pi/2\), because ray-trace path is shrinking.
draw_error <- function(a,t) {
p0 <- ell_p(a,t)
n0 <- ell_n(a,t)
df <- tibble(th=seq(-.99*pi/2,0,pi/500),
tr=map(th,tri_angle2,a=a,p=p0,norm=n0),
d2=map_dbl(tr,"d2")%>%sqrt,
d_adj=map_dbl(tr,"d_adj")) %>%
select(-tr)
df %>%
gather("key","value",d2,d_adj) %>%
ggplot(aes(th,value,group=key,color=key)) +
geom_line()
}
draw_error(2,pi/4)
Optimization functions
# tri_angle_d2 <- function(th,a,p,norm) tri_angle2(th,a,p,norm)$d2
tri_angle_d_adj <- function(th,a,p,norm) tri_angle2(th,a,p,norm)$d_adj
opt_th <- function(a,p,norm,tol=1e-4)
optimize(tri_angle_d_adj,
interval=c(-.99,-.05)*pi/2, # there's a direct back-and-forth hit at t in 8~9 degrees
tol=tol,a=a,p=p,norm=norm)
# global
#opt_th_global <- function(a,p,norm,tol=1e-4) optim(c(-.99*pi/2,0),function(x) #tri_angle_d_adj(x,a=a,p=p,norm=norm),method="BFGS")Calculate optimum \(\theta\) for this triangle
Investigate optimum \(\theta\) for various values of \(t\) from -90 to 0 degrees
df_opt_th <- tibble(t=seq(-90,90,.5),
opt_th=map(t%>%to_rad,
~opt_th(2,ell_p(2,.x),ell_n(2,.x)))) %>%
mutate(min=map_dbl(opt_th,"minimum"),
obj=map_dbl(opt_th,"objective")) %>%
select(-opt_th)
df_opt_th %>%
gather("key","value",min,obj) %>%
ggplot(aes(t,value,group=key,color=key)) +
geom_line()
Some values are causing minimizer to be unstable
summary(df_opt_th)
#> t min obj
#> Min. :-90 Min. :-1.0522 Min. :1.531e-07
#> 1st Qu.:-45 1st Qu.:-1.0177 1st Qu.:3.886e-05
#> Median : 0 Median :-0.8932 Median :1.012e-04
#> Mean : 0 Mean :-0.7813 Mean :1.158e-04
#> 3rd Qu.: 45 3rd Qu.:-0.5983 3rd Qu.:1.731e-04
#> Max. : 90 Max. :-0.1319 Max. :4.259e-04Values are between 8 and 9 degrees
Orbit perimeter is constant but area is not
df_per_area <- tibble(t=seq(0,359,5),
tris=t%>%map(get_opt_tri,a=2),
sides=tris%>%map(get_tri_sides),
per=sides%>%map_dbl(get_tri_per),
area=sides%>%map_dbl(get_tri_area))
df_per_area %>%
mutate(per) %>%
gather(key,value,per,area) %>%
ggplot(aes(t,value,group=key,color=key)) +
geom_line(size=1) +
scale_y_log10(breaks=scales::log_breaks(10))
Fixed point on boundary, vary output angle wrt normal
show_variable_angle <- function(a=2,t_deg=45,
th_min=-52,th_max=-50,th_step=.5) {
t <- t_deg %>% to_rad
p = ell_p(a,t)
n = ell_n(a,t)
n_tip = p + n*.1
arrow15 <- arrow(angle = 15, type = "closed",length=unit(10,"points"))
df_ell <- tibble(th = seq(th_min,th_max,th_step),
tri_angle = th %>% to_rad %>%
map(~tri_angle_df(.x,a,p,n))) %>%
unnest(tri_angle) # expands to the side!
# Animate
df_ell %>%
ggplot() +
# background
geom_point(x=0,y=0) +
geom_path(aes(px,py),data=ell_poly(a)) +
#geom_segment(x=p[1],y=p[2],
# xend=n_tip[1],yend=n_tip[2],
# color="#0000ff",arrow=arrow15) +
# moving
geom_path(aes(x=p_x,y=p_y,group=th)) +
geom_point(aes(x=p_x,y=p_y,group=th),color="blue") +
# labs(title = 'th: {th_range}') +
geom_point(x=p[1],y=p[2],color="red") +
coord_fixed() +
facet_wrap(~th) # transition_manual(row)
}
show_variable_angle()
Calculate incenter locus
Note: 150 degrees base vertex motion sweeps 360
get_incenter_locus <- function(a,deg_step=1)
seq(0,360,deg_step) %>%
map(~attr(get_opt_tri(.x,a),"incenter")) %>%
map_dfr(~tibble(x=.x[1],y=.x[2]))Compute / retrieve incenter locus as data frame
fname_incenter <- "data/incenter_locus.csv"
if (file.exists(fname_incenter)) {
df_incenter <- read_csv(fname_incenter)
} else {
df_incenter <- get_incenter_locus(2) %>%
mutate_all(~round(.,4))
# avoid multiple windings
max_row <- df_incenter %>%
{which(.$x<0&.$y>0&lead(.$y)<0)} # y zero crossing
df_incenter %>%
head(max_row) %>%
write_csv(fname_incenter)
}
#> Parsed with column specification:
#> cols(
#> x = col_double(),
#> y = col_double()
#> )Show orbit triangle
show_opt_tri <- function(t_deg,a,incenter_locus=NULL) {
df_ell_bound <- ell_poly(a)
df_ell_foci <- ell_foci(a)
df_opt_tri <- get_opt_tri(t_deg,a)
incenter <- attr(df_opt_tri,"incenter")
arrow15 <- arrow(angle = 15, type = "closed",length=unit(10,"points"))
ggplot(df_opt_tri) +
# fixed
geom_point(x=0,y=0,shape=4) +
geom_path(aes(px,py),df_ell_bound) +
#geom_point(x=p0[1],y=p0[2]) +
geom_point(aes(fx,fy),df_ell_foci) +
geom_polygon(aes(p_x,p_y),
linetype="dotted",
color="blue",fill=NA) +
geom_segment(aes(x=p_x,y=p_y,xend=ntip_x,yend=ntip_y),arrow=arrow15) +
geom_point(aes(p_x,p_y),color="blue") +
geom_point(aes(p_x,p_y),df_opt_tri%>%slice(1),color="red") +
{
if(is.null(incenter_locus))
geom_blank()
else
geom_path(aes(x,y)
,data=incenter_locus
,color="green"
,linetype=2
#,shape="."
)
} +
geom_point(x=incenter[1],y=incenter[2],
color="green",size=3) +
coord_fixed() +
labs(title="triangular orbit",
subtitle=glue("a={a}, t={t_deg}°"),x="",y="") #+
# theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
}
show_opt_tri(45,2,df_incenter)
Compute animation
Save frames to files
fnames_png <- list.files(path = "./pics",
pattern = "tri_.*\\.png$",
full.names = T)
if (length(fnames_png)>0)
file.remove(fnames_png)
tic()
plan("multiprocess")
deg_vec <- seq(0,360,1)
walk_vec <- deg_vec %>%
future_iwalk(~{
fname<-sprintf("./pics/tri_%04d.png",.y)
p_tri<-show_opt_tri(.x,2,df_incenter)
fname%>%ggsave(p_tri,dpi = 100)
# fname%>%image_read%>%image_trim%>%image_write(fname)
})
toc()Combine frames into animated .gif
fnames_png <- list.files(path = "./pics",
pattern = "tri_.*\\.png$",
full.names = T)
fnames_png %>%
map(image_read) %>% # reads each path file
map(image_trim) %>% # trims blank borders
image_join() %>% # joins image
# image_animate(fps=5) %>% # animates, can opt for number of loops
# image_write_video("tris_a=2.mp4",framerate=10)
image_write_gif("pics/tris_a=2.gif",delay=.1)
file.remove(fnames_png)