! motion of a billiard in stadium ! H. Nakanishi program billiards option nolet call initialize(r,alpha,x0,y0,vx0,vy0,nstep) call calculate(r,alpha,x0,y0,vx0,vy0,nstep) end ! initialize variables sub initialize(r,alpha,x0,y0,vx0,vy0,nstep) option angle degrees print "Billiards in a stadium; radius of circular part = 1" r = 1 v0 = 1 input prompt "elongation factor alpha -> ": alpha input prompt "initial x position -> ": x0 input prompt "initial y position -> ": y0 input prompt "initial angle (degrees) -> ": theta0 vx0 = v0*cos(theta0) vy0 = v0*sin(theta0) input prompt "how many bounces? -> ": nstep aspect = 1.333 set window -(1+alpha)*r*aspect,(1+alpha)*r*aspect,-(1+alpha)*r,(1+alpha)*r set background color "white" set color "black" clear plot -(1+alpha)*r*aspect,0;(1+alpha)*r*aspect,0 plot 0,(1+alpha)*r;0,-(1+alpha)*r set color "green" plot -alpha*r,r;alpha*r,r; x = alpha*r for i=1 to 100 x = x+r/100 dx = i*r/100 plot x,sqr(r^2-dx^2); next i for i=1 to 100 x = x-r/100 dx = r-i*r/100 plot x,-sqr(r^2-dx^2); next i plot alpha*r,-r;-alpha*r,-r; x = -alpha*r for i=1 to 100 x = x-r/100 dx = i*r/100 plot x,-sqr(r^2-dx^2); next i for i=1 to 100 x = x+r/100 dx = r-i*r/100 plot x,sqr(r^2-dx^2); next i plot set color "green" plot area: -alpha*r-0.01,0.01;-alpha*r+0.01,0.01;-alpha*r+0.01,-0.01;-alpha*r-0.01,0-0.01 plot area: alpha*r-0.01,0.01;alpha*r+0.01,0.01;alpha*r+0.01,-0.01;alpha*r-0.01,0-0.01 set color "red" plot area: x0-0.01,y0+0.01;x0+0.01,y0+0.01;x0+0.01,y0-0.01;x0-0.01,y0-0.01 end sub ! use exact method sub calculate(r,alpha,x0,y0,vx0,vy0,nstep) set color "blue" plot x0,y0; vx1 = vx0 vy1 = vy0 x1 = x0 y1 = y0 for i=1 to nstep if vy1>0 then x = (r-y1)*vx1/vy1+x1 if -alpha*r <= x and x <= alpha*r then y = r vx = vx1 vy = -vy1 else if x>alpha*r then a = 1+(vx1/vy1)^2 b = 2*vx1/vy1*(x1-alpha*r-vx1/vy1*y1) c = (x1-alpha*r-vx1/vy1*y1)^2-r^2 y = (-b+sqr(b^2-4*a*c))/(2*a) x = sqr(r^2-y^2)+alpha*r nx = (x-alpha*r)/sqr((x-alpha*r)^2+y^2) ny = y/sqr((x-alpha*r)^2+y^2) vperpx = (vx1*nx+vy1*ny)*nx vperpy = (vx1*nx+vy1*ny)*ny vparx = vx1-vperpx vpary = vy1-vperpy vx = vparx-vperpx vy = vpary-vperpy else a = 1+(vx1/vy1)^2 b = 2*vx1/vy1*(x1+alpha*r-vx1/vy1*y1) c = (x1+alpha*r-vx1/vy1*y1)^2-r^2 y = (-b+sqr(b^2-4*a*c))/(2*a) x = -sqr(r^2-y^2)-alpha*r nx = (x+alpha*r)/sqr((x+alpha*r)^2+y^2) ny = y/sqr((x+alpha*r)^2+y^2) vperpx = (vx1*nx+vy1*ny)*nx vperpy = (vx1*nx+vy1*ny)*ny vparx = vx1-vperpx vpary = vy1-vperpy vx = vparx-vperpx vy = vpary-vperpy end if else if vy1<0 then x = -(r+y1)*vx1/vy1+x1 if -alpha*r <= x and x <= alpha*r then y = -r vx = vx1 vy = -vy1 else if x>alpha*r then a = 1+(vx1/vy1)^2 b = 2*vx1/vy1*(x1-alpha*r-vx1/vy1*y1) c = (x1-alpha*r-vx1/vy1*y1)^2-r^2 y = (-b-sqr(b^2-4*a*c))/(2*a) x = sqr(r^2-y^2)+alpha*r nx = (x-alpha*r)/sqr((x-alpha*r)^2+y^2) ny = y/sqr((x-alpha*r)^2+y^2) vperpx = (vx1*nx+vy1*ny)*nx vperpy = (vx1*nx+vy1*ny)*ny vparx = vx1-vperpx vpary = vy1-vperpy vx = vparx-vperpx vy = vpary-vperpy else a = 1+(vx1/vy1)^2 b = 2*vx1/vy1*(x1+alpha*r-vx1/vy1*y1) c = (x1+alpha*r-vx1/vy1*y1)^2-r^2 y = (-b-sqr(b^2-4*a*c))/(2*a) x = -sqr(r^2-y^2)-alpha*r nx = (x+alpha*r)/sqr((x+alpha*r)^2+y^2) ny = y/sqr((x+alpha*r)^2+y^2) vperpx = (vx1*nx+vy1*ny)*nx vperpy = (vx1*nx+vy1*ny)*ny vparx = vx1-vperpx vpary = vy1-vperpy vx = vparx-vperpx vy = vpary-vperpy end if else if vx1>0 then y = y1 x = alpha*r+sqr(r^2-y1^2) nx = (x-alpha*r)/sqr((x-alpha*r)^2+y^2) ny = y/sqr((x-alpha*r)^2+y^2) vperpx = (vx1*nx+vy1*ny)*nx vperpy = (vx1*nx+vy1*ny)*ny vparx = vx1-vperpx vpary = vy1-vperpy vx = vparx-vperpx vy = vpary-vperpy else if vx1<0 then y = y1 x = -alpha*r-sqr(r^2-y1^2) nx = (x+alpha*r)/sqr((x+alpha*r)^2+y^2) ny = y/sqr((x+alpha*r)^2+y^2) vperpx = (vx1*nx+vy1*ny)*nx vperpy = (vx1*nx+vy1*ny)*ny vparx = vx1-vperpx vpary = vy1-vperpy vx = vparx-vperpx vy = vpary-vperpy else set color "red" print "billiard is stopped!" exit sub end if end if x1 = x y1 = y vx1 = vx vy1 = vy plot x,y; next i plot end sub