Race.m

Needs["RungeKutta`",":NumericalMath:RungeKutta.m"]

BeginPackage["Race`","RungeKutta`"]

Usage Messages

Race::usage = 
			"Race[curvelist, (options)] simulates beads sliding down 
wire frames given by the curves in curvelist(limited to 5).
These curves should start at the origin with the orientation
that the positive y-axis is down. If more than one curve is to be 
raced, they must be delimited by set brackets. A curve may 
be entered in the form {y,{x,xmin,xmax},options} where y is 
an expression in the variable x or as {x,y,{t,tmin,tmax},options} 
where x and y are expressions in the parameter t. Useful 
options are: Animation -> True, the default is False which 
returns a plot of the curve, length and time; FrictionCoeff -> number,
the default is zero; and Name -> curve label. The notebook
RaceImplementation has a full account of all detatils and options; 
see Message[Race::notes] for a brief trouble shooting message.";

Race::notes =
			"The curve must be physically meaningful, i.e. the initial
slope must be strictly greater than m, where m is the coefficient
of friction and the ending point must satisfy ycoord >= m xcoord
(the default value for µ is 0). Use exact values such as 0, Pi, Sqrt[2],
ArcSin[1/2],etc. for the range values so Mathematica can recognize
that the integral giving the time for the bead to slide down the curve
is infinite at the starting point (or ending point when this is to be
on the line y=mx) and not just a very large number; otherwise, invalid
results may occur. Mathematica attempts to carry out a series expansion
symbolically for the curve, and even when exact values are used, it may
get hung up with more complicated curves. This is especially true when
friction is included and attempting to end the curve on the line of
repose y=mx. It may be necessary to move the ending point down slightly
from this line (or see SeriesExpansion option in the notebook
RaceImplementation for an alternative).";

CycloidEqns::usage = "CycloidEqns[a,b], for a>0 and b>=0, finds the 
constants r and tfinal in the parameterization, 
x=r(t-Sin[t]), y=r(1-Cos[t]), 0<=t<=tfinal, that gives the cycloid 
starting at the origin with a vertical tangent and ending at the 
point (a,b). It returns {r,tfinal}.";

Options

Options[Race] = {Animation -> False, Name -> Automatic,
				 AccuracyCheck -> False, Dimension -> Automatic,
				 AnimationFrames -> Automatic, FrictionCoeff -> 0,
				 FreeFall -> Automatic, SlopeTolerance -> 150,
				 ImplicitStepSize -> Automatic, ReposeLine -> False,
				 SeriesExpansion -> True, DiffExpr -> Automatic,
				 ComparisonCurve -> "Cycloid"};

Options[Graph] = {PlotStyle -> Automatic, Prolog -> {},
				  Epilog -> {}, PlotPoints -> 25,
				  PlotColor -> Automatic, PlotDivision -> 20.,
				  MaxBend -> 10.};

Begin["`Private`"]

Internal validity checking messages

Race::complex = 
			"The problem of finding the time for a bead to slide down a curve 
involves an equation with an outermost square root. Mathematica 
automatically applies the simplification rule, Sqrt[-x]->Sqrt[x]I, where 
I is the imaginary unit; this is wrong for x<0! This rule will be applied
to any curve that has an outermost minus sign and will lead to undesired 
complex output. This minus sign can usually be hidden by pulling it inside
some other part of the function defining your curve; e.g. y=-1/2 x(x-5) 
could be written as y=1/2 x(5-x). Also, remember that the positive y-axis 
is oriented downward. If the time printed in the Race::time message 
appears to be real and correct except off in sign, then the program will 
probably correct this. Any other errors are probably due to an invalid 
curve or a failure to numerically handle this example.";

Race::time = 
			"Complex output in evaluation of the time: `1` 
