feat: expose XFOIL Type 2/3 polar modes (variable Re)

crates/rustfoil-python/src/lib.rs

@@ -3,7 +3,7 @@ use pyo3::types::PyDict;
 use rayon::prelude::*;
 use rustfoil_core::{naca, point, Body, CubicSpline, PanelingParams, Point};
 use rustfoil_xfoil::oper::{solve_body_oper_point, AlphaSpec};
-use rustfoil_xfoil::XfoilOptions;
+use rustfoil_xfoil::{ReType, XfoilOptions};
 
 struct FaithfulResult {
     alpha_deg: f64,
@@ -17,10 +17,19 @@ struct FaithfulResult {
     x_tr_lower: f64,
     cd_friction: f64,
     cd_pressure: f64,
+    reynolds_eff: f64,
     success: bool,
     error: Option<String>,
 }
 
+fn re_type_from_int(v: u8) -> ReType {
+    match v {
+        2 => ReType::Type2,
+        3 => ReType::Type3,
+        _ => ReType::Type1,
+    }
+}
+
 fn points_from_flat(coords: &[f64]) -> Vec<Point> {
     coords.chunks(2).map(|c| point(c[0], c[1])).collect()
 }
@@ -34,7 +43,7 @@ fn flat_from_points(pts: &[Point]) -> Vec<(f64, f64)> {
 /// Returns a dict with keys: cl, cd, cm, converged, iterations, residual,
 /// x_tr_upper, x_tr_lower, cd_friction, cd_pressure, alpha_deg, success, error.
 #[pyfunction]
-#[pyo3(signature = (coords, alpha_deg, reynolds=1.0e6, mach=0.0, ncrit=9.0, max_iterations=100))]
+#[pyo3(signature = (coords, alpha_deg, reynolds=1.0e6, mach=0.0, ncrit=9.0, max_iterations=100, re_type=1))]
 fn analyze_faithful(
     py: Python<'_>,
     coords: Vec<f64>,
@@ -43,6 +52,7 @@ fn analyze_faithful(
     mach: f64,
     ncrit: f64,
     max_iterations: usize,
+    re_type: u8,
 ) -> PyResult<Py<PyDict>> {
     if coords.len() < 6 || coords.len() % 2 != 0 {
         let d = PyDict::new(py);
@@ -67,6 +77,7 @@ fn analyze_faithful(
         mach,
         ncrit,
         max_iterations,
+        re_type: re_type_from_int(re_type),
         ..Default::default()
     };
 
@@ -84,6 +95,7 @@ fn analyze_faithful(
             d.set_item("cd_friction", r.cd_friction)?;
             d.set_item("cd_pressure", r.cd_pressure)?;
             d.set_item("alpha_deg", r.alpha_deg)?;
+            d.set_item("reynolds_eff", r.reynolds_eff)?;
             d.set_item("success", true)?;
             d.set_item("error", py.None())?;
         }
@@ -242,13 +254,15 @@ fn solve_one_faithful(body: &Body, alpha_deg: f64, options: &XfoilOptions) -> Fa
             converged: r.converged, iterations: r.iterations, residual: r.residual,
             x_tr_upper: r.x_tr_upper, x_tr_lower: r.x_tr_lower,
             cd_friction: r.cd_friction, cd_pressure: r.cd_pressure,
+            reynolds_eff: r.reynolds_eff,
             success: true, error: None,
         },
         Err(e) => FaithfulResult {
             alpha_deg, cl: 0.0, cd: 0.0, cm: 0.0,
             converged: false, iterations: 0, residual: 0.0,
             x_tr_upper: 1.0, x_tr_lower: 1.0,
             cd_friction: 0.0, cd_pressure: 0.0,
+            reynolds_eff: options.reynolds,
             success: false, error: Some(format!("{e}")),
         },
     }
@@ -267,6 +281,7 @@ fn faithful_result_to_pydict(py: Python<'_>, r: &FaithfulResult) -> PyResult<Py<
     d.set_item("x_tr_lower", r.x_tr_lower)?;
     d.set_item("cd_friction", r.cd_friction)?;
     d.set_item("cd_pressure", r.cd_pressure)?;
+    d.set_item("reynolds_eff", r.reynolds_eff)?;
     d.set_item("success", r.success)?;
     d.set_item("error", r.error.as_deref())?;
     Ok(d.into())
@@ -276,7 +291,7 @@ fn faithful_result_to_pydict(py: Python<'_>, r: &FaithfulResult) -> PyResult<Py<
 ///
 /// Returns a list of dicts (same schema as analyze_faithful), one per alpha.
 #[pyfunction]
-#[pyo3(signature = (coords, alphas, reynolds=1.0e6, mach=0.0, ncrit=9.0, max_iterations=100))]
+#[pyo3(signature = (coords, alphas, reynolds=1.0e6, mach=0.0, ncrit=9.0, max_iterations=100, re_type=1))]
 fn analyze_faithful_batch(
     py: Python<'_>,
     coords: Vec<f64>,
@@ -285,6 +300,7 @@ fn analyze_faithful_batch(
     mach: f64,
     ncrit: f64,
     max_iterations: usize,
+    re_type: u8,
 ) -> PyResult<Vec<Py<PyDict>>> {
     let err_msg = if coords.len() < 6 || coords.len() % 2 != 0 {
         Some("Invalid coordinates".to_string())
@@ -298,6 +314,7 @@ fn analyze_faithful_batch(
                 converged: false, iterations: 0, residual: 0.0,
                 x_tr_upper: 1.0, x_tr_lower: 1.0,
                 cd_friction: 0.0, cd_pressure: 0.0,
+                reynolds_eff: reynolds,
                 success: false, error: Some(msg.clone()),
             };
             faithful_result_to_pydict(py, &r)
@@ -315,6 +332,7 @@ fn analyze_faithful_batch(
                     converged: false, iterations: 0, residual: 0.0,
                     x_tr_upper: 1.0, x_tr_lower: 1.0,
                     cd_friction: 0.0, cd_pressure: 0.0,
+                    reynolds_eff: reynolds,
                     success: false, error: Some(msg.clone()),
                 };
                 faithful_result_to_pydict(py, &r)
@@ -324,6 +342,7 @@ fn analyze_faithful_batch(
 
     let options = XfoilOptions {
         reynolds, mach, ncrit, max_iterations,
+        re_type: re_type_from_int(re_type),
         ..Default::default()
     };
 
@@ -421,7 +440,7 @@ fn analyze_inviscid_batch(
 /// ue_upper, ue_lower, x_tr_upper, x_tr_lower, converged, iterations,
 /// residual, success, error.
 #[pyfunction]
-#[pyo3(signature = (coords, alpha_deg, reynolds=1.0e6, mach=0.0, ncrit=9.0, max_iterations=100))]
+#[pyo3(signature = (coords, alpha_deg, reynolds=1.0e6, mach=0.0, ncrit=9.0, max_iterations=100, re_type=1))]
 fn get_bl_distribution(
     py: Python<'_>,
     coords: Vec<f64>,
@@ -430,6 +449,7 @@ fn get_bl_distribution(
     mach: f64,
     ncrit: f64,
     max_iterations: usize,
+    re_type: u8,
 ) -> PyResult<Py<PyDict>> {
     use rustfoil_xfoil::oper::{build_state_from_coords, solve_operating_point_from_state, coords_from_body, AlphaSpec};
     use rustfoil_xfoil::XfoilOptions;
@@ -455,6 +475,7 @@ fn get_bl_distribution(
     let body_coords = coords_from_body(&body);
     let options = XfoilOptions {
         reynolds, mach, ncrit, max_iterations,
+        re_type: re_type_from_int(re_type),
         ..Default::default()
     };
 

crates/rustfoil-xfoil/src/config.rs

@@ -1,3 +1,19 @@
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum ReType {
+    /// Type 1: constant Re (default). `Re_eff = Re`.
+    Type1,
+    /// Type 2: fixed `Re * sqrt(CL)`. `Re_eff = Re / sqrt(|CL|)`.
+    Type2,
+    /// Type 3: fixed `Re * CL`. `Re_eff = Re / |CL|`.
+    Type3,
+}
+
+impl Default for ReType {
+    fn default() -> Self {
+        Self::Type1
+    }
+}
+
 #[derive(Debug, Clone, Copy, PartialEq)]
 pub enum OperatingMode {
     PrescribedAlpha,
@@ -19,6 +35,7 @@ pub struct XfoilOptions {
     pub tolerance: f64,
     pub wake_length_chords: f64,
     pub operating_mode: OperatingMode,
+    pub re_type: ReType,
 }
 
 impl Default for XfoilOptions {
@@ -31,6 +48,7 @@ impl Default for XfoilOptions {
             tolerance: 1.0e-4,
             wake_length_chords: 1.0,
             operating_mode: OperatingMode::PrescribedAlpha,
+            re_type: ReType::Type1,
         }
     }
 }
@@ -46,6 +64,17 @@ impl XfoilOptions {
         self
     }
 
+    pub fn effective_reynolds(&self, cl: f64) -> f64 {
+        if cl.abs() < 1e-6 {
+            return self.reynolds;
+        }
+        match self.re_type {
+            ReType::Type1 => self.reynolds,
+            ReType::Type2 => self.reynolds / cl.abs().sqrt(),
+            ReType::Type3 => self.reynolds / cl.abs(),
+        }
+    }
+
     pub fn validate(&self) -> Result<(), &'static str> {
         if self.reynolds <= 0.0 {
             return Err("Reynolds number must be positive");
@@ -65,3 +94,64 @@ impl XfoilOptions {
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn opts(re_type: ReType, reynolds: f64) -> XfoilOptions {
+        XfoilOptions {
+            reynolds,
+            re_type,
+            ..Default::default()
+        }
+    }
+
+    #[test]
+    fn type1_returns_nominal() {
+        let o = opts(ReType::Type1, 1e6);
+        assert_eq!(o.effective_reynolds(0.5), 1e6);
+        assert_eq!(o.effective_reynolds(-1.2), 1e6);
+    }
+
+    #[test]
+    fn type2_scales_by_sqrt_cl() {
+        let o = opts(ReType::Type2, 1e6);
+        let re = o.effective_reynolds(0.25);
+        // Re / sqrt(0.25) = 1e6 / 0.5 = 2e6
+        assert!((re - 2e6).abs() < 1.0);
+    }
+
+    #[test]
+    fn type3_scales_by_cl() {
+        let o = opts(ReType::Type3, 1e6);
+        let re = o.effective_reynolds(0.5);
+        // Re / 0.5 = 2e6
+        assert!((re - 2e6).abs() < 1.0);
+    }
+
+    #[test]
+    fn negative_cl_uses_abs() {
+        let o = opts(ReType::Type2, 1e6);
+        let re_pos = o.effective_reynolds(0.5);
+        let re_neg = o.effective_reynolds(-0.5);
+        assert!((re_pos - re_neg).abs() < 1.0);
+    }
+
+    #[test]
+    fn near_zero_cl_returns_nominal() {
+        let o = opts(ReType::Type2, 1e6);
+        assert_eq!(o.effective_reynolds(0.0), 1e6);
+        assert_eq!(o.effective_reynolds(1e-8), 1e6);
+
+        let o3 = opts(ReType::Type3, 1e6);
+        assert_eq!(o3.effective_reynolds(0.0), 1e6);
+    }
+
+    #[test]
+    fn default_is_type1() {
+        assert_eq!(ReType::default(), ReType::Type1);
+        let o = XfoilOptions::default();
+        assert_eq!(o.re_type, ReType::Type1);
+    }
+}

crates/rustfoil-xfoil/src/lib.rs

@@ -20,7 +20,7 @@ pub mod state_ops;
 pub mod update;
 pub mod wake_panel;
 
-pub use config::{OperatingMode, XfoilOptions};
+pub use config::{OperatingMode, ReType, XfoilOptions};
 pub use error::{Result, XfoilError};
 pub use oper::{
     build_state_from_coords, coords_from_body, solve_body_oper_point, solve_coords_oper_point,

crates/rustfoil-xfoil/src/oper.rs

@@ -140,6 +140,16 @@ pub fn solve_operating_point_from_state(
             cl_inv, cm_inv
         );
     }
+    // Seed state.cl from inviscid solution so Type 2/3 have an initial CL estimate.
+    {
+        let (cl_inv, _cm_inv) = compute_panel_forces_from_gamma(
+            &state.panel_x,
+            &state.panel_y,
+            &state.qinv,
+            state.alpha_rad,
+        );
+        state.cl = cl_inv;
+    }
     stfind(state);
     iblpan(state);
     xicalc(state);
@@ -163,14 +173,15 @@ pub fn solve_operating_point_from_state(
     }
 
     for iter in 1..=options.max_iterations {
-        let mut assembly = setbl(state, options.reynolds, options.ncrit, options.mach, iter);
+        let re_eff = options.effective_reynolds(state.cl);
+        let mut assembly = setbl(state, re_eff, options.ncrit, options.mach, iter);
         let solve = blsolv(state, &mut assembly, iter);
         update(
             state,
             &assembly,
             &solve,
             options.mach,
-            options.reynolds,
+            re_eff,
             iter,
         );
         if let OperatingMode::PrescribedCl { .. } = state.operating_mode {
@@ -179,7 +190,7 @@ pub fn solve_operating_point_from_state(
         }
         if is_debug_active() {
             // Match XFOIL: dump BL state after UPDATE, before QVFUE/GAMQV/STMOVE.
-            emit_full_state(state, iter, options.mach, options.reynolds);
+            emit_full_state(state, iter, options.mach, re_eff);
         }
         qvfue(state);
         gamqv(state);
@@ -190,7 +201,7 @@ pub fn solve_operating_point_from_state(
         if std::env::var("RUSTFOIL_DISABLE_STMOVE").is_err() {
             stmove(state);
         }
-        update_force_state(state, options.mach, options.reynolds);
+        update_force_state(state, options.mach, re_eff);
         state.iterations = iter;
         state.residual = solve.rms;
         state.converged = solve.rms <= options.tolerance;
@@ -211,6 +222,7 @@ pub fn solve_operating_point_from_state(
         }
     }
 
+    let final_re_eff = options.effective_reynolds(state.cl);
     Ok(XfoilViscousResult {
         alpha_deg: state.alpha_rad.to_degrees(),
         cl: state.cl,
@@ -223,8 +235,9 @@ pub fn solve_operating_point_from_state(
         residual: state.residual,
         cd_friction: state.cdf,
         cd_pressure: state.cdp,
-        x_separation: separation_x(state, crate::state::XfoilSurface::Upper, options.mach, options.reynolds)
-            .or_else(|| separation_x(state, crate::state::XfoilSurface::Lower, options.mach, options.reynolds)),
+        x_separation: separation_x(state, crate::state::XfoilSurface::Upper, options.mach, final_re_eff)
+            .or_else(|| separation_x(state, crate::state::XfoilSurface::Lower, options.mach, final_re_eff)),
+        reynolds_eff: final_re_eff,
     })
 }
 

crates/rustfoil-xfoil/src/result.rs

@@ -14,6 +14,7 @@ pub struct XfoilViscousResult {
     pub cd_friction: f64,
     pub cd_pressure: f64,
     pub x_separation: Option<f64>,
+    pub reynolds_eff: f64,
 }
 
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]

crates/rustfoil-xfoil/tests/support.rs

@@ -549,6 +549,7 @@ pub fn run_workflow_case(alpha_deg: f64) -> rustfoil_xfoil::result::XfoilViscous
             tolerance: 1.0e-4,
             wake_length_chords: 1.0,
             operating_mode: OperatingMode::PrescribedAlpha,
+            ..Default::default()
         },
     )
     .expect("run workflow case")