14,608,247 members

Rate this:
20 Feb 2018CPOL
Calculate a wheel profile matching a road profile.

## 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:

1. A formula to draw the wheel.
2. A method to find dimensions of the wheel (ie the right H height of the wheel, see formula).
3. A mechanism to store the expressions of the road (=ground).
4. A control to evaluate a mathematic expression.
5. 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=H-y(x) \Rightarrow dr=-dy \end{aligned}

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}{H-y(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 (H-y) : 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}{H-y(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 :

1. 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}
2. Calculate points : see CalculatePoints() method.

CalculateHeight() method (details are in comments):

/// <summary>
/// <para>Calculates the height of the wheel hub in order to the wheel length matches n periods.</para>
/// <para>Dichotomy search.</para>
/// </summary>
/// <param name="ppl">Curve point list.</param>
/// <param name="nbPeriod">Number of period.</param>
/// <param name="maxGround">Maximum of the road function</param>
/// <returns></returns>
public static double CalculateHeight(
FunctionType typefunction,
List<DefinitionFunction> pfeGroudFunctions,
int nbPeriod,
double maxGround)
{

// Initial interval to begin search : maxGround<height<nL
// L : ground curve length on one period.
// interval min for the dichotomy
double heightmin = maxGround;
// interval max for the dichotomy = nbPeriod*[length of the curve for one period].
double heightmax = nbPeriod * CalculateCurveLength(typefunction, pfeGroudFunctions);

var h = ConstantKeys.CALCULATION_STEP;
double height = 0;          // current calculated height (=H)
double total_theta = 0;     // current angle
var isHauteurOK = false;    // flag indicating if height was found
double iteration = 0;       // number of iteration

// Expression compiler component
var myfonc = FonctionFactory.CreateFonction(typefunction, pfeGroudFunctions);

// Minimum of the definition interval of the ground function
var defmin = myfonc.Interval.Min;
// Maximum corresponding to nbPeriod in order to calculate all the points of the wheel.
var defmax = nbPeriod * (myfonc.Interval.Max - defmin) + defmin;

while (!isHauteurOK && iteration < ConstantKeys.DICHOTOMY_ITERATION_MAX)
{
iteration++;

// Height to test = average of the bounds
height = (heightmin + heightmax) / 2;

total_theta = 0;

var firstPoint = myfonc.Evaluate(defmin); // first point
var pointPrev = firstPoint; // previous point

// In the wheel point list, property X stores the angle, the property Y stores the radius.
// we go through the entire definition interval
for (var t = defmin+h; t <= defmax; t+=h)
{
var point = myfonc.Evaluate(t);     // current point

if (radius < 0) break;              // problem if length<0!
{
// formula dƟ = dx / (height-y(x)).
total_theta += (point.X - pointPrev.X) / radius;
}

pointPrev = point;
}

// Two double values can't be exactly equal, we will use h as precision.
if (total_theta >= (Math.PI * 2 - h) && total_theta <= (Math.PI * 2 + h))
{
isHauteurOK = true;
}
else
{
// we redefine the dichotomy interval
{
if (total_theta > Math.PI * 2)
{
// if angle > 2*pi then we have to increase the height because angle is inversely proportional to height
// according to formula dƟ=dx/height.
// => interval min is increased
heightmin = height;
}
else
{
// if angle < 2*pi then we have to decrease the height because angle is inversely proportional to height
// according to formula dƟ=dx/height.
// => interval max is decreased.
heightmax = height;
}
}
else
{
heightmin = height;
}
}

} // while end

if (isHauteurOK)
{
// just for information
TraceMessage.TraceInfoMessage(log,
String.Format(Resources.MsgNumberIteration_0_1, iteration, nbPeriod));
return height;
}

return -1;
}

/// <summary>
/// Calculate the curve length on one period.
/// </summary>
/// <param name="typefunction">Function type.</param>
/// <returns></returns>
public static double CalculateCurveLength(FunctionType typefunction, List<DefinitionFunction> dfs)
{
// Expression compiler component
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);
// pythagore formula
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;

// 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 X-Axis
: ConstantKeys.CALCULATION_STEP;

// First point
var point = myfonc.Evaluate(x);
var xprev = point.X;    // previous point

// Other points
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.
{
// increment of theta according to formula dƟ=dx/H
total_theta += (point.X - xprev) / radius;
}

// Radius for the nex point

// Backup the current point
xprev = point.X;

// Add the point to the list if the radius is positive.
{
}

// 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="x-1" />
<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="(t-sin(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 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
{
}
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;  // see <a href="https://www.codeproject.com/Articles/19768/Flee-Fast-Lightweight-Expression-Evaluator">"Walkthrough: Creating and Evaluating an Expression"</a>

private double intervalMin = double.MaxValue;
private double intervalMax = double.MinValue;

public AbstractFonction(FunctionType typeFunction, IEnumerable<DefinitionFunction> expressions)
{
TypeFunction = typeFunction;
eo = new ExpressionOwner();

// Define the context of our expression
ec = new ExpressionContext(eo);

eDynamics = new Dictionary<Interval, IGenericExpression<double>>();

foreach (var expr in expressions)
{
}

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 no function defined then return NaN
if (eDynamics.Count == 0) return new Point(x,double.NaN);

// intermediate variable to calculate the function outside of its period.
var t = x;
// Ground functions are periodic.
t = SetVariableInInterval(x);

eo.X = t;

// If just one function is defined
if (eDynamics.Count == 1)
{
return new Point(x,eDynamics.First().Value.Evaluate());
}

// For piecewise continuous function
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 user-modifiable as well). The mouse-wheel can also be used for zooming. There is also a context menu that allows you to un-zoom and un-pan 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 X-axis and Y-axis when zooming or panning in GraphComponent class:

private void Zgc_ZoomEvent(ZedGraphControl sender, ZoomState oldState, ZoomState newState)
{
// Zoom event raised by zedgraph.
// zoom event is raised when user pan or zoom the first pane in the MasterPane.
// When zooming, we restrict the minimum of the axis to begin at 0 or lower.

Pane.XAxis.Scale.Min = Math.Min(0, Pane.XAxis.Scale.Min);
Pane.YAxis.Scale.Min = Math.Min(0, Pane.YAxis.Scale.Min);

// Redraw curves
}

In the MyZedGraph control, I added method to make coordinate system orthonormal and this method is called in OnResize event.

public void InitScaleOrthonormal()
{
// First pane in the Master pane
var pane = this.GraphPane;
// The rectangle in which the curves are displayed
var rect = pane.Chart.Rect;

// if window is minimize
if (rect.Width < 0 || rect.Height < 0) { this.Refresh(); return; }

// ratio
var ratio = rect.Height / rect.Width;

// We calculate the maximum on the Y Axis to have an orthonormal coordinate system relative to the X axis dimensions.
// Ymin is not recalculate because it can come from a result of a zoom event.
pane.YAxis.Scale.Max = pane.YAxis.Scale.Min + ratio * (pane.XAxis.Scale.Max - pane.XAxis.Scale.Min);

// Indicates that axis changed
pane.AxisChange();

this.Refresh();
}

## 6. Some interesting solutions

Description Screen capture
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).

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.

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 grid.
Choice of interface languages (French or English).
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.
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.
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/Flee-Fast-Lightweight-Expression-Evaluator

https://flee.codeplex.com/

https://archive.codeplex.com/?p=flee

https://sourceforge.net/projects/zedgraph/

http://zedgraph.sourceforge.net/linesamples.html

A flexible charting library for .NET : https://www.codeproject.com/Articles/5431/A-flexible-charting-library-for-NET

https://msdn.microsoft.com/en-us/library/system.configuration.configurationelementcollection(v=vs.110).aspx : ConfigurationElementCollection Class

## History

Version 1.0 of the program.

## Share

 Software Developer (Senior) France
I have a Master's degree in Electronics Engineering. My job is software developer since 1996. I like programming and I use the c#, javascript, vba, php languages. I made websites in php and ASP.NET, programs in vb.net and c#, plugins to Office, mobile project (for bar code reader) and webparts for SharePoint.
I like painting and drawing, and sculpting clay.

 First Prev Next
 Interesting Chrris Dale18-Mar-18 11:48 Chrris Dale 18-Mar-18 11:48
 Great! Simply great! Friedel H. Hartmann 22-Feb-18 7:21 Friedel H. Hartmann 22-Feb-18 7:21
 Re: Great! Simply great! Christian Dalle23-Feb-18 4:05 Christian Dalle 23-Feb-18 4:05
 I like it phil.o22-Feb-18 6:55 phil.o 22-Feb-18 6:55
 Re: I like it Christian Dalle23-Feb-18 4:11 Christian Dalle 23-Feb-18 4:11
 Missing Images in Article Saineshwar Bageri20-Feb-18 5:58 Saineshwar Bageri 20-Feb-18 5:58
 Re: Missing Images in Article Christian Dalle20-Feb-18 9:37 Christian Dalle 20-Feb-18 9:37
 Last Visit: 15-Aug-20 4:36     Last Update: 15-Aug-20 4:36 Refresh 1