Wheel-Road couple

• Download source - 1.5 MB
• Download Wheelprofile_bin.zip - 862.5 KB

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:

And we have (|\textbf{e}_\textbf{r}|=|\textbf{e}_\textbf{θ}|=1 \mbox{ and } \textbf{e}_\textbf{r}\textbf{e}_\textbf{θ}=0)

Road:

Rolling without sliding condition:

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
    double radius = 0;          // current radius
    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

            radius = height - pointPrev.Y;      // previous radius
            if (radius < 0) break;              // problem if length<0!
            if (radius > 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(radius > 0)
            {
                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
            {
                 // we have to increase the radius because radius is negative.
                 heightmin = height;
            }
        }

    } // while end

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

    // height was not found.
    return -1;
}

/// <summary>
/// Calculate the curve length on one period.
/// </summary>
/// <param name="typefunction">Function type.</param>
/// <param name="dfs">Road functions.</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;
}

/// <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 X-Axis
              : 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="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 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;  // 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);
        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).

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
   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()
{
   // 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
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/Flee-Fast-Lightweight-Expression-Evaluator

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/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.