Vectors, Projectiles, and More: Doing Physics in .net Part III

Having just started getting into C#, I can already say that it’s a pleasure to work with – especially when used with its tightly coupled IDE, Visual Studio. Of course, being a physics/astronomy/computer science geek, I got straight to work doing some projectile physics in C# to see how it turned out. To get myself started on this blog, I’m going to write a couple of posts about doing physics in .net:

• Units of Measure: Doing Physics in .net Part I
• Physics Library: Doing Physics in .net Part II
• Vectors, Projectiles, and More: Doing Physics in .net Part III

In part 2 I started out an F# physics library containing the units we want to use, and some unit-safe math functions called Units.fs (note the edit added to part 2 regarding the sqrt function – you don't need to write one, F# has its own built-in, unit-safe sqrt). This will be used by the rest of the library and should remain the first source file in your project (in Visual Studio, the order the files appear in is the order they are compiled in. Now we want to start fleshing out our library and making it into something useful. I'll start by moving the whole projectile class into the F# library.

But wait. I'm pretty much just writing this whole thing over again in F#, right? The point of this is to interoperate between the two languages, after all. Really, this is the proper thing to do – we're going to write the projectile physics as a class in F#. The class can be extended by C# or other .net languages if you want to add qualitative properties to the projectile (such as, making a red ball that you want to send flying through the air – you just attach the projectile physics to your otherwise normal red ball object).

So here's our C# class, translated into F#:

//Projectiles.fs
#light

namespace FsPhysicsDemo

module Projectiles =
    [<Literal>]
    let g = -9.808<m/s^2>

open Projectiles
open UnitSafeFunctions

type Projectile =
    val mutable xi : float<m>
    val mutable yi : float<m>
    val mutable vxi : float<m/s>
    val mutable vyi : float<m/s>

    new() = {
        xi = 0.0<m>; yi = 0.0<m>;
        vxi = 0.0<m/s>; vyi = 0.0<m/s>;
    }
    
    member this.xf (t : float<s>) =
        this.xi + this.vxi * t
    member this.yf (t : float<s>) =
        this.yi + this.vyi * t + g * sqr t

Here we've moved our g constant out of the class, and into a module. As for the Projectile class, you see that we opened the Projectiles and UnitSafeFunctions modules just above them. This is because we use the g constant, and our unit-safe square function in the last member function in the class. We don't have to open anything in order to use the units however, since those weren't defined in a module, but rather just in the namespace (our Projectile class is also defined at namespace level – types can be in a namespace).

That's not bad, now we have our C# class, but a little more condensed and unit-aware. We're not really harnessing much power from functional programming right now though—but what more can we do? Well, those x and y values all belong together. In C# you'd probably use the Point class in System.Drawing, or make your own class. In a functional programming language (or one that supports them nicely, such as Python) you could use a Tuple:

type Projectile =
    val mutable xi : float<m> * float<m>
    val mutable vi : float<m/s> * float<m/s>

    new() = { xi = (0.0<m>, 0.0<m>); vi = (0.0<m/s>, 0.0<m/s>); }

    member this.xf (t : float<s>) =
        let x = (fst this.xi) + (fst this.vi) * t
        let y = (snd this.xi) + (snd this.vi) * t + g * sqr t
        (x, y)

This works well and is more succinct and appropriate: x and y coordinates belong together. But I still thought we could do better—I see two issues with this approach. The first is that, while C# supports Tuples, it doesn't do so as nicely as F# does:

using System;
using FsPhysicsDemo;

namespace CsPhysicsDemo
{
    class Program
    {
        static void Main(string[] args)
        {
            Projectile ball = new Projectile()
            {
                xi = new Tuple<double, double>(42.0, 42.0),
                vi = new Tuple<double, double>(12.0, 24.0),
            };

            Tuple<double, double> xf = ball.xf(6.0);
            Console.WriteLine("Final X: " + xf.Item1);
            Console.WriteLine("Final Y: " + xf.Item2);
        }
    }
}

The second issue is that this function only takes Cartesian coordinates. This makes sense for setting the initial position of an object, but the initial velocity is frequently expressed using polar coordinates: that is, a magnitude and direction. We should accommodate this by creating a Vector type. This type will look vastly different from the Projectile type, though: we're going to use two very important concepts in functional programming:

• A Discriminated Union: Like an abstract parent type with a known set of subclasses – a DU has no C# equivalent and is very useful
• Pattern Matching: Like using a switch block in C#, except it makes sure all possible values of the test expression are being accounted for.