for the curve, `2`,`3` If the time appears to be real and correct or just 
off in sign, this will probably be	corrected (See Message[Race::complex]
for more information). If not, you may want to abort an animation by
choosing interrupt in the Action menu since this indicates either an
invalid curve or a failure to numerically handle this example.";

Race::initptfail = 
			"The starting point, x = `1` and y = `2` for the curve, `3`, `4` 
is not at the origin. If this point is approx. at the origin, then this 
message indicates that not enough accuracy has been retained and there 
would likely be problems in obtaining valid numerical results. The 
problem is sensitive with respect to this starting point (see 
Message[Race::notes]).";

Race::initptwarning = 
			"The range variable xmin = `1` for the curve, `2`, `3` is not an
exact value. Since finding the time in this problem can be sensitive with
respect to the initial point, this may give invalid results and/or
drastically slow down the computations. If this happens, try using exact
values for the range variables (e.g. Pi,Sqrt[2],ArcSin[1/2],etc.see
Message[Race::notes]).";

Race::initslopefail = 
			"Invalid initial slope, `1` for the curve, `2`, `3` 
for the bead to begin sliding. If friction is ignored, the initial 
slope must be positive (remember the positive y-axis is to be oriented 
downward!). In general, this slope must be strictly greater than the 
friction coefficient.";

Race::endptfail = 
			"The ending point x = `1` and y = `2` for the curve, `3`, `4` 
does not satisfy, yend >= m xend, where m is the coefficient of
friction (if it is approx.true then this message indicates that not
enough accuracy has been retained in the calculation of the final
point - use exact range values or keep higher accuracy with computed
values).";

Race::length =
			"Mathematica's message indicates either a failure to guarantee
accuracy or to numerically compute the length of curve, `1`, `2` 
This will not affect the results for the time or an animation of the
curve and can usually be corrected by increasing the working precision,
if an infinite result occured (but will greatly slow down computations),
or increase MaxRecursion if an accuracy failure (see NIntegrate
options under SeriesExpansion option in RaceImplementation notebook).";

Race::endptwarning = 
			"The end point x = `1` and y = `2` for the curve, `3`, `4` 
is approx. on the line of repose y = m x, where m is the coefficient of 
friction. But, the range variable xmin = `5` 
for this curve is not an exact value and finding the time in this case
can be sensitive with respect to this final point; this may give invalid
results and drastically slow down the computations. If this happens, try
using exact values for the range variables (e.g. Pi,Sqrt[2],ArcSin[1/2],
etc.see Message[Race::notes]).";

Race::timefail =
			"Unable to determine correct time for the curve, `1`, `2` due to
complex output. Try setting SeriesExpansion->False (see
RaceImplementation concerning this option) or moving the end point off
the line y=µx slightly, where µ is the coefficient of friction.";

Race::initptinfy =
			"Failure to numerically compute time for the curve, `1`, `2`  Try
increasing WorkingPrecision or decreasing AccuracyGoal (in
RaceImplementation notebook under Numerical Options for computing the
time, see SeriesExpansion option concerning NIntegrate settings for the
initial point; currently WorkingPrecision is at default value,
AccuracyGoal is 4 and MaxRecursion is 10).";

Race::endptinfy =
			 "Failure to numerically compute time for the curve, `1`, `2`  Try
increasing WorkingPrecision or decreasing AccuracyGoal (in
RaceImplementation notebook under Numerical Options for computing the
time, see SeriesExpansion option concerning NIntegrate settings for the
end point; currently WorkingPrecision is 30, AccuracyGoal is 4 and
MaxRecursion is 10)";

CycloidEqns::Endpt = 
	"The end point a = `1` and b = `2` fails to satisfy a>0 and b>=0.";

Derivative rule for the absolute value function so it can be used symbolically to help prevent accidental complex output.

protected=Unprotect[Abs,Sign]
	Abs/: Abs':= Sign
	Sign[0]=Indeterminate
Protect[Release[protected]];

