Simple physics in F# and WPF

0_README.md

Simple physics in F# and WPF

This is an example on how to use WPF and F# to do some simple physics using Verlet Integration.

The program creates a system of particles and constraints. Particles have inertia and is affected by gravity but their motion is also contrained by the constraints.

I tried to annotate the source code to help guide a developer familiar with languages like C#. If you have suggestions for how to improve it please leave a comment below.

How to run

1. WPF requires a Windows box
2. Install dotnet: https://dotnet.microsoft.com/en-us/download
3. Create a folder named for example: FsPhysics
4. Create file in the folder named: FsPhysics.fsproj and copy the content of 1_FsPhysics.fsproj below into that file
5. Create file in the folder named: Program.fs and copy the content of 2_Program.fs below into that file
6. Launch the application in Visual Studio or through the command line dotnet run from the folder FsPhysics
7. Should look like the image in this tweet.
8. Expand the particle system at lines 140+

1_FsPhysics.fsproj

<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>WinExe</OutputType>
    <TargetFramework>net6.0-windows</TargetFramework>
    <UseWPF>true</UseWPF>
  </PropertyGroup>

  <ItemGroup>
    <Compile Include="Program.fs" />
  </ItemGroup>

</Project>

2_Program.fs

// Hi!. The particle system is defined at row: 140+

// `open` are F# version of C# `using`
open System
open System.Diagnostics
open System.Globalization
open System.Numerics
open System.Windows
open System.Windows.Media
open System.Windows.Media.Animation

open FSharp.Core.Printf

type V1 = float32
type V2 = Vector2

// F# inline is sometimes used for performance but often
//  it's used to get access to more advanced generics than
//  supported by .NET CLR
//  x here can be any type that supports conversion to float32
let inline v1 x   = float32 x
let inline v2 x y = V2 (float32 x, float32 y)
let v2_0          = v2 0.F 0.F

// Define a particle record
//  Mass is stored in difference ways to avoid recomputing it
//  Current is the current position
//  Previous is the previous position
//  The speed then implicitly is Current-Previous
//  This representation is used in something called Verlet Integration
//  Verlet Integration avoids needing to update the speed vector
//  when computing the constraints
//  It doesn't produce accurate physics but it looks believable
//  which is good enough for this program
//  Verlet Integration is described with some detail here:
//    https://en.wikipedia.org/wiki/Verlet_integration
type Particle =
  {
    Mass              : V1
    SqrtMass          : V1
    InvertedMass      : V1
    mutable Current   : V2
    mutable Previous  : V2
  }

  // Verlet step moves the particle with inertia and gravity
  member x.Step (gravity : V1) =
    // InvertedMass of 0 means this is a fixed particle of infinite
    //  mass. These particles don't move
    if x.InvertedMass > 0.F then
      let c = x.Current
      let g = v2 0.0F gravity
      x.Current   <- g + c + (c - x.Previous)
      x.Previous  <- c

// Makes a particle given mass, position x,y and velocity vx,vy
let inline mkParticle mass x y vx vy : Particle =
  let m = v1 mass
  let c = v2 x y
  let v = v2 vx vy
  {
    Mass          = m
    InvertedMass  = 1.F/m
    SqrtMass      = sqrt m
    Current       = c
    Previous      = c - v
  }

// Makes a fix particle position x,y
//  a fix particle has infinite mass and doesn't move
//  used as an anchor point for other particles and constraints
let inline mkFixParticle x y = mkParticle infinityf x y 0.F 0.F

// Defines a constraint which is either a stick or a rope
//  a stick tries to make sure that the distance between two particles
//  are the Length value
//  a rope tries to makes sure that distance between two particles
//  are at most the Length value
type Constraint =
  {
    IsStick         : bool
    Length          : V1
    Left            : Particle
    Right           : Particle
  }

  // After the verlet step most constraints are "over stretched"
  //  Relax moves the two particles so that the constraint is "relaxed"
  //  again. This will in turn make other constraints "over stretched"
  //  but it turns out applying Relax over and over moves the system
  //  to a relaxed state
  member x.Relax () =
    // Bunch of math but the intent is this:
    //  compute the distance between the two particles in the constraint
    //  if the distance is not the right distance
    //  then move the two particles towards or away from eachother
    //  so that the distance is correct
    //  The comparitive mass of the particles is used to make sure
    //  that a small particle moves more than the bigger one it's 
    //  connected to
    let l  = x.Left
    let r  = x.Right
    let lc = l.Current
    let rc = r.Current

    let diff  = lc - rc
    let len   = diff.Length ()
    let ldiff = len - x.Length
    let test  = if x.IsStick then abs ldiff > 0.F else ldiff > 0.F
    if test then
      let imass = 0.5F/(l.InvertedMass + r.InvertedMass)
      let mdiff = (imass*ldiff/len)*diff
      let loff  = l.InvertedMass * mdiff
      let roff  = r.InvertedMass * mdiff

      l.Current <- lc - loff
      r.Current <- rc + roff

// Makes a stick constraint between two particles
let inline mkStick (l : Particle) (r : Particle) : Constraint =
  {
    IsStick   = true
    Length    = (l.Current - r.Current).Length ()
    Left      = l
    Right     = r
  }

// Makes a rope constraint between two particles
//  allows making the rope a bit longer than the initial distance
let inline mkRope extraLength (l : Particle) (r : Particle) : Constraint =
  {
    IsStick   = false
    Length    = (1.0F + abs (float32 extraLength))*(l.Current - r.Current).Length ()
    Left      = l
    Right     = r
  }

