Introduction
When I was younger, I saw a program on a Casio handheld computer (year 80) that drew a wheel profile corresponding to a road profile. The principle is to roll the wheel on the road so that the hub of the wheel moves horizontally. In the screen capture, yellow line is the path of the wheel, red line displays the position of the hub. The example taken was a crenel road. The program drew three wheels, one at the beginning of the road, one in the middle and one at the end. The road function was entered in assembler directly into the program memory area.
So, the goal of this project is to do at minima the same with additional features (see below).
To realize this project, we must go through these steps:
 A formula to draw the wheel.
 A method to find dimensions of the wheel (ie the right H height of the wheel, see formula).
 A mechanism to store the expressions of the road (=ground).
 A control to evaluate a mathematic expression.
 A control to draw the curves.
Point of interest
 Configuration classes to store function expressions.
 Ciloci.Flee.dll, an expression compiler C#: use for writing function expressions.
 Animation of the wheel.
 Dichotomy method to find dimensions of the wheel to match a number of road period.
1. A little mathematics to find formula
To start, we need a formula to draw the wheel. Below, we explain how we find the formula.
Wheel:
(Bold characters are vectors).
Coordinate system ((\textbf{e}_\textbf{r},\textbf{e}_\textbf{θ})) is linked to the wheel.
Polar coordinates for infinitesimal shifting: (d\textbf{M}=dr\textbf{e}_\textbf{r}+rdθ\textbf{e}_\textbf{θ})
Infinitesimal length on the wheel:
$\begin{aligned}d\textbf{M}.d\textbf{M}=ds_2^2=(dr\textbf{e}_\textbf{r}+rdθ\textbf{e}_\textbf{θ})^2\end{aligned}$
And we have (\textbf{e}_\textbf{r}=\textbf{e}_\textbf{θ}=1 \mbox{ and } \textbf{e}_\textbf{r}\textbf{e}_\textbf{θ}=0)
Therefore $\begin{aligned} ds_2^2=(dr\textbf{e}_\textbf{r}+rdθ\textbf{e}_\textbf{θ})^2=dr^2+r^2dθ^2 \end{aligned}$
With $\begin{aligned} r=Hy(x) \Rightarrow dr=dy \end{aligned}$
Road:
An infinitesimal length on the road: $\begin{aligned} ds_1^2=dx^2+dy^2=dx^2(1+\frac{dy^2}{dx^2})=dx^2(1+(\frac{dy}{dx})^2) \end{aligned}$
Rolling without sliding condition:
$\begin{aligned} ds_1=ds_2 \Rightarrow dx^2(1+(\frac{dy}{dx})^2) &= dr^2+r^2d\theta^2\\ 1+(\frac{dy}{dx})^2 &= \frac{dr^2}{dx^2}+r^2\frac{d\theta^2}{dx^2}\\ 1+(\frac{dy}{dx})^2 &= (\frac{dr}{dx})^2+r^2(\frac{d\theta}{dx})^2\\ 1+(\frac{dy}{dx})^2 &= (\frac{dy}{dx})^2+r^2(\frac{d\theta}{dx})^2\\ 1 &= r^2(\frac{d\theta}{dx})^2\\ \frac{d\theta}{dx} &= \frac{1}{r} \end{aligned}$
Finally, we have (\boxed{d\theta=\frac{dx}{Hy(x)}})
With this formula, we can construct step by step the wheel curve.
Notes:
 The road can be drawn above or below the X axis: it means that it can be take negative values.
 The road function is a piecewise continuous function.
 The road is supposed to be periodic. Thus, the road functions definition set is defined on an interval of length T (=period)
 Obviously, radius of the wheel must be positive (Hy) : it means that H is greater than the sup(y(x)) (maximum of y(x)).
 Wheel curve is entirely drawn when θ goes through the interval [0,2π]
2. How we calculate points of the wheel?
The formula (d\theta=\frac{dx}{Hy(x)}) assumes that the quantities H and y(x) are known. y(x) is the function of the road but H can take all values between sup(y(x)) and infinity. If we take any random value, the wheel curve may not be closed. Therefore, we are going to calculate H to match a certain number of period of the road: in the configuration file, the number of period corresponds to nbPeriodWheel property (see Configuration) and in the screen capture it corresponds to Number of period (in the Wheel control group).
The WheelComponent
class contains methods to calculate points of the wheel curve. The steps are :
 Calculate H with the dichotomy method : see
