farsite: match the reference at flank and back

CHANGELOG.md

@@ -39,6 +39,9 @@ history, so they summarise each release rather than list every change. The
 
 ### Fixed
 
+- The Farsite model now matches the FARSITE reference (mbedward/farsite) at
+  the flanks and back of the fire, not only at the head. On slopes it is an
+  approximation. This changes Farsite results. ([#270])
 - The Farsite model reads the wind in m/s and the slope as a tangent, like the
   other Rothermel models, and slope only speeds up spread uphill. It used to
   read the wind as if it were in mph and double the slope, so downhill spread
@@ -343,3 +346,4 @@ repository stays easy to return to.
 [#172]: https://github.com/forefireAPI/forefire/issues/172
 [#265]: https://github.com/forefireAPI/forefire/issues/265
 [#242]: https://github.com/forefireAPI/forefire/pull/242
+[#270]: https://github.com/forefireAPI/forefire/pull/270

src/propagation/FarsitePropagationModel.cpp

@@ -29,6 +29,8 @@ namespace libforefire
 
 		size_t idx_slope;
 		size_t idx_wind;
+		size_t idx_windU;
+		size_t idx_windV;
 
 		size_t idx_m_ones;
 		size_t idx_m_liveh;
@@ -88,6 +90,8 @@ namespace libforefire
 		/* defining the properties needed for the model */
 		idx_slope = registerProperty("slope");
 		idx_wind = registerProperty("normalWind");
+		idx_windU = registerProperty("windU");
+		idx_windV = registerProperty("windV");
 
 		idx_m_ones = registerProperty("moist.ones");
 		idx_m_liveh = registerProperty("moist.liveh");
@@ -183,7 +187,7 @@ namespace libforefire
 		//-----------------------------------------------------------------------------
 		// 2) Read actual values from valueOf
 		//-----------------------------------------------------------------------------
-		double windspd = adjustementWind * valueOf[idx_wind]; // wind times cos(slope) here, m/s
+		double windspd = adjustementWind * valueOf[idx_wind]; // along the normal, times cos(slope), m/s
 		double slope = valueOf[idx_slope];
 		// Except with splines, ForeFire puts the slope in the normal's z, so
 		// models get sin(slope) and the wind times cos(slope). FARSITE wants the
@@ -197,6 +201,20 @@ namespace libforefire
 			slope = s / c;
 			windspd /= c;
 		}
+		// FARSITE spreads each point as an ellipse along the wind, so the rate
+		// is computed with the whole wind and then taken at the angle between
+		// the node's normal and the wind. On a slope this is an approximation:
+		// the node only sees the slope along its normal, so the ellipse follows
+		// the wind alone.
+		double windU = valueOf[idx_windU];
+		double windV = valueOf[idx_windV];
+		double windNorm = adjustementWind * std::hypot(windU, windV);
+		double cosWind = 1.0;
+		if (windNorm > 1e-9)
+			cosWind = std::max(-1.0, std::min(1.0, windspd / windNorm));
+		else
+			windNorm = windspd;
+		windspd = windNorm;
 		// moinstures variables
 		double m_ones = valueOf[idx_m_ones];  
 		double m_liveh = valueOf[idx_m_liveh]; 
@@ -492,6 +510,24 @@ namespace libforefire
 			rate_of_spread = rateo * (1 + phiew);
 		}
 
+		// Ellipse: FARSITE's Mechanix::ellipse (Anderson 1983 length-to-breadth,
+		// Alexander 1985 head-to-back), wind in mph. The rate normal to the
+		// front at angle theta from the wind is Richards (1990):
+		// c cos(theta) + sqrt(a^2 cos^2(theta) + b^2 sin^2(theta)).
+		double windMph = std::max(windspd, 0.0) / 0.44704;
+		double lb = 0.936 * std::exp(0.1147 * windMph)
+				  + 0.461 * std::exp(-0.0692 * windMph) - 0.397;
+		lb = std::max(1.0, std::min(8.0, lb));
+		double part = std::sqrt(lb * lb - 1.0);
+		double head = rate_of_spread;
+		double back = head * (lb - part) / (lb + part);
+		double ea = (head + back) / 2.0;
+		double ec = (head - back) / 2.0;
+		double eb = ea / lb;
+		double sin2 = 1.0 - cosWind * cosWind;
+		rate_of_spread = ec * cosWind
+					   + std::sqrt(ea * ea * cosWind * cosWind + eb * eb * sin2);
+
 
 		rate_of_spread = rate_of_spread/60.0;
 

tests/python/test_farsite_reference.py

@@ -23,8 +23,9 @@
           is in how ForeFire derives the normal wind and slope it hands over.
 
 The test fails when the head rate of spread differs beyond REL_TOL in either
-comparison. Flank, back and length-to-breadth are only reported: ForeFire has
-no ellipse, so a node not facing the wind spreads at the no-wind rate.
+comparison, or the flank or back one does on flat ground. On a slope they are
+only reported: the model sees the slope along each node's normal, not the
+whole slope vector FARSITE turns its ellipse with.
 
 Every case runs in its own interpreter, because a second ``ForeFire()`` in a
 process inherits the first one's parameters (see test_moisture_invariants.py).
@@ -267,9 +268,13 @@ def main(argv=None):
             "head(model)": rel(h["ROS"], mo["head"]),
             "head(pipeline)": rel(h["ROS"], pi["head"]),
         }
+        if c["slope"] == 0:
+            errs["flank"] = rel(flank["ROS"], pi["flank"])
+            errs["back"] = rel(back["ROS"], pi["back"])
         # The spline normals do not point at the head: report that gap only.
         checked = [k for k in errs
-                   if not (k == "head(pipeline)" and args.scheme == "splines")]
+                   if not (k in ("head(pipeline)", "flank", "back")
+                           and args.scheme == "splines")]
         bad = [k for k in checked if errs[k] > REL_TOL]
         if bad:
             failures.append((c, bad))