Fix PHOT reflecting off thin walls of particles incorrectly

When PHOT fails to move (do_move or eval_move return "no move"), it looks for
a surface (a contour of boundaries, as reported by is_boundary) along its path
and reflects off (or refracts into, see below) it, using get_normal_interp to
find the point of incidence and get_normal to deduce the surface normal.
get_normal is given the point and angle of incidence, and attempts to traverse
the surface the point belongs to by running two "surface scout" processes.

These processes remember their own position and "heading", a subset of the
eight cardinal directions on the grid. They are initialized with the point of
incidence and a heading that includes all directions whose dot product with
the angle of incidence is non-negative (see direction_to_map). They then
perform a few iterations (SURF_RANGE).

In each iteration, the processes check all eight neighbours of the cell they
are on and select the first neighbouring cell they find that is both a
boundary (as reported by is_boundary) and that is within their heading. They
then move to this neighbouring cell and update their heading by discarding
directions that are not similar enough to (differ by more than 45 degrees
from) the one that took them where they are now (see find_next_boundary). If
they find no such neighbour, they stop.

Continuing the militaristic line of thinking introduced by the term "surface
scout", you can imagine the two processes as two paratroopers who arrive from
above, land on a horizontal surface, and one starts going left, while the
other starts going right. They initially expect the surface they land on to be
close to horizontal, but are also prepared for not too erratic changes in its
angle as they go. Changes too erratic (imagine a precipice) scare them and
force them to stop.

Once the processes finish, an imaginary line segment is drawn between the
cells they ended up on. If the line segment is long enough (estimated by j,
and compared against NORMAL_MIN_EST), get_normal returns a normal that is
perpendicular to it. If it is too short, get_normal gives up and returns
nothing (which results in the PHOT being killed).

This amounts to our paratroopers attempting to get the "lay of the land" by
walking away from where they landed and comparing where they end up. They also
know that if they are still relatively close to each other at the end of their
walk, their measurement is probably wrong and their mission should be aborted.

The bug this commit fixes is that get_normal returns bogus surface normals
when it encounters thin walls of particles, defined as walls exactly two
layers of particles thick. One-layer walls are not really walls, as movement
code allows particles to penetrate these, and three-layer and thicker walls
are too thick for the bug to manifest.

The bug manifests for two-layer walls because the "left" scout process is
drawn to the side of the wall opposite to the one with the point of incidence.
This is because scout processes check neighbours in a clockwise order, and
always select the first suitable neighbour they find. As particles on the
other side of the wall are both boundaries and are within the heading of the
processes, they also qualify as suitable neighbours, so whether a scout
process selects the correct side of the wall depends on the order in which
neighbours are checked.

Essentially, the paratroopers look at their immediate surroundings in a
clockwise order. The right paratrooper always finds the ground and knows where
to step. The left paratrooper finds the Upside Down from Stranger Things and
teleports there.

This bug also affects refraction into and out of thin walls, but since these
walls are thin, the path the PHOT takes inside them is rather short and the
incorrect angle of travel is difficult to see. Furthermore, upon exit, the
same normal deduction bug causes the PHOT to take a path whose angle is almost
identical to that of the path that took it to the wall, so much so that it is
also difficult to see over shorter distances.

The solution is to have the left scout process check neighbours in reverse
order, so that it prefers the right side of the wall over the wrong one. This
does not affect its behaviour when facing thicker walls, but fixes its
behaviour when facing two-layer walls.

The changes in this commit also make find_next_boundary interact with
is_blocking directly to detect a change between the blocking trait of
immediate neighbours. This makes more sense than relying on is_boundary
because find_next_boundary is meant to find a transition from non-blocking to
blocking neighbours within the current heading, rather than to find any
boundary particle. The difference is subtle but important.

src/simulation/Simulation.cpp

@@ -2987,25 +2987,35 @@ int Simulation::is_boundary(int pt, int x, int y)
 	return 1;
 }
 
-int Simulation::find_next_boundary(int pt, int *x, int *y, int dm, int *em)
+int Simulation::find_next_boundary(int pt, int *x, int *y, int dm, int *em, bool reverse)
 {
 	static int dx[8] = {1,1,0,-1,-1,-1,0,1};
 	static int dy[8] = {0,1,1,1,0,-1,-1,-1};
 	static int de[8] = {0x83,0x07,0x0E,0x1C,0x38,0x70,0xE0,0xC1};
-	int i, ii, i0;
 
 	if (*x <= 0 || *x >= XRES-1 || *y <= 0 || *y >= YRES-1)
+	{
 		return 0;
+	}
 
-	if (*em != -1) {
-		i0 = *em;
-		dm &= de[i0];
-	} else
-		i0 = 0;
+	if (*em != -1)
+	{
+		dm &= de[*em];
+	}
 
-	for (ii=0; ii<8; ii++) {
-		i = (ii + i0) & 7;
-		if ((dm & (1 << i)) && is_boundary(pt, *x+dx[i], *y+dy[i])) {
+	unsigned int mask = 0;
+	for (int i = 0; i < 8; ++i)
+	{
+		if ((dm & (1U << i)) && is_blocking(pt, *x + dx[i], *y + dy[i]))
+		{
+			mask |= (1U << i);
+		}
+	}
+	for (int i = 0; i < 8; ++i)
+	{
+		int n = (i + (reverse ? 1 : -1)) & 7;
+		if (((mask & (1U << i))) && !(mask & (1U << n)))
+		{
 			*x += dx[i];
 			*y += dy[i];
 			*em = i;
@@ -3039,9 +3049,9 @@ int Simulation::get_normal(int pt, int x, int y, float dx, float dy, float *nx,
 	j = 0;
 	for (i=0; i<SURF_RANGE; i++) {
 		if (lv)
-			lv = find_next_boundary(pt, &lx, &ly, ldm, &lm);
+			lv = find_next_boundary(pt, &lx, &ly, ldm, &lm, true);
 		if (rv)
-			rv = find_next_boundary(pt, &rx, &ry, rdm, &rm);
+			rv = find_next_boundary(pt, &rx, &ry, rdm, &rm, false);
 		j += lv + rv;
 		if (!lv && !rv)
 			break;
@@ -3052,7 +3062,6 @@ int Simulation::get_normal(int pt, int x, int y, float dx, float dy, float *nx,
 
 	if ((lx == rx) && (ly == ry))
 		return 0;
-
 	ex = float(rx - lx);
 	ey = float(ry - ly);
 	r = 1.0f/hypot(ex, ey);

src/simulation/Simulation.h

@@ -129,7 +129,7 @@ public:
 
 	int is_blocking(int t, int x, int y);
 	int is_boundary(int pt, int x, int y);
-	int find_next_boundary(int pt, int *x, int *y, int dm, int *em);
+	int find_next_boundary(int pt, int *x, int *y, int dm, int *em, bool reverse);
 	void photoelectric_effect(int nx, int ny);
 	unsigned direction_to_map(float dx, float dy, int t);
 	int do_move(int i, int x, int y, float nxf, float nyf);