Fix and add tests for angleToNb mapping (follow-up to #207)

Cell2Fire/Cells.cpp

@@ -137,34 +137,7 @@ Cells::initializeFireFields(std::vector<std::vector<int>>& coordCells,
         {
             int a = -1 * coordCells[nb - 1][0] + coordCells[this->id][0];
             int b = -1 * coordCells[nb - 1][1] + coordCells[this->id][1];
-
-            int angle = -1;
-            if (a == 0)
-            {
-                if (b >= 0)
-                    angle = 270;
-                else
-                    angle = 90;
-            }
-            else if (b == 0)
-            {
-                if (a >= 0)
-                    angle = 180;
-                else
-                    angle = 0;
-            }
-            else
-            {
-                // TODO: check this logi
-                double radToDeg = 180 / M_PI;
-                // TODO: i think all the negatives and abs cancel out
-                double temp = std::atan(b * 1.0 / a) * radToDeg;
-                if (a > 0)
-                    temp += 180;
-                if (a < 0 && b > 0)
-                    temp += 360;
-                angle = temp;
-            }
+            double angle = std::fmod(std::atan2(-b, -a) * 180.0 / M_PI + 360.0, 360.0);
             this->angleDict[nb] = angle;
             if (availSet.find(nb) != availSet.end())
             {
@@ -442,14 +415,24 @@ Cells::manageFire(int period,
     msg_list_aux.push_back(0);
     std::vector<int> msg_list;
 
+    /* Matias V. proposed change:
     int head_angle = wdf_ptr->waz - 90;
     if (head_angle < 0)
         head_angle += 360;
-
-    head_angle = std::round(head_angle / 45.0) * 45.0;
+    New logic: Convert meteorological azimuth to grid angle
+        angleToNb uses: 0°=South, 90°=East, 180°=North, 270°=West
+        Meteorological: 0°=North, 90°=East, 180°=South, 270°=West
+        Conversion: grid = (180 - met + 360) mod 360
+     */
+    int head_angle = (static_cast<int>(180.0 - wdf_ptr->waz) + 360) % 360;
+
+    // Snap to nearest 45° increment for 8-direction lookup
+    head_angle = static_cast<int>(std::round(head_angle / 45.0) * 45.0);
+    if (head_angle >= 360)
+        head_angle = 0;
 
     int head_cell = angleToNb[head_angle];  // head cell for slope calculation
-    if (head_cell <= 0)                       // solve boundaries case
+    if (head_cell <= 0)                     // solve boundaries case
     {
         head_cell = this->realId;  // as it is used only for slope calculation, if
                                    // it is a boundary cell, it uses the

Cell2Fire/makefile

@@ -25,7 +25,7 @@ TEST_TARGET = test_cell2fire
 
 # Test source files
 TEST_DIR = ../test/unit_tests
-TEST_FILES = $(TEST_DIR)/*.cpp FuelModelKitral.cpp FuelModelUtils.cpp
+TEST_FILES = $(TEST_DIR)/*.cpp FuelModel*.cpp Cells.cpp
 
 # Build rules
 all: $(TARGET)

Cell2Fire/makefile.test_angleToNb

@@ -0,0 +1,27 @@
+# Makefile for angleToNb unit test
+# Usage: make -f makefile.test_angleToNb
+
+CXX = g++
+CXXFLAGS = -m64 -fPIC -fno-strict-aliasing -fexceptions -fopenmp -DNDEBUG -DIL_STD -std=c++14 -O3
+CATCH2_V3_LIBS := $(shell echo "int main(){}" | $(CXX) -x c++ - -lCatch2Main -lCatch2 -o /dev/null 2>/dev/null && echo "-lCatch2Main -lCatch2" || echo "")
+LDFLAGS = -lm -lpthread $(CATCH2_V3_LIBS)
+
+# Source files needed for the test
+TEST_SRC = ../test/unit_tests/test_angleToNb.cpp
+SRCS = Cells.cpp FuelModel*.cpp
+
+# Target executable
+TARGET = test_angleToNb
+
+all: $(TARGET)
+
+$(TARGET): $(TEST_SRC) $(SRCS)
+	$(CXX) $(CXXFLAGS) -o $(TARGET) $(TEST_SRC) $(SRCS) $(LDFLAGS)
+
+run: $(TARGET)
+	./$(TARGET)
+
+clean:
+	rm -f $(TARGET)
+
+.PHONY: all run clean

test/unit_tests/README_angleToNb.md

@@ -0,0 +1,113 @@
+# angleToNb Unit Test
+
+## Purpose
+
+This test validates the `angleToNb` mapping and wind direction conversion logic that fixes issue #207: "slope effect not working on kitral due to head cell calculation".
+
+## What It Tests
+
+1. **angleToNb mapping**: Verifies that each of the 8 cardinal and diagonal directions correctly maps to the appropriate neighbor cell
+2. **Meteorological to grid angle conversion**: Tests the conversion from meteorological azimuth (0° = North) to grid coordinates (0° = East)
+3. **16 wind directions**: Validates that all 16 compass directions (every 22.5°) correctly round to the nearest 45° increment
+4. **Head cell selection**: Simulates the actual use case of selecting the correct neighbor for slope calculations based on wind direction
+5. **Boundary handling**: Tests edge and corner cells that don't have all 8 neighbors
+6. **Angle snapping**: Verifies the rounding mechanism that snaps arbitrary angles to 45° increments
+
+## The Bug and The Fix
+
+### Original Bug
+The original code directly used the meteorological wind azimuth as an index:
+```cpp
+int head_cell = angleToNb[wdf_ptr->waz];  // WRONG!
+```
+
+This failed because:
+- **Meteorological convention**: 0° = North, 90° = East, 180° = South, 270° = West
+- **Grid convention**: 0° = East, 90° = North, 180° = West, 270° = South
+
+### The Fix
+```cpp
+int head_angle = wdf_ptr->waz - 90;        // Convert meteorological → grid
+if (head_angle < 0)
+    head_angle += 360;                      // Keep in [0, 360) range
+
+head_angle = std::round(head_angle / 45.0) * 45.0;  // Snap to nearest 45°
+
+int head_cell = angleToNb[head_angle];      // Correct lookup
+```
+
+## Building and Running
+
+### Option 1: Using the dedicated makefile
+```bash
+cd Cell2Fire
+make -f makefile.test_angleToNb
+./test_angleToNb
+```
+
+### Option 2: Using the main makefile
+```bash
+cd Cell2Fire
+make tests
+./test_cell2fire "[angleToNb]"  # Run only angleToNb tests
+```
+
+## Test Output
+
+You should see output like:
+```
+All tests passed (XX assertions in X test cases)
+```
+
+## Test Cases
+
+### 1. Basic 8-direction mapping
+Tests that the 8 cardinal and diagonal directions map correctly to neighbors.
+
+### 2. Meteorological to grid conversion
+Tests conversion for all cardinal and diagonal directions.
+
+### 3. 16 wind directions
+Tests all 16 compass directions to ensure proper rounding:
+- N, NNE, NE, ENE, E, ESE, SE, SSE, S, SSW, SW, WSW, W, WNW, NW, NNW
+
+### 4. Head cell selection with 16 winds
+Simulates actual usage with wind from all 16 directions, verifying the correct neighbor is selected for slope calculations.
+
+### 5. Boundary cells
+Tests corner and edge cells that have fewer than 8 neighbors.
+
+### 6. Angle snapping
+Validates the mathematical rounding to 45° increments.
+
+## Grid Layout
+
+The test uses a 5×5 grid with the center cell at position 13:
+```
+ 1  2  3  4  5
+ 6  7  8  9 10
+11 12 13 14 15
+16 17 18 19 20
+21 22 23 24 25
+```
+
+Neighbors of cell 13:
+```
+NW(7)  N(8)   NE(9)
+W(12)  C(13)  E(14)
+SW(17) S(18)  SE(19)
+```
+
+## Dependencies
+
+- Catch2 testing framework
+- C++14 or later
+- Standard math library
+
+## Related Files
+
+- Source: `Cell2Fire/Cells.cpp` (lines 124-180: `initializeFireFields()`)
+- Source: `Cell2Fire/Cells.cpp` (line 451: fixed `manageFire()`)
+- Header: `Cell2Fire/Cells.h` (line 102: `angleToNb` declaration)
+- Issue: #207
+- Pull Request: #208

test/unit_tests/nb_debug.py

@@ -0,0 +1,80 @@
+#!python3
+import numpy as np
+from matplotlib.pyplot import plot, show
+
+# testing head_cell functions
+
+
+@np.vectorize
+def func1(x):
+    return (int(180 - x) + 360) % 360
+
+
+@np.vectorize
+def func2(x):
+    a = x - 90
+    if a > 0:
+        a += 360
+    return a
+
+
+circle = np.linspace(0, 360, num=3600, endpoint=False)
+
+# for c in circle:
+#     print(f"func({c}) = {func(c)}")
+
+plot(circle, func1(circle))
+plot(circle, func2(circle))
+show()
+
+# testing nb to angle
+
+
+def nb_angle1(base, neighbor):
+    """original code"""
+    a, b = base - neighbor
+    if a == 0:
+        if b >= 0:
+            return 270
+        else:
+            return 90
+    elif b == 0:
+        if a >= 0:
+            return 180
+        else:
+            return 0
+    # print("1")
+    return nb_angle2(base, neighbor)
+
+
+def nb_angle2(base, neighbor):
+    """only the second part of original code"""
+    a, b = base - neighbor
+    tmp = np.atan(b / a) * 180 / np.pi
+    # print(a,b,tmp)
+    if a > 0:
+        # print("a>0 +180")
+        # print("2")
+        tmp += 180
+    if a < 0 and b > 0:
+        # print("a<0 b>0 +360")
+        # print("3")
+        tmp += 360
+    return tmp
+
+
+def nb_angle3(base, neighbor):
+    """suggested fix"""
+    a, b = base - neighbor
+    return (np.arctan2(-b, -a) * 180 / np.pi + 360) % 360
+
+
+base = np.array([2, 2])
+neighbors = np.meshgrid(np.arange(5), np.arange(5))
+
+print("base, neighbor, angle1, angle2, angle3".upper())
+for nb in zip(neighbors[0].flatten(), neighbors[1].flatten()):
+    angle1 = nb_angle1(base, np.array(nb))
+    angle2 = nb_angle2(base, np.array(nb))
+    angle3 = nb_angle3(base, np.array(nb))
+    print(f"base: {base} neighbor: {list(map(int,nb))} angle: {angle1:6.2f}, {angle2:6.2f}, {angle3:6.2f}")