First lets look at the basic DU structure:

//Vectors.fs
#light

namespace FsPhysicsDemo

type Vector<[<Measure>] 'u> =
    | Cartesian of float<'u> * float<'u>
    | Polar of float<'u> * float

module useAVector =
    let a = Vector<m>.Cartesian(2.0<m>, 2.0<m>)
    let b = Vector<m>.Polar(2.83<m>, 0.393)

Now we can access these values by pattern matching – and when using pattern matching, we can't address a Cartesian vector without also addressing a Polar vector:

This is a good safety, but we won't be using it in this manner for Vectors. They would be unwieldly when dealing with the two vectors contained in our Projectile class:

    member this.xf (t : float<s>) =
        match xi with
        | Cartesian(x, y) -> 
            match vi with
            | Cartesian(x, y) -> //...
            | Polar(m, d) -> //...
        | Polar(m, d) -> 
            match vi with
            | Cartesian(x, y) -> //...
            | Polar(m, d) -> //...

So we'll have to re-think our Vector type. For every Cartesian vector there is a corresponding Polar vector, and vice-versa. So we can convert between the two—but that would mean we still have to match for each type of vector, so we know whether or not we need to convert it in a function. Also, that would mean each function written that manipulates vectors would have to choose which type to manipulate, and convert any other types to that type. This is a minor nitpick, but still seems wrong.

So how should we make the Vector type work? If this were a C# class, we might think of each subtype in the DU as different constructors to the same class. This works because Cartesian and Polar vectors are just two different ways of expressing the same data—so each constructor can sufficiently set up the object. From X and Y we can get magnitude and direction, and vice-versa—so we can just calculate all these values when the object is created, and then we don't have to match different types of vectors everywhere we want to use them:

open UnitSafeFunctions

type Vector<[<Measure>] 'u> =
    | Cartesian of float<'u> * float<'u>
    | Polar of float<'u> * float
    member this.x =
        match this with
        | Cartesian(x, _) -> x
        | Polar(m, d) -> m * cos d
    member this.y =
        match this with
        | Cartesian(_, y) -> y
        | Polar(m, d) -> m * sin d
    member this.magnitude =
        match this with
        | Cartesian(x, y) -> sqrt (sqr x + sqr y)
        | Polar(m, _) -> m
    member this.direction =
        match this with
        | Cartesian(x, y) -> tan (y / x)
        | Polar(_, d) -> d

Note that we could easily have made member this.cartesian = … and only had two members, but in the interest of saving lines wherever vectors are to be used, I opted to use a few extra lines here. Now our projectile class looks like so:

type Projectile =
    val mutable xi : Vector<m>
    val mutable vi : Vector<m/s>

    new() = { xi = Vector<m>.Polar(0.0<m>, 0.0); vi = Vector<m/s>.Polar(0.0<m/s>, 0.0); }

    member this.xf (t : float<s>) =
        Vector<m>.Cartesian(
            this.xi.x + this.vi.x * t,
            this.xi.y + this.vi.y * t + g * sqr t)

    member this.vf (t : float<s>) =
        Vector<m/s>.Cartesian(this.vi.x, this.vi.y + g * t)

And from C#:

class Program
{
    static void Main(string[] args)
    {
        Projectile ball = new Projectile()
        {
            xi = Vector.NewCartesian(42, 42),
            vi = Vector.NewPolar(5, Math.PI / 8),
        };

        Vector xf = ball.xf(6);
        Console.WriteLine("Final X: " + xf.x);
        Console.WriteLine("Final Y: " + xf.y);
        Vector vf = ball.vf(6);
        Console.WriteLine("Final Velocity: " + vf.magnitude);
        Console.WriteLine("Final Direction: " + vf.direction);
    }
}

Gorgeous. Now we can extend the Projectile class from C# and do whatever we want with it, and still retain our easy-to-verify math. We can also tell the difference between a Cartesian and polar vector in a strongly-typed way, and – although the units are lacking on the C# side, any further F# code that wants to use a Projectile (even a C# subclass of Projectile!) will have compile-time units of measure to ensure its correctness:

All-in-all, F# is a great fit for physics and physical equations. It's powerfully expressive nature makes it easy to write and bind together physical concepts and equations with little overhead, and its running on the .NET CLR lets you couple these abilities with rich, mature C#.