use Matrix math

Author: monkstone

processing_app/library/vecmath/vec2d/circumcircle.md

@@ -33,3 +33,26 @@ c = Vec2D.new(200, 200)
 (a - b).cross(b - c) == 0 # the area of the triangle is zero, so a, b, c are collinear
 
 ```
+
+Also because we were able to separate the logic, we were able to confidently re-factor the Barbara Almeida sketch to use Matrix math to determine the circumcenter
+
+### Matrix Math ###
+
+For detailed workings see [Circumcircle at Mathworld Wolfram.com][circumcircle]
+
+
+a = {{x<sub>1</sub> y<sub>1</sub> 1}, {x<sub>2</sub> y<sub>2</sub> 1}, {x<sub>3</sub> y<sub>3</sub> 1}}
+
+bx = -{{x<sub>1</sub><sup>2</sup> + y<sub>1</sub><sup>2</sup> y<sub>1</sub> 1}, {x<sub>2</sub><sup>2</sup> + y<sub>2</sub><sup>2</sup> y<sub>2</sub> 1}, {x<sub>3</sub><sup>2</sup> + y<sub>3</sub><sup>2</sup> y<sub>3</sub> 1}}
+
+by = {{x<sub>1</sub><sup>2</sup> + y<sub>1</sub><sup>2</sup> x<sub>1</sub> 1}, {x<sub>2</sub><sup>2</sup> + y<sub>2</sub><sup>2</sup> x<sub>2</sub> 1}, {x<sub>3</sub><sup>2</sup> + y<sub>3</sub><sup>2</sup> x<sub>3</sub> 1}}
+
+xo = -bx / 2 * a
+
+yo = -by / 2 * a
+
+
+[circumcircle]:http://mathworld.wolfram.com/Circumcircle.html
+
+
+

processing_app/library/vecmath/vec2d/library/circle/lib/circumcircle.rb

@@ -1,40 +1,43 @@
-# frozen_string_literal: true
-PBisector = Struct.new(:vector, :angle) # perpendicular bisector
-Vect = Struct.new(:x, :y, :z) # for calculation of center
 # Circumcircle from 3 points
+require 'matrix'
+
 class Circumcircle
-  include Math
   attr_reader :center, :radius, :points
   def initialize(points)
     @points = points
   end
 
   def calculate
-    ab = bisector(points[0], points[1]) # find 2 midpoints
-    bc = bisector(points[1], points[2])
-    @center = circumcenter(ab, bc)
+    @center = Vec2D.new(-(bx / am), -(by / am))
     @radius = center.dist(points[2]) # points[2] = c
   end
 
-  def bisector(a, b)
-    midpoint = (a + b) / 2.0 # middle of ab (or bc)
-    theta = atan2(b.y - a.y, b.x - a.x) # slope of ab (or bc)
-    PBisector.new(midpoint, theta - PI / 2)
+  private
+
+  # Matrix math see matrix_math.md and in detail
+  # http://mathworld.wolfram.com/Circumcircle.html
+
+  def am
+    2 * Matrix[
+      [points[0].x, points[0].y, 1],
+      [points[1].x, points[1].y, 1],
+      [points[2].x, points[2].y, 1]
+    ].determinant
+  end
+
+  def bx
+    -Matrix[
+      [points[0].x * points[0].x + points[0].y * points[0].y, points[0].y, 1],
+      [points[1].x * points[1].x + points[1].y * points[1].y, points[1].y, 1],
+      [points[2].x * points[2].x + points[2].y * points[2].y, points[2].y, 1]
+    ].determinant
   end
 
-  def circumcenter(pb1, pb2)
-    # equation of the first bisector (ax - y =  -b)
-    a0 = tan pb1.angle
-    v0 = pb1.vector
-    a1 = tan pb2.angle
-    v1 = pb2.vector
-    eq0 = Vect.new(a0, -1, -1 * (v0.y - v0.x * a0))
-    eq1 = Vect.new(a1, -1, -1 * (v1.y - v1.x * a1))
-    # calculate x and y coordinates of the circumcenter
-    ox = (eq1.y * eq0.z - eq0.y * eq1.z) /
-         (eq0.x * eq1.y - eq1.x * eq0.y)
-    oy = (eq0.x * eq1.z - eq1.x * eq0.z) /
-         (eq0.x * eq1.y - eq1.x * eq0.y)
-    Vec2D.new(ox, oy)
+  def by
+    Matrix[
+      [points[0].x * points[0].x + points[0].y * points[0].y, points[0].x, 1],
+      [points[1].x * points[1].x + points[1].y * points[1].y, points[1].x, 1],
+      [points[2].x * points[2].x + points[2].y * points[2].y, points[2].x, 1]
+    ].determinant
   end
 end

processing_app/library/vecmath/vec2d/library/simple_circle/circle.rb

@@ -0,0 +1,3 @@
+require_relative 'lib/circumcircle'
+require_relative 'lib/t_points'
+require_relative 'lib/triangle_point'

processing_app/library/vecmath/vec2d/library/simple_circle/lib/circumcircle.rb

@@ -1,40 +1,43 @@
-# frozen_string_literal: true
-PBisector = Struct.new(:vector, :angle) # perpendicular bisector
-Vect = Struct.new(:x, :y, :z) # for calculation of center
 # Circumcircle from 3 points
+require 'matrix'
+
 class Circumcircle
-  include Math
   attr_reader :center, :radius, :points
   def initialize(points)
     @points = points
   end
 
   def calculate
-    ab = bisector(points[0], points[1]) # find 2 midpoints
-    bc = bisector(points[1], points[2])
-    @center = circumcenter(ab, bc)
+    @center = Vec2D.new(-(bx / am), -(by / am))
     @radius = center.dist(points[2]) # points[2] = c
   end
 
-  def bisector(a, b)
-    midpoint = (a + b) / 2.0 # middle of ab (or bc)
-    theta = atan2(b.y - a.y, b.x - a.x) # slope of ab (or bc)
-    PBisector.new(midpoint, theta - PI / 2)
+  private
+
+  # Matrix math see matrix_math.md and in detail
+  # http://mathworld.wolfram.com/Circumcircle.html
+
+  def am
+    2 * Matrix[
+      [points[0].x, points[0].y, 1],
+      [points[1].x, points[1].y, 1],
+      [points[2].x, points[2].y, 1]
+    ].determinant
+  end
+
+  def bx
+    -Matrix[
+      [points[0].x * points[0].x + points[0].y * points[0].y, points[0].y, 1],
+      [points[1].x * points[1].x + points[1].y * points[1].y, points[1].y, 1],
+      [points[2].x * points[2].x + points[2].y * points[2].y, points[2].y, 1]
+    ].determinant
   end
 
-  def circumcenter(pb1, pb2)
-    # equation of the first bisector (ax - y =  -b)
-    a0 = tan pb1.angle
-    v0 = pb1.vector
-    a1 = tan pb2.angle
-    v1 = pb2.vector
-    eq0 = Vect.new(a0, -1, -1 * (v0.y - v0.x * a0))
-    eq1 = Vect.new(a1, -1, -1 * (v1.y - v1.x * a1))
-    # calculate x and y coordinates of the circumcenter
-    ox = (eq1.y * eq0.z - eq0.y * eq1.z) /
-         (eq0.x * eq1.y - eq1.x * eq0.y)
-    oy = (eq0.x * eq1.z - eq1.x * eq0.z) /
-         (eq0.x * eq1.y - eq1.x * eq0.y)
-    Vec2D.new(ox, oy)
+  def by
+    Matrix[
+      [points[0].x * points[0].x + points[0].y * points[0].y, points[0].x, 1],
+      [points[1].x * points[1].x + points[1].y * points[1].y, points[1].x, 1],
+      [points[2].x * points[2].x + points[2].y * points[2].y, points[2].x, 1]
+    ].determinant
   end
 end

processing_app/library/vecmath/vec2d/library/simple_circle/lib/particle.rb

@@ -0,0 +1,113 @@
+# The Particle System
+
+class ParticleSystem
+  include Processing::Proxy
+
+  attr_reader :particles, :particle_shape
+
+  def initialize(width, height, sprite, n)
+    @particles = []
+    # The PShape is a group
+    @particle_shape = create_shape(GROUP)
+
+    # Make all the Particles
+    n.times do |i|
+      particles << Particle.new(width, height, sprite)
+      # Each particle's PShape gets added to the System PShape
+      particle_shape.add_child(particles[i].s_shape)
+    end
+  end
+
+  def update
+    particles.each do |p|
+      p.update
+    end
+  end
+
+  def set_emitter(x, y)
+    particles.each do |p|
+      # Each particle gets reborn at the emitter location
+      p.rebirth(x, y) if (p.dead?)
+    end
+  end
+
+  def display
+    shape(particle_shape)
+  end
+end
+
+# An individual Particle
+
+class Particle
+  include Processing::Proxy
+
+  GRAVITY = Vec2D.new(0, 0.1)
+
+  attr_reader :center, :velocity, :lifespan, :s_shape, :part_size,
+  :boundary_x, :boundary_y, :sprite
+
+  def initialize width, height, sprite
+    @sprite = sprite
+    @boundary_x = Boundary.new(0, width)
+    @boundary_y = Boundary.new(0, height)
+    part_size = rand(10..60)
+    # The particle is a textured quad
+    @s_shape = create_shape
+    s_shape.begin_shape(QUAD)
+    s_shape.no_stroke
+    s_shape.texture(sprite)
+    s_shape.normal(0, 0, 1)
+    s_shape.vertex(-part_size / 2.0, -part_size / 2.0, 0, 0)
+    s_shape.vertex(+part_size / 2.0, -part_size / 2.0, sprite.width, 0)
+    s_shape.vertex(+part_size / 2.0, +part_size / 2.0, sprite.width, sprite.height)
+    s_shape.vertex(-part_size / 2.0, +part_size / 2.0, 0, sprite.height)
+    s_shape.end_shape
+    # Initialize center vector
+    @center = Vec2D.new
+    # Set the particle starting location
+    rebirth(width / 2.0, height / 2.0)
+  end
+
+  def rebirth(x, y)
+    theta = rand(-PI .. PI)
+    speed = rand(0.5 .. 4)
+    # A velocity with random angle and magnitude
+    @velocity = Vec2D.from_angle(theta)
+    @velocity *= speed
+    # Set lifespan
+    @lifespan = 255
+    # Set location using translate
+    s_shape.reset_matrix
+    s_shape.translate(x, y)
+    # Update center vector
+    center.x, center.y = x, y
+  end
+
+  # Is it off the screen, or its lifespan is over?
+  def dead?
+    return true if lifespan < 0
+    return true if boundary_y.exclude? center.y
+    return true if boundary_x.exclude? center.x
+    false
+  end
+
+  def update
+    # Decrease life
+    @lifespan = lifespan - 1
+    # Apply gravity
+    @velocity += GRAVITY
+    s_shape.set_tint(color(255, lifespan))
+    # Move the particle according to its velocity
+    s_shape.translate(velocity.x, velocity.y)
+    # and also update the center location
+    @center += velocity
+  end
+end
+
+# unusually in this case we are looking for excluded values
+
+Boundary = Struct.new(:lower, :upper) do
+  def exclude? val
+    true unless (lower...upper).cover? val
+  end
+end

processing_app/library/vecmath/vec2d/library/simple_circle/lib/t_points.rb

@@ -0,0 +1,33 @@
+# frozen_string_literal: true
+require 'forwardable'
+MAX_POINT = 3
+# A collection of a maximum of 3 points in the processing world
+# includes a collinearity test using Vec2D
+class TPoints
+  extend Forwardable
+  def_delegators(:@points, :each, :map, :size, :shift, :clear, :[])
+  include Enumerable
+
+  attr_reader :points
+
+  def initialize
+    @points = []
+  end
+
+  def <<(pt)
+    points << pt
+    shift if size > MAX_POINT
+  end
+
+  def collinear?
+    full? ? (vec[0] - vec[1]).cross(vec[1] - vec[2]).zero? : false
+  end
+
+  def vec
+    points.map { |point| point.pos }
+  end
+
+  def full?
+    points.length == MAX_POINT
+  end
+end

processing_app/library/vecmath/vec2d/library/simple_circle/lib/triangle_point.rb

@@ -0,0 +1,45 @@
+# frozen_string_literal: true
+# particle and triangle point
+class TPoint
+  include Processing::Proxy
+  attr_reader :pos, :vel, :accel, :xbound, :ybound
+
+  def initialize(position)
+    @pos = position
+    @vel = Vec2D.new
+    @accel = Vec2D.random
+    @xbound = Boundary.new(0, width)
+    @ybound = Boundary.new(0, height)
+  end
+
+  def direction(angle)
+    # direction of the acceleration is defined by the new angle
+    @accel = Vec2D.from_angle(angle)
+    # magnitude of the acceleration is proportional to the angle between
+    # acceleration and velocity
+    dif = accel.angle_between(vel)
+    dif = map1d(dif, 0..PI, 0.1..0.001)
+    @accel *= dif
+  end
+
+  def update
+    @vel += accel
+    @vel.set_mag(1.5) { vel.mag > 1.5 }
+    @pos += vel
+    check_bounds
+  end
+
+  private
+
+  def check_bounds
+    @vel.x *= -1 if xbound.exclude? pos.x
+    @vel.y *= -1 if ybound.exclude? pos.y
+  end
+end
+
+# we are looking for excluded values
+Boundary = Struct.new(:lower, :upper) do
+  def exclude?(val)
+    true unless (lower...upper).cover? val
+  end
+end