// Creates a small system of particles and constraints
//  That demonstrate connected swinging particles
let particles = 
  // 3 particles, one is fixed
  [|
    //            PosX  PosY
    mkFixParticle 0. -400.
    //            Mass  PosX  PosY  SpeedX  SpeedyY
    mkParticle    100.  400.  -400. 0.      0.
    mkParticle    200.  800.  -400. 0.      0.
  |]
let constraints =
  // Helper function to allow us creating stick constraints
  //  by using particle indices
  let mkStick l r = mkStick particles.[l] particles.[r]
  // Two constraints
  //  connected first to second then second to third
  [|
    //        Left  Right
    mkStick   0     1
    mkStick   1     2
  |]

// Creates a CanvasElement class that will act like a canvas for us
//  We override the OnRender method to draw graphics. In order to make the graphics
//  animate we have a time animation that invalidates the element which forces a redraw
type CanvasElement () =
  class
    // This is how in F# we inherit, this is typically not done as much
    //  as in C# but in order to be part of WPF Visual tree we need to 
    //  inherit UIElement
    inherit UIElement ()

    // Declaring a DependencyProperty member for Time
    //  This is WPF magic but it's created so that we can create
    //  an "animation" of the time value.
    //  This will help use do smooth updates.
    //  Nothing like web requestAnimationFrame in WPF AFAIK
    static let timeProperty =
      let pc = PropertyChangedCallback CanvasElement.TimePropertyChanged
      let md = PropertyMetadata (0., pc)
      DependencyProperty.Register ("Time", typeof<float>, typeof<CanvasElement>, md)

    // Freezing resources prevents updates of WPF Resources
    //  Can help WPF optimize rendering
    //  #Freezable is like C# constraint : where T : Freezable
    let freeze (f : #Freezable) = 
      f.Freeze ()
      f

    // Helper function to create pens
    let makePen thickness brush = 
      Pen (Thickness = thickness, Brush = brush) |> freeze

    let particlePen = makePen 2. Brushes.White
    let stickPen    = makePen 2. Brushes.Yellow
    let ropePen     = makePen 2. Brushes.GreenYellow

    // More WPF dependency property magic
    //  Not very interesting but this becomes member function in the class
    static member TimePropertyChanged (d : DependencyObject) (e : DependencyPropertyChangedEventArgs) =
      let g = d :?> CanvasElement
      // Whenever time change we invalidate the entire canvas element
      g.InvalidateVisual ()

    // Idiomatically WPF Dependency properties should be readonly 
    //  static fields. However, F# don't allow us to declare that
    //  Luckily it seems static readonly properties works fine
    static member TimeProperty = timeProperty

    // Gets the Time dependency property
    member x.Time = x.GetValue CanvasElement.TimeProperty :?> float

    // Create an animation that animates a floating point from 0 to 1E9
    //  over 1E9 seconds thus the time. This animation is then hooked onto the Time property
    //  Basically more WPF magic
    member x.Start () =
      // Initial time value
      let b   = 0.
      // End time, application animation stops after approx 30 years
      let e   = 1E9
      let dur = Duration (TimeSpan.FromSeconds (e - b))
      let ani = DoubleAnimation (b, e, dur) |> freeze
      // Animating Time property
      x.BeginAnimation (CanvasElement.TimeProperty, ani);

    // Finally we get to the good stuff!
    //  dc is a DeviceContext, basically a canvas we can draw on
    override x.OnRender dc =
      // Get the current time, will change over time (hohoh)
      let time  = x.Time
      // This is the size of the canvas in pixels
      let rs    = x.RenderSize

      let center= v2 (0.5*rs.Width) (0.5*rs.Height)

      // Apply the verlet step to all particles
      for p in particles do
        p.Step 0.1F

      // Relax all constraints 5 times
      //  If you relax less times the system becomes more "bouncy"
      //  More times makes it more "stiff"
      for _ = 1 to 5 do
        for c in constraints do
          c.Relax ()

      // inline here allows us to create helper function that
      //  uses a local variable without the overhead of creating
      //  a new function object
      //  Creating a bunch of objects during drawing can lead
      //  to GC which we like to avoid
      let inline toPoint (p : Particle) =
        let pos = p.Current + center
        Point (float pos.X, float pos.Y)

      // Draw all constraints
      for c in constraints do
        let pen = if c.IsStick then stickPen else ropePen
        dc.DrawLine (pen , toPoint c.Left, toPoint c.Right)

      // Draw all particles
      for p in particles do
        let r, b = 
          if p.InvertedMass = 0.F then 
            10., Brushes.White
          else
            let r = 3.0F + p.SqrtMass |> float
            r, Brushes.Black
        dc.DrawEllipse (b, particlePen, toPoint p, r, r)

  end

// Tells F# that this method is the main entry point
[<EntryPoint>]
// More 1990s magic! Basically in Windows there's a requirement that
//  UI controls runs in something called a Single Threaded Apartment.
//  So we tell .NET that the thread that calls main should be in a
//  Single Threaded Apartment.
//  Basically MS idea in 1990s on how to solve the problem of writing
//  multi threaded applications.
//  The .NET equivalent to apartments could be SynchronizationContext
[<STAThread>]
let main argv =
  // Sets up the main window
  let window  = Window (Title = "FsPhysics", Background = Brushes.Black)
  // Creates our canvas
  let element = CanvasElement ()
  // Makes our canvas the content of the Window
  window.Content <- element
  // Starts the time animation
  element.Start ()
  // Shows the Window
  window.ShowDialog () |> ignore
  0