FilterOptions[command_Symbol, opts___]:=
	Block[{keywords=First /@ Options[command] },
		Sequence @@ Select[{opts},MemberQ[keywords,First[#]]& ]
	]

Single function form for Race

Race[f_, {x_,xmin_,xmax_},options___Rule]:=
		Race[{x,f,{x,xmin,xmax},options},options] /; !SameQ[Head[f],List] && NumberQ[N[xmin]] && NumberQ[N[xmax]]

Race[x_, y_, {t_,tmin_,tmax_},options___Rule]:=
		Race[{x,y,{t,tmin,tmax},options},options] /; !SameQ[Head[x],List] && !SameQ[Head[y],List] && NumberQ[N[tmin]] && NumberQ[N[tmax]]

Main definition Block for Race

Race[curves__List,options___Rule]:=

Block[{ len=Length[{curves}],curvelist={curves},numdata,graphdata,
								tnumdata,tgraphdata,g=32.174,r,tend,comparisons,comptime,
								percents,place,
								compcurve=ComparisonCurve/.{options}/.Options[Race],
					   colornames={"blue","red","magenta","green",
					  										   "brown"},
					   colors={RGBColor[0,0,.6],RGBColor[1,0,0],
					           RGBColor[0.967, 0.048, 0.962],RGBColor[.14,.56,.14],
					           RGBColor[.8,.5,.2]},
					   colorrules,curvenames,cases,animationlist={},xdim,ydim,ymax,
					   dim,xend,yend,labeldim,i,j,n,slowest,diff,centerlist,timestep,
					   frames=AnimationFrames/.{options}/.Options[Race],
					   frametext,time,time1,time2,point,legend,positions 
					  },
		
		colorrules=(Color->#)& /@ Take[colors,len];
		
		numdata=Check[numerical /@ curvelist,Return[],Race::initptfail,
												Race::initslopefail,Race::endptfail,Race::timefail,
												Race::initptinfy,Race::endptinfy];
		tnumdata=Transpose[numdata];
	 
	 graphdata=Graph /@ (curvelist=Release[Thread[
	 																					  Hold[Append][curvelist,colorrules] ]]);
	 tgraphdata=Transpose[graphdata];

			 
	 (*   LABELING
	   First, set up dimensions for lengths of axes to accomodate the
	   number of curves to appear in the legends *)
	 
	 {xend,yend}=numdata[[1,1]];
	 xdim=1.1xend;
	 ymax=Max[ tnumdata[[2]] ];
	 dim=Max[xend,ymax];
	 ydim=1.1ymax + .07(len+1)dim;
	 labeldim=.07len dim;
	 comparisons=Table["Slower",{len}];
	 	 
	 (*Determine curve names for labeling by replacing Automatic (the
	   default name) with "Curve#", where # is the position of Automatic 
	   in the curvenames list*)
	 
	 curvenames=tgraphdata[[1]];
	 cases=Cases[curvenames,_String];
	 (cases=Insert[cases,SequenceForm["curve",#],#])& /@ Flatten[
	 														Position[curvenames,Automatic] ];
	 curvenames=cases;
	 
	 legend = 
	 	ColumnForm[ Thread[
	 				SequenceForm[" ",#1,"-",#2,", length= ",#3]& [
	 							curvenames,Take[colornames,len],N[tnumdata[[3]],4] ] 
	 	]];
	 	
		point=SequenceForm["(",N[xend,3],",",N[yend,3],")"];
		
		time1=SequenceForm["Times(comparison w/",compcurve,")"];
		
		Which[SameQ[compcurve,"Cycloid"],
										{r,tend}=CycloidEqns[xend,yend];
										comptime=tend Sqrt[r/g],
								SameQ[compcurve,"Line"],comptime=Sqrt[2(xend^2+yend^2)/(g yend)],
								True,comptime=Min[ tnumdata[[4]] ]
							];
		
		percents=(100(#-comptime)/comptime)& /@ tnumdata[[4]];
				
		If[SameQ[compcurve,"FastestRaced"],
							place=Flatten[Position[ tnumdata[[4]],comptime ]][[1]];
							comparisons[[place]]="Fastest";
							percents[[place]]=0
				];
		If[SameQ[compcurve,"Line"],
					If[Negative[ percents[[#]] ],comparisons[[#]]="Faster"]& /@ Range[len]
				];

		time2=ColumnForm[ Thread[
									SequenceForm["  ",#1,":",#2,"(",#3,"% ",#4,")"]& [
											curvenames,N[tnumdata[[4]],6],Chop[N[Abs[percents],3]],comparisons ]
								]];
		time=ColumnForm[{time1,time2}];
		frametext[j_]:=ColumnForm[{ SequenceForm["Frame"],
																														SequenceForm[" #",j] }];
		
		(* Determine which curves require animation *)
			
		AppendTo[animationlist,curvelist[[#]] ]& /@ (positions=Flatten[ 
																													Position[Transpose[graphdata][[2]],True] ]); 
	 
		If[Length[animationlist]==0,
				(* THEN no animation required; show a static plot of the curves 
				   indicating numerical data *)
				Show[Transpose[graphdata][[3]], Graphics[{ Point[{xend,-yend}],
					Text[legend,{0,.02dim},{-1,-1}],
					Text["y-axis",{0,-ydim},{-1,1}],
					Text[point,{xend + .02dim,-(yend+.02dim)},{-1,1}],
					Text[time,{.2dim,-ydim+.02dim},{-1,-1}] }],
					FilterOptions[Graphics,options],
					PlotRange -> {{-.1dim,xdim}, {-ydim,labeldim}},
					Ticks -> None, AxesLabel -> {"x-axis         ",""},
					AspectRatio -> (ydim+labeldim)/(xdim+.1dim),
					DisplayFunction -> $DisplayFunction]; Return[] 
		];
		
		(* ELSE pass animation list to 'centers' for RungeKutta routine *)
		
		(*1st determine the number of graphic frames to animate and based on 
		  the slowest curve  set a fixed timestep for finding
		  the centers of the beads (Note: The RungeKutta routine uses 
		  adaptive step size control, so this timestep doesn't limit accuracy*)
		
		slowest=Max[ Transpose[numdata][[4]] ];
		If[SameQ[frames,Automatic],
				(* THEN - this algorithm sets frames=7 for times <.3sec. incrementing
				   1 frame every .1 sec.*)
				If[(diff=slowest -.2)>0,frames=7+Floor[diff/.1],frames=7 ]
		];
				
		timestep=slowest/frames;
	  
		centerlist=centers /@ Release[ Thread[ Hold[Prepend] [animationlist,
																Thread[List[Transpose[numdata][[4,#]]& /@ positions,
																																									frames,timestep]] ]
																									]];
		
		Do[
				If[n<frames+1,
							(* THEN - display graphs and beads but not final times *)
							Show[Transpose[graphdata][[3]],
												Graphics[{ PointSize[.045],
														Table[{centerlist[[i,1]],Point[centerlist[[i,2,n]] ]},
																{i,Length[animationlist]}] }],
												Graphics[{ Point[{xend,-yend}],
														Text[legend,{0,.02dim},{-1,-1}],
														Text["y-axis",{0,-ydim},{-1,1}],
														Text[point,{xend + .02dim,-(yend+.02dim)},{-1,1}],
														Text[frametext[n],{-.2dim,.02dim},{-1,-1}] }],
												FilterOptions[Graphics,options],
												PlotRange -> {{-.2dim,xdim}, {-ydim,labeldim}},
												Ticks -> None, AxesLabel -> {"x-axis         "," "},
												AspectRatio -> (ydim+labeldim)/(xdim+.2dim),
												DisplayFunction -> $DisplayFunction],
								(* ELSE - for the last frame, include the curve times legend *)
								Show[Transpose[graphdata][[3]],
												Graphics[{ PointSize[.045],
														Table[{centerlist[[i,1]],Point[centerlist[[i,2,n]] ]},
																{i,Length[animationlist]}] }],
												Graphics[{ Text[legend,{0,.02dim},{-1,-1}],
														Text["y-axis",{0,-ydim},{-1,1}],
														Text[point,{xend + .02dim,-(yend+.02dim)},{-1,1}],
														Text[frametext[n],{-.2dim,.02dim},{-1,-1}],
														Text[time,{.2dim,-ydim+.02dim},{-1,-1}] }],
												FilterOptions[Graphics,options],
												PlotRange -> {{-.2dim,xdim}, {-ydim,labeldim}},
												Ticks -> None, AxesLabel -> {"x-axis         "," "},
												AspectRatio -> (ydim+labeldim)/(xdim+.2dim),
												DisplayFunction -> $DisplayFunction]
						],
				{n,frames+1}];
				
				Print[SequenceForm[
											"Animation Setting: ",N[1/timestep,4],"frames/sec."]]
		
(* end block for race *) ]

Subsidiary function, numerical, for calculation of times and lengths of curves

numerical[{f_,{x_,xmin_,xmax_},opts___Rule}]:=
		numerical[{x,f,{x,xmin,xmax},opts}]

numerical[{x_,y_,{theta_,thetamin_,thetamax_},opts___Rule}]:=
			
Block[{px=Function[theta,Release[x]], py=Function[theta,Release[y]],
							theta0=N[thetamin],thetaf=N[thetamax],dpx,dpy,numerator,integrand,
							dim=N[Dimension /.{opts}/.Options[Race] ],g=N[32174/1000],
		   		m=FrictionCoeff/.{opts}/.Options[Race],initslope,length,max,
		   		xend,yend,endtest,repose=False,time,time1,time2,time3=0,
		   		series=SeriesExpansion/.{opts}/.Options[Race],integral,
		   		diffexpr=DiffExpr/.{opts}/.Options[Race]},
		   		
		 (* Check if starting point is approx. at the origin *)
		 
		 If[(Chop[px[thetamin]//N]!=0) || (Chop[py[thetamin]//N]!=0),
	  				Message[ Race::initptfail,px[thetamin],py[thetamin],x,y ];
	  				Return[]
	  		];
	  If[(!SameQ[Precision[thetamin],Infinity]) && series,
	  					Message[Race::initptwarning,thetamin,x,y]
	  		];
	  												
		 dpx=px'; dpy=py';
		 
		 (* Check if initial slope > m *)
		 
		 If[NumberQ[
	  				(initslope=(Limit[dpy[theta]/dpx[theta],theta->thetamin]//N))
	  										] && (initslope <= N[m]),
	  			Message[Race::initslopefail,initslope,x,y];Return[]
	  ];
	  
	  (* Check that final point satisfies y >= m x *)
	  
	  {xend,yend}=N[{px[thetaf],py[thetaf]}];
	  If[((endtest=N[yend- m xend]) < 0) && (Chop[endtest]!=0),
	  		 Message[Race::endptfail,xend,yend,x,y]; Return[]
	  		];
		 
		 numerator=Sqrt[dpx[theta]^2 + dpy[theta]^2];
		 If[SameQ[diffexpr,Automatic],
		 			integrand=numerator/Sqrt[2g(py[theta]-m px[theta])],
		 		 (* ELSE,integrand is the reciprocal of this expr. but remove
		 		 			Absolute values if necessary *)
		 		 If[ FreeQ[diffexpr,Abs[_]],
		 		 				integrand=1/diffexpr,
		 		 				integrand=1/Apply[Sequence,diffexpr,{1}]
		 		 ]
		 ];
		
			Check[
					(length=NIntegrate[ numerator,{theta,theta0,thetaf},
																AccuracyGoal->2,
																Release[ FilterOptions[NIntegrate,opts] ],
																MaxRecursion->10 ]),
					(Message[Race::length,x,y];length="length"),
					 Power::infy,NIntegrate::notnum,NIntegrate::conv ];
		
			(* Check if the ending point is on the line of repose - which 
			   would give another singularity to avoid in integration of time *)
			
			If[Chop[endtest]==0,
					 repose=True;
					 If[(!SameQ[Precision[thetamax],Infinity]) && series,
					 			Message[Race::endptwarning,xend,yend,x,y,thetamax]
					 ]
			];
			
			If[series,
			   (* THEN, use series expansions around singularities *)
					 time1=N[((integral=Integrate[
																										 Normal[Series[integrand,{theta,thetamin,8}]],
										 														 theta]) /. theta->thetamin + 1/1000) -
										 			  (integral /. theta->thetamin) ];
			
					 If[repose,
						   (*THEN*)
								 time2 = NIntegrate[integrand,{theta,theta0+.001,thetaf-.001}];
								 time3 = N[ ((integral=Integrate[
																											Normal[Series[integrand,{theta,thetamax,8}]],
										 																theta]) /. theta->thetamax) -
										 										(integral /. theta->thetamax-1/1000) ];
							  If[Chop[Im[time3]]==0,
							 			 time3=Abs[Re[time3]],
							 			 If[Chop[Re[time3]]==0,
							 					  time3=Abs[Im[time3]],
							 					  (* ELSE both Re and Im parts are non-zero and we can't 
							 					     determine the correct time! *)
							 					  Message[Race::timefail,x,y];Return[]
							 			 ]
							 	],
							 		
							 (*ELSE*)
								time2 = NIntegrate[integrand,{theta,theta0 + .001,thetaf}]
						];
						
						If[SameQ[ Head[ (time=time1+time2+time3) ],Complex ],
						  (*THEN take real part and print warning message*)
						   Message[Race::time,N[time,19],x,y ];
						   time=Abs[Re[time1+time2]]+time3
				  ],
						
				(* ELSE, series is false so use NIntegrate over the whole interval *)
						Check[
									(time1=NIntegrate[integrand,{theta,theta0,theta0+.001},
																		Release[ FilterOptions[NIntegrate,opts] ],
																		AccuracyGoal->4,MaxRecursion->10 ]),
									(Message[Race::initptinfy,x,y];Return[]),
									Power::infy,NIntegrate::notnum ];
						If[repose,
								 (* THEN *)
								 time2=NIntegrate[integrand,{theta,theta0+.001,thetaf-.001}];
								 Check[
											(time3=NIntegrate[integrand,{theta,thetaf-.001,thetaf},
																				Release[ FilterOptions[NIntegrate,opts] ],
																				AccuracyGoal->4,WorkingPrecision->30,
																				MaxRecursion->10]),
											(Message[Race::endptinfy,x,y];Return[]),
												Power::infy,NIntegrate::notnum ],
									(* ELSE *)
									time2=NIntegrate[integrand,{theta,theta0+.001,thetaf},
																	Release[ FilterOptions[NIntegrate,opts] ] ]
						];
						
						If[SameQ[ Head[ (time=time1+time2+time3) ],Complex ],
						  (*THEN take real part and print warning message*)
						   Message[Race::time,N[time,19],x,y ];
						   time=Abs[Re[time1+time2+time3]]
				  ]
					
			];
					

(* calculate the max. y-coordinate for the curve to determine relative 
				placement of labels and suppress messages to console*)
		
		If[SameQ[dim,Automatic],
			
				(* THEN *)	
				Off[FindMinimum::brkfail];Off[FindMinimum::regext];Off[FindMinimum::convFM];
				If[(* FindMinimum returns a real number *)
							SameQ[
									Head[
										(max=First[FindMinimum[-py[theta],{theta,(theta0+thetaf)/2,theta0,thetaf}]])
													], Real],
								(* Then *)
								dim=Max[yend,Abs[max]],
								(* ELSE *)
								dim=yend
						];
					On[FindMinimum::brkfail];On[FindMinimum::regext];On[FindMinimum::convFM]
		];
		
		Return[{{xend,yend},dim,length,time}]
	
	(* end of Block for numerical *) ] /; ( NumberQ[N[thetamin]] && NumberQ[N[thetamax]] )

Subsidiary function, Graph, for a plotting the curves

Graph[{f_,{x_,xmin_,xmax_},opts___}]:=
	Graph[{x,f,{x,xmin,xmax},opts}]

Graph[{x_,y_,{theta_,thetamin_,thetamax_},opts___}]:=
	
Block[{px=Function[theta,Release[x]], py=Function[theta,Release[y]],
		color=Color/.{opts},animation=Animation/.{opts}/.Options[Race],
		name=Name/.{opts}/.Options[Race],xend,m,label,repose={},
		reposeline=ReposeLine/.{opts}/.Options[Race]},
		
		If[reposeline,
				 xend=px[thetamax]//N;
				 m=FrictionCoeff/.{opts}/.Options[Race];
				 label=SequenceForm["y=",m,"x"];
				 repose={RGBColor[0,0,0],Dashing[{.01,.01}],
				 								Line[{ {0,0},{1.1xend,-1.1m xend} }],
				 								Text[label,{1.1xend/2,-1.1m xend/2},{-1,-1}]}
			];
		
		(* generate graphs but not displayed *)
		
		curveplot=ParametricPlot[{px[theta],-py[theta]},
															{theta,thetamin,thetamax},
															Release[ FilterOptions[Graph,opts] ],
															PlotStyle->{ {color,repose} },
															DisplayFunction -> Identity];
		
		Return[{name,animation,curveplot}]
		
(* end of block for Graph *) ]

Subsidiary function, centers, for simulation of a bead sliding down the curves

centers[{{time_,frames_,timestep_},f_,{x_,xmin_,xmax_},opts___}]:=
		cntr[time,frames,timestep,x,f,{x,xmin,xmax},opts]

centers[{{time_,frames_,timestep_},x_,y_,{theta_,thetamin_,thetamax_},
									opts___}]:=
		cntr[time,frames,timestep,x,y,{theta,thetamin,thetamax},opts]

cntr[curvetime_,frames_,timestep_,x_,y_,
										{theta_,thetamin_,thetamax_},opts___]:=

Block[{px=Function[theta,Release[x]],	py=Function[theta,Release[y]],
							theta0=N[thetamin],thetaf=N[thetamax],g=N[32174/1000],
							clist,color=Color/.{opts},dpx,dpy,diffexpr,
							m =FrictionCoeff/.{opts}/.Options[Race],
							slope=N [SlopeTolerance/.{opts}/.Options[Race] ],
							stepsize=N[ ImplicitStepSize/.{opts}/.Options[Race] ],
							ff=FreeFall/.{opts}/.Options[Race],impltheta,thetadrop,guess,
							theta1,implicit,freefall,tvert,vertsteps,rkstep,
							constant=Ceiling[Chop[frames-curvetime/timestep,10^-15]],i,k,
							diffexpr=DiffExpr/.{opts}/.Options[Race]},
						 
		dpx=px'; dpy=py';
		If[SameQ[diffexpr,Automatic],
					diffexpr=Sqrt[Abs[
							(2g(py[theta]-m  px[theta]))/(dpx[theta]^2 + dpy[theta]^2) ]]
		];
	
		(* START-UP ROUTINES: IMPLICIT and FREE-FALL *)
		
		implicit:=(
	   		
	   		If[SameQ[stepsize,Automatic],stepsize=1/1000];
	   		guess=stepsize/(10curvetime)(thetaf-theta0) + theta0;
	     impltheta=1/6(4 theta + (stepsize)diffexpr) /.
	   										FindRoot[theta==stepsize(diffexpr/.theta->theta/4),
	   														{theta,guess},AccuracyGoal->12];
	   			
	   			theta1=Last[
	   					RungeKutta[{N[diffexpr]},{t,theta},{stepsize,impltheta},
	   											timestep,10^-7] ][[2]];
	   				
	   			clist={{0,0}};
	   			AppendTo[ clist,{px[theta1]//N,-py[theta1]//N} ];
	   			rkstep=1
	   );
				
			freefall:=(
	   	
	   		thetadrop=theta/.FindRoot[dpy[theta]/dpx[theta]==slope,
	   														{theta,(thetaf-theta0)/10^6},AccuracyGoal->12,
	   															MaxIterations->20];
	   		
	   		tvert=Sqrt[2py[thetadrop]/g]//N;
	   		vertsteps=Floor[tvert/timestep];
	   		clist=Table[{0,-g/2(k timestep)^2},{k,0,vertsteps}];
	   						
	   		theta1=Last[
	   				RungeKutta[{N[diffexpr]},{t,theta},{tvert,thetadrop},
	   					(vertsteps+1)timestep,10^-7] ][[2]];
	   		AppendTo[ clist,{px[theta1]//N,-py[theta1]//N} ];
	   		rkstep=vertsteps+1
	  );	   	
	   	
	 (* Determine start-up routine based on - free-fall if a vertical 
	 		 initial tangent otherwise use an implicit step *)
	 		   		
		If[SameQ[ff,Automatic], ( Off[Limit::nlm];
	   	  		If[NumberQ[Limit[dpy[theta]/dpx[theta],theta->theta0]//N],
	   	  					ff=False,ff=True];
	   	  		On[Limit::nlm] )
	   	 ];
	   	 
	 If[ff,freefall,implicit];

		
		(* Do a loop of 1 step rungekutta runs to handle the adaptive stepsize
		   control so the beads are animated at equal time steps *)
		
		
		Do[
						theta1=Last[
														RungeKutta[{N[diffexpr]},{t,theta},{i timestep,theta1},
																		(i+1)timestep,10^-7] ][[2]];
						AppendTo[clist,{px[theta1]//N,-py[theta1]//N}],
		{i,rkstep,frames-constant-1}];
		
		If[AccuracyCheck/.{opts}/.Options[Race],
					Print[SequenceForm["Computed end pt. ",N[{px[theta1],py[theta1]}],
									", Error=",Abs[{px[theta1],py[theta1]}-{px[thetaf],py[thetaf]}]//N]]
		];	
		
		(* Now, Append the constant centers at the end point *)
		clist=clist~Join~Table[{px[thetaf],-py[thetaf]},{constant}];
		
		Return[{color,clist}]

(* end of block for centers *) ]

Subsidiary function, CycloidEqns, for finding the parameterization of a cycloid as a comparison curve

CycloidEqns[a_,b_]:=
	Block[{guess,r,tfinal,t},
	
		If[( !Positive[N[a]] || Negative[N[b]] ),
			Message[CycloidEqns::Endpt,a,b];Return[]
		];
	
		guess := 2.2(a/b) /; (N[b] > 0 && N[a] <= N[b Pi/2]);
		guess := Pi(2-(Pi b)/(2a)) /; N[b/a] < N[2/Pi]; 
	
		tfinal = t /. FindRoot[ a(1 - Cos[t]) == b(t - Sin[t]),
						{t, guess}];
	
		r = a/(tfinal - Sin[tfinal]);
	
		{r,tfinal}
	
	] /; NumberQ[a//N] && NumberQ[b//N] (* End Block for CycloidEqns *)

End[]

Protect[Race,CycloidEqns];

EndPackage[]