Monday, November 02, 2015


It's hard to find a newspaper without one, and I personally find that Sudoku solutions are an interesting problem. Randomly filling a grid of 81 cells with numbers from 1-9 wouldn't take long, but the end result probably wouldn't be a Sudoku solution, where each row, each column, and each box, all need to have the numbers 1-9 exactly once.

Instead of doing it randomly, we could try a constrained approach: only pick numbers from the set that would form an allowable Sudoku solution. We could start off by defining the set of all values in row M, all values in column N, and all values in box MN. The union of all three sets would indicate the values that have already been picked; therefore they're the numbers we couldn't re-use if we were to make a valid Sudoku solution. \( R(m) \cup C(n) \cup B(m,n) \) is the numbers that have been used, \( (R(m) \cup C(n) \cup B(m,n))' \) are the numbers available to us (this last set is the equivalent of \(R(m)' \cap C(n)' \cap B(m,n)'\)).

Let's assume we're looking at an empty cell, in an empty grid. \(R(m)=\emptyset, C(n)=\emptyset, B(m,n)=\emptyset\), therefore we can choose anything from \(\emptyset'\). The world is our oyster. Let's pick 1. We could have picked anything, but 1's a good place to start. In the adjacent cell, \(R(m)\) and \(B(m,n)\) would both yield \(\{1\}\) so we'd be forced to pick anything from \(\{1\}'\). Let's pick 2. Continuing in ascending order (and increasing level of constraint) we'd finally reach the 9th cell where \(R(m)=\{1,2,3,4,5,6,7,8\}\) and \(B(m,n)=\{7,8\}\). The only number we can pick is 9. Our top row is complete.

Continuing at the first column of the second row, and proceeding in the same fashion, we keep filling cells. It's too easy, something's got to break.

Ah. On the second row, at the seventh colum, we reach Box(2,7). We've already written \(\{4,5,6,1,2,3\}\) into the row, and \(\{7,8,9\}\) into the box. What can we do now?

The answer is to backtrack. We didn't have to choose the lowest value from \(R(m)' \cap C(n)' \cap B(m,n)'\) every time; in fact - more often than not - we had a wide array of options. So, we backtrack. That involves reverting the operation on the current cell (leaving it blank), moving to the previous cell, making sure it's blank (but we remembered which value(s) we'd already tried), and then trying another value. We don't need to limit ourselves to going back just one cell, either. In fact, in this case, we need to revert the last three cells. A second row that looks like \(\{4,5,6,7,8,9,1,2,3\}\) is a valid partial solution to a Sudoku problem. We use backtracking even more on the third row, producing \(\{7,8,9,1,2,3,4,5,6\}\). Still: a valid partial solution.

It turns out that our set-based constraints, plus the backtracking ability, give us the power to solve any Sudoku problem, even when the starting point is ambiguous (in the case of the empty grid).

1 2 3 4 5 6 7 8 9
4 5 6 7 8 9 1 2 3
7 8 9 1 2 3 4 5 6
2 1 4 3 6 5 8 9 7
3 6 5 8 9 7 2 1 4
8 9 7 2 1 4 3 6 5
5 3 1 6 4 2 9 7 8
6 4 2 9 7 8 5 3 1
9 7 8 5 3 1 6 4 2

Saturday, October 03, 2015


var forEach = new TransformManyBlock<IEnumerable<int>, int>(x => x);

Saturday, May 17, 2014


Satellites in a geostationary orbit appear to maintain a fixed position relative to a position on the surface of the Earth. To accomplish this, they need to rotate at the same angular speed that the Earth spins, which makes one full rotation every sidereal day. There are 23.93446122 hours (or 86,164.0604 minutes) in one sidereal day. A full revolution of 2π radians in 86,164.0604 seconds gives us \(\omega = 7.29212 \times 10^{-5} \ radians \cdot s^{-1}\). This is our target angular speed if we want to keep a satellite overhead the same spot.

Earth's standard gravitational parameter is its mass M x G (the gravitational constant), which is roughly μ = 398,600.4418. The standard gravitational parameters for planets are better known than either M or G individually (the best way we have at our disposal to weigh extremely massive objects is by observing their gravitational effect on other bodies like orbiting satellites). Acceleration due to gravity is \(\mu \over r^2\). The equatorial radius of Earth is 6,378.14 km. If we plug in that figure, we'd see acceleration on the surface is the familiar 0.009798285 km/s. 422 km above the surface (in the orbital sphere of the ISS), it's just 0.008619904 km/s (87.97%). The further out we go, the lower the acceleration and the longer the orbital period.

To find out how far away we have to place a geostationary satellite, consider that it will need to be in a circular orbit with an orbital period of one sidereal day. In one second, it will move \(7.29212 \times 10^{-5} radians\); therefore \(r \times (1 - cos(7.29212 \times 10^{-5}))\) will be the distance that it falls towards the Earth in that time. We know from the SUVAT equation \(s = ut + {1 \over 2}at^2\) that in a second, with no initial velocity, that the distance travelled is half the acceleration. We now have an identity \(r (1 - \cos(7.29212 \times 10^{-5})) = {\mu \over r^2}\). Rearrange it and solve: \[r = \sqrt[3]{\mu \over {2 - 2 cos(7.29212 \times 10^{-5})}}\] It turns out that the altitude of 42164.15974 km above the center of the Earth (35786.02274 km above the equatorial surface) is where you'll find geostationary satellites - and they can only exist in that one orbital plane.