CalculateHeight()
method. We need interval for dichotomy, so we considered that H is above the road curve and less than the length of the curve on nbperiod: $\begin{aligned} sup(f(t))\leq H\leq \int_{t_{min}}^{t_{nbperiod}}f'(t)dt\end{aligned}$  Calculate points : see
CalculatePoints()
method.
CalculateHeight()
method (details are in comments):
public static double CalculateHeight(
FunctionType typefunction,
List<DefinitionFunction> pfeGroudFunctions,
int nbPeriod,
double maxGround)
{
double heightmin = maxGround;
double heightmax = nbPeriod * CalculateCurveLength(typefunction, pfeGroudFunctions);
var h = ConstantKeys.CALCULATION_STEP;
double height = 0;
double total_theta = 0;
double radius = 0;
var isHauteurOK = false;
double iteration = 0;
var myfonc = FonctionFactory.CreateFonction(typefunction, pfeGroudFunctions);
var defmin = myfonc.Interval.Min;
var defmax = nbPeriod * (myfonc.Interval.Max  defmin) + defmin;
while (!isHauteurOK && iteration < ConstantKeys.DICHOTOMY_ITERATION_MAX)
{
iteration++;
height = (heightmin + heightmax) / 2;
total_theta = 0;
var firstPoint = myfonc.Evaluate(defmin);
var pointPrev = firstPoint;
for (var t = defmin+h; t <= defmax; t+=h)
{
var point = myfonc.Evaluate(t);
radius = height  pointPrev.Y;
if (radius < 0) break;
if (radius > 0)
{
total_theta += (point.X  pointPrev.X) / radius;
}
pointPrev = point;
}
if (total_theta >= (Math.PI * 2  h) && total_theta <= (Math.PI * 2 + h))
{
isHauteurOK = true;
}
else
{
if(radius > 0)
{
if (total_theta > Math.PI * 2)
{
heightmin = height;
}
else
{
heightmax = height;
}
}
else
{
heightmin = height;
}
}
}
if (isHauteurOK)
{
TraceMessage.TraceInfoMessage(log,
String.Format(Resources.MsgNumberIteration_0_1, iteration, nbPeriod));
return height;
}
return 1;
}
public static double CalculateCurveLength(FunctionType typefunction, List<DefinitionFunction> dfs)
{
var myfonc = FonctionFactory.CreateFonction(typefunction, dfs);
var defmin = myfonc.Interval.Min;
var defmax = myfonc.Interval.Max;
double h = ConstantKeys.CALCULATION_STEP;
double L = 0;
var pointPrev = myfonc.Evaluate(defmin);
for (var t = defmin+h; t <= defmax+h; t += h)
{
var point = myfonc.Evaluate(t);
L += Math.Sqrt((point.X  pointPrev.X)*(point.X  pointPrev.X) + (point.Y  pointPrev.Y)*(point.Y  pointPrev.Y));
pointPrev = point;
}
return L;
}
CalculatePoints()
method (details are in comments):
/// <summary>
/// <para>Calculate the wheel points. The curve is closed and the period is 2*PI.</para>
/// <para>The first point of the wheel is theta = 0 corresponding to x=xmin.</para>
/// </summary>
/// <param name="ProfileGround">Ground profile.</param>
/// <param name="hauteur">Height of the wheel hub.</param>
/// <returns>ProfileFunctionWheel object</returns>
public static ProfileFunctionWheel CalculatePoints(ProfileFunctionExt ProfileGround, double height)
{
var ret = new ProfileFunctionWheel()
{
TypeFunction = Enums.FunctionType.Polar,
Name = ConstantKeys.CURVE_WHEEL_NAME,
TypeCurve = Enums.CurveType.Wheel,
Period = 2*Math.PI
};
// init variables
double total_theta = 0;
double radius = 0;
// Function of the ground
var myfonc = FonctionFactory.CreateFonction(ProfileGround.TypeFunction, ProfileGround.DefinitionFunctions);
// Lower bound of the function interval
var defxmin = myfonc.Interval.Min;
// distance from xmin
double x = defxmin;
// step on X axis.
// 
// defxmin defxmin+h defxmin+2h
var h = ProfileGround.TypeFunction == FunctionType.Cartesian
? ProfileGround.Points[1].X  ProfileGround.Points[0].X // distance between two road curve points on XAxis
: ConstantKeys.CALCULATION_STEP;
// First point
var point = myfonc.Evaluate(x);
radius = height  point.Y;
var xprev = point.X; // previous point
ret.Points.Add(new PointPair(total_theta, radius));
// Other points
// theta = h/(previous radius)
while (total_theta <= (2 * Math.PI))
{
// We move along the X axis by increment of h
x += h;
// Function evaluation
point = myfonc.Evaluate(x);
// radius must be positive : if it's negative, we leave the point and go ahead.
// here is the previous radius.
if (radius > 0)
{
// increment of theta according to formula dƟ=dx/H
total_theta += (point.X  xprev) / radius;
}
// Radius for the nex point
radius = height  point.Y;
// Backup the current point
xprev = point.X;
// Add the point to the list if the radius is positive.
if (radius > 0)
{
// Add current radius with previous total_theta
ret.Points.Add(total_theta, radius);
}
// Prevent from infinite loop or to much points.
if (ret.Points.Count > ConstantKeys.MAX_NB_POINTS_WHEEL) break;
}
ret.Centroid = CalculateCentroid(ret.Points);
return ret;
}
Here is a schema which shows the algorithm above:
3. Configuration
Configuration classes stores the expressions of the road. The classes reflect the structure of the application configuration file.
Extract of the App.config file (custom tags are defined in the appconfig.xsd file):
<configSections>
<section name="DefinitionFunctionsSection" type="WheelProfile.Configuration.DefinitionFunctionsConfigurationSection,WheelProfile" />
</configSections>
<DefinitionFunctionsSection xmlns="http://tempuri.org/appconfig.xsd">
<DefinitionFunctions>
<DefinitionFunction name="Triangle" period="2" typefunction="Cartesian" nbPeriodWheel="7">
<items>
<item min="0" max="T/2" expression="x1" />
<item min="T/2" max="T" expression="x+1" />
</items>
</DefinitionFunction>
<DefinitionFunction name="Cycloid" period="2*pi" typefunction="Parametric" nbPeriodWheel="3">
<items>
<item min="0" max="T" expression="(tsin(t))" />
<item min="0" max="T" expression="2*sin(t/2)*sin(t/2)" />
</items>
</DefinitionFunction>
[...]
</DefinitionFunctions>
</DefinitionFunctionsSection>
Note :
 The appconfig.xsd is used to have completion in Visual Studio for custom tags when modifying app.config file.
 Parametric function use variable