Saturday, May 10, 2014

Gravity / ISS

Let's say we're interested in finding how fast the International Space Station has to move in order not to fall back to earth (i.e. to stay in orbit). The earth has an equatorial radius of about 6384km and the satellite is about 422km above the equator when it's overhead on its inclined orbital path. For the sake of simplicity we will ignore the oblate spheroid shape of the earth, and the replace the elliptical orbit with a circular one, 6806km from the center of the earth.

On Earth's surface the acceleration due to gravity is \(9.8 m \cdot s^{-2}\). The further you go away, the lower the force of gravity. 422km above the surface we're told that the acceleration is 89% of surface gravity - it's \(8.722 m \cdot s^{-2}\).

Using the SUVAT equation \(s = ut + {1 \over 2}at^2\) we can see that an object dropped from 422km (i.e. \(g = 8.772 m \cdot s^{-2}\)) would fall 4.361m in the first second.

In the same second, we know that the ISS traces out the circular orbit (i.e. it doesn't crash into the Earth). The angular distance it travels (we could use the unit circle for visual confirmation) is \(\arccos({{r - s} \over r})\) or \(\arccos({{6806000 - 4.361} \over 6806000})\) or \(\arccos({6805995.639 \over 6806000})\) or simply 0.001132041 radians.

If it travels 0.001132041 radians in a second, it will take 5550.316851 seconds to perform a complete revolution of 2π. 5550 seconds is 92½ minutes, a figure that's very close to the published figure (on Wikipedia) and that's quite amazing given the rough estimates we've made to get this far. A circle of radius 6806km has a circumference of 13612km. Divide that circumference by the orbital period in seconds and you get \(7704 m \cdot s^{-1}\). You could also try \(r \sin(\theta)\) or \(6806000 \sin(0.001132041)\) which gives the same result. It's moving pretty quickly, but would have to be even quicker if the orbit was any closer to the surface!

Wednesday, April 16, 2014

Newtonian Reflector

Newtonian reflector telescopes use two mirrors – a figured primary mirror (the objective) which focuses incoming light, and a flat secondary mirror which redirects that light at an angle so that it can be viewed without getting your head stuck in the tube. The secondary mirror is only practically important – a CCD sensor could be placed directly into the path of the light focused by the objective mirror and it would capture the light without requiring an additional reflection. However, the subject of this post is the primary mirror, particularly: it’s shape.

A 2D parabola is defined by the quadratic equation \(y = ax^2 + bx + c\). If we make some simplifying assumptions (such as: its open side faces up; it “rests” on the x-axis; the focus is on the y-axis at height p) then our equation reduces to \(y = 4px^2\). The first order derivative (e.g. slope) of this equation is \({dy \over dx} = 2px\). That means, for a given focal length (and radius from the centre of the mirror), we can estimate the height of the reflective surface from the x-axis, as well as the slope at that point. These two computed values show where light is reflected when it encounters the surface of the objective mirror. Remember that incident light is reflected about the normal vector to the surface:
\(R = D – (2D \cdot N)N\)

To enable a 3D version (a paraboloid) we can take advantage of the knowledge that we only need rotate the parabola (i.e. the shape is rotationally symmetric about the vertical axis). In 3D, I chose to rotate around the z-axis.

It then becomes rather straight forward to take a solid angle of light, to compute the angles at which each incident ray encounters the parabolic reflector surface (parallax will affect closer light sources more than further away sources) and then note the points at which two or more rays converge to a single point. It turns out that the parabolic shape is ONLY able to focus light that is travelling parallel to the reflector’s primary axis. Closer sources converge further back than the stated focal length; sources at infinity converge exactly at the parabola’s focal point. Off-axis sources result in a complex out-of-focus shape when we attempt to capture them on a flat focal plane like a CCD. Looking at the distribution of focal points for a given field of view it was difficult to imagine a single transformation that might enable coma to be completely (and correctly) removed, though clearly a retail market for coma-correcting lenses abounds.

Thursday, January 24, 2013


Ever wonder what makes those striking day-glo colors shine brighter than everything around them? Our eyes can only see visible light (the colors of the rainbow) but the sun emits electromagnetic waves across a vastly wider spectrum (which we cannot see with our eyes, of course, and not just because it's not a good idea to stare at the sun). But just because we cannot see them, doesn't mean they aren't there. When photons from the ultraviolet region of the spectrum bump into an object, some of the energy can be reflected. Because we can't see them, they go unnoticed. However, certain chemical combinations transform higher-energy ultraviolet photons into lower-energy visible light photons (a process called fluorescence). When those fluorescent objects start emitting more visible light than the non-fluorescent objects around them, they give the perception of being brighter: day-glo. Interesting fact: whiter-than-white washing detergent is known to employ fluorescence to make clothes appear brighter.

Wednesday, October 17, 2012


If you're reading this then it's already too late you've probably come across scalars and vectors and pondered their differences. The pair of terms is used in many contexts as diverse as operations performed by a CPU. I think the names go back as far as matrix mathematics, and - if you'll listen - I'll tell you why...

Matrices can be "multiplied" in two ways:

  1. two matrices are multiplied and the result is an arrangement of dot products;
  2. one matrix and one scalar are multiplied and every value in the matrix is scaled by the scalar.

Now, every time an algorithm requires to scale a collection of values by some number, will I name that variable "factor" or "scalar"?

Which would you choose?