t
in the formula instead of x
. The first expression is for x(t) and the second for y(t). The intervals have the same bounds. The program does not deal with the case of piecewise parametric functions.
We store the properties:
Name
: the name of the functions that represents the road. Period
: the period of the functions. You can use any of the constants defined in Math
class or a formula (like 2*pi). It is important not to be wrong, otherwise we may not find all the wheel profiles.
Example: sin(2x) has a period of pi and a wheel corresponding to one period is an ellipse. If you set the period to 2pi the road is well drawn but the wheel profile corresponding to one period is the peanut!
T=pi, number period of wheel=1  T=2pi, number period of wheel=1 
 
Typefunction
: it can be Cartesian or Parametric. (Polar is reserved for the wheel). If functions are parametric, the min and max for the two functions must be the same (see Cycloid in the app.config file). nbPeriodWheel
: default nbPeriod value for drawing the first wheel. Min
: Minimum of the interval for the function. You can use a formula like T/2, 2*pi. T represents the property Period
. Max
: Maximum of the interval for the function. You can use a formula like T/2, 2*pi. Expression
: Mathematic function of the road.
To store the properties, we define some classes:
ConfigurationClass.cs
: Abstract class that defines some useful methods. MyConfiguration.cs
: Singleton configuration class which maps the .config file. MyConfigurationSection.cs
: Base class for configuration section classes. It defines namespace "xmlns" for using my own schema (appconfig.xsd). DefinitionFunctionsConfigurationSection.cs
: Class to define section "DefinitionFunctions". Inherites from MyConfigurationSection. DefinitionFunctionsCollection.cs
: Contains multiples classes: ItemDefinitionFunctionsElement
, DefinitionFunctionsElement
, DefinitionFunctionsCollection
and ItemDefinitionFunctionsCollection
Extract of MyConfiguration.cs
:
public class MyConfiguration : ConfigurationClass
{
private static readonly log4net.ILog log = log4net.LogManager.GetLogger(MethodBase.GetCurrentMethod().DeclaringType);
private const string DefinitionFunctionsSectionAttribute = "DefinitionFunctionsSection";
private const string DefaultDelayThreadAttribute = "DefaultDelayThread";
private const string QuickDelayThreadAttribute = "QuickDelayThread";
private static MyConfiguration _current;
private static readonly object SyncRoot = new Object();
public static MyConfiguration Current
{
get
{
if (_current == null)
{
lock (SyncRoot)
{
if (_current == null)
_current = new MyConfiguration();
}
}
return _current;
}
}
public DefinitionFunctionsConfigurationSection DefinitionFunctionsSection { get; private set; }
public int DefaultDelayThread { get; private set; }
public int QuickDelayThread { get; private set; }
private MyConfiguration()
{
DefinitionFunctionsSection = (DefinitionFunctionsConfigurationSection)ConfigurationManager.GetSection(DefinitionFunctionsSectionAttribute);
try
{
DefaultDelayThread = Convert.ToInt32(ConfigurationManager.AppSettings[DefaultDelayThreadAttribute]);
QuickDelayThread = Convert.ToInt32(ConfigurationManager.AppSettings[QuickDelayThreadAttribute]);
}
catch (Exception ex)
{
TraceMessage.TraceErrorMessage(log, ex.Message);
}
}
4. Expression compiler
"Flee is a .NET library that allows you to parse and evaluate arbitrary expressions. The main feature that distinguishes it from other expression evaluators is that it uses a custom compiler to convert expressions into IL and then, using Lightweight CodeGen, emits the IL to a dynamic method at runtime." [code project article]
Indeed, this evaluator is very fast. I create an abstract class AbstractFonction
which implements IFonction
interface. This abstract class contains all the mechanism to use the Ciloci.Flee component.
protected ExpressionContext ec = null;
protected Dictionary<Interval, IGenericExpression<double>> eDynamics = null;
protected ExpressionOwner eo = null;
private double intervalMin = double.MaxValue;
private double intervalMax = double.MinValue;
public AbstractFonction(FunctionType typeFunction, IEnumerable<DefinitionFunction> expressions)
{
TypeFunction = typeFunction;
eo = new ExpressionOwner();
ec = new ExpressionContext(eo);
ec.Imports.AddType(typeof(Math));
eDynamics = new Dictionary<Interval, IGenericExpression<double>>();
foreach (var expr in expressions)
{
eDynamics.Add((Interval)expr, ec.CompileGeneric<double>(expr.Expression));
}
intervalMin = Math.Min(intervalMin, expressions.First().Min);
intervalMax = Math.Max(intervalMax, expressions.Last().Max);
Period = intervalMax  intervalMin;
Interval = new Interval { Min = intervalMin, Max = intervalMax };
}
public FunctionType TypeFunction { get; private set; }
public double Period { get; private set; }
public Interval Interval { get; private set; }
public virtual Point Evaluate(double t) {}
Two other classes inherit from AsbtractFonction
:
CartesianFonctionComponent
: used for cartesian function y(x). ParametricFonctionComponent
: used for parametric function x(t), y(t).
A static class FonctionFactory
instantiates one or the other class depending on the function type. It returns IFonction
interface.
To evaluate a function, the Evaluate()
method is used (extract from CartesianFonctionComponent.Evaluate()
):
public override Point Evaluate(double x)
{
if (eDynamics.Count == 0) return new Point(x,double.NaN);
var t = x;
t = SetVariableInInterval(x);
eo.X = t;
if (eDynamics.Count == 1)
{
return new Point(x,eDynamics.First().Value.Evaluate());
}
foreach(var edyn in eDynamics)
{
var key = edyn.Key;
if(t >= key.Min && t < key.Max)
{
return new Point(x,edyn.Value.Evaluate());
}
}
return new Point(x, double.NaN);
}
5. ZedGraph
I compile ZedGraph with .NET 4.5. It's not the most appropriate control for the program but it's the better I found. In further version, I will use Research.DataDynamic
in a WPF program (but I need to learn WPF!)
I extend ZedGraph to add some methods and initialization of the principal pane.
What I found interesting is the feature of zooming and panning in the zoom event (ZoomEvent
). "To pan, either click with the middle mouse button or hold down the shift key, and left click to drag the graph around (the zoom/pan key combinations are usermodifiable as well). The mousewheel can also be used for zooming. There is also a context menu that allows you to unzoom and unpan to prior states, to restore the scale to full auto mode, to show point value tool tips, and to copy the graph (bitmap form) to the clipboard. Event options are provided to allow for customization of pan, zoom, point values, etc."[see article].
I forced the minimum on the Xaxis and Yaxis when zooming or panning in GraphComponent
class:
private void Zgc_ZoomEvent(ZedGraphControl sender, ZoomState oldState, ZoomState newState)
{
Pane.XAxis.Scale.Min = Math.Min(0, Pane.XAxis.Scale.Min);
Pane.YAxis.Scale.Min = Math.Min(0, Pane.YAxis.Scale.Min);
UIThreadInvoke(RefreshCurves);
}
In the MyZedGraph
control, I added method to make coordinate system orthonormal and this method is called in OnResize
event.
public void InitScaleOrthonormal()
{
var pane = this.GraphPane;
var rect = pane.Chart.Rect;
if (rect.Width < 0  rect.Height < 0) { this.Refresh(); return; }
var ratio = rect.Height / rect.Width;
pane.YAxis.Scale.Max = pane.YAxis.Scale.Min + ratio * (pane.XAxis.Scale.Max  pane.XAxis.Scale.Min);
pane.AxisChange();
this.Refresh();
}
6. Some interesting solutions
Description  Screen capture 
Road is a "reverse chain".
nbPeriod=4.
Wheel is a square.
 
Road is a periodic "reverse semi circles".
nbPeriod=1.
Wheel is a circle with hub on the circumference.
The curve is not closed because of the calculus precision of the wheel diameter (= hub height).
 
Road is a sinewave.
nbPeriod=1.
Wheel is an ellipse and the hub is at one of the focus.
If we increase the number of period (=2), the wheel becomes a "peanut" and the hub is at the center.
 
Road is a crenel.
 
Road is a "reverse cycloid".
nbPeriod=1.
Wheel is a cardioid.
 
User Interface
When user select a road name in the list, the program displays the road curve, the wheel curve and the expressions of the road functions.
Below, the list of commands:
Commands  Description 
Toolbar:  
 Displays/Hides the road. 
 Displays/Hides the grid. 
 Choice of interface languages (French or English). 
Road Equation area:  
 Save the current functions of the road (DefinitionFunction item), the number of period (see Wheel area) to the configuration file. 
 Delete the current road functions (DefinitionFunction item) from the configuration file. 
ComboBox  List all the road functions (DefinitionFunction items). 
GridView  Displays the functions. You can add a row, delete a row. An empty row is discarded. 
Period T=  Period of the road. The road is supposed to be periodic. 
Create Road  Button to generate the road and wheel curves. 
Context menu  Datagrid context menu to delete a row.

Wheel area:  
Number of period:  The number of period (period=T) corresponding to n tour of the wheel. 
Wheel height:  The height of the wheel hub. User can click the upbutton or downbutton to modify the wheel and see that if the height is not equal to a certain value, the curve is not closed. 
Calculate  Button to generate the wheel curve corresponding to the number of period defined. 
Animation:  Controls to animate the wheel. 
 Play the animation. (Default speed) 
 Pause the animation. 
 Stop the animation. The curve return to the start position. 
 Increase/decrease speed and play forward (when the animation is running). 
 Increase/decrease speed and play backward (when the animation is running). 
Graphic area:  
pan  To pan, click and drag with the middle mouse button 
zoom  To zoom, use the mouse wheel. if the wheel curve is bigger than the graphic area, you can zoom out. 
Bibliography
https://www.mathcurve.com/
Ciloci.Flee : expression parser:
https://www.codeproject.com/Articles/19768/FleeFastLightweightExpressionEvaluator
https://flee.codeplex.com/
https://archive.codeplex.com/?p=flee
A simple C# library for graph plotting (for precisionTimer.cs): https://www.codeproject.com/script/Articles/ViewDownloads.aspx?aid=32836
https://sourceforge.net/projects/zedgraph/
http://zedgraph.sourceforge.net/linesamples.html
A flexible charting library for .NET : https://www.codeproject.com/Articles/5431/AflexiblechartinglibraryforNET
https://msdn.microsoft.com/enus/library/system.configuration.configurationelementcollection(v=vs.110).aspx : ConfigurationElementCollection Class
History
Version 1.0 of the program.