Replace the C++ implementation with the Rust port

The Rust port (formerly gipfel-rs) is now the official successor to the
FLTK/GSL C++ code base. This squashes the port into a single commit and
removes the C++ sources (they remain in the git history below).

Layout:
  crates/gipfel-core  headless core: hills, panorama geometry, projections
                      (maxima-generated derivatives), pure-Rust LM fit,
                      gipfel: JPEG comment + EXIF I/O, profiles, stitching
  crates/gipfel-cli   command-line tools (export, info, stitch, profiles)
  crates/gipfel-app   Cosmic DE (libcosmic) GUI, binary name `gipfel`
  codegen/            maxima scripts + regen.sh for the projection code
  tests/              golden files + C++ reference harnesses (historical)

The original README is kept for now; it will be updated to describe the
Rust build. gipfel.dat, AUTHORS, ChangeLog and NEWS are unchanged.

16 files changed, 3709 insertions, 0 deletions

diff --git a/crates/gipfel-core/Cargo.toml b/crates/gipfel-core/Cargo.toml
new file mode 100644
index 0000000..d8ddfd9
--- /dev/null
+++ b/crates/gipfel-core/Cargo.toml
@@ -0,0 +1,13 @@
+[package]
+name = "gipfel-core"
+version = "0.1.0"
+edition = "2021"
+description = "Port of gipfel's computational core (hills, panorama, projections, least-squares fit, metadata)"
+license = "GPL-2.0-or-later"
+
+[dependencies]
+kamadak-exif.workspace = true
+quick-xml.workspace = true
+image.workspace = true
+tiff.workspace = true
+dirs.workspace = true
diff --git a/crates/gipfel-core/examples/genjpg.rs b/crates/gipfel-core/examples/genjpg.rs
new file mode 100644
index 0000000..ae12b94
--- /dev/null
+++ b/crates/gipfel-core/examples/genjpg.rs
@@ -0,0 +1,24 @@
+// creates a test JPEG with a gipfel comment: genjpg <out.jpg>
+use gipfel_core::metadata::ImageMetaData;
+
+fn main() {
+    let out = std::env::args().nth(1).unwrap();
+    let (w, h) = (1600u32, 900u32);
+    let img = image::DynamicImage::new_rgb8(w, h).grayscale();
+    img.save(&out).unwrap();
+
+    let mut md = ImageMetaData::new();
+    md.longitude = 11.15;
+    md.latitude = 47.36;
+    md.height = 2100.0;
+    md.direction = 20.0;
+    md.pitch = 0.0;
+    md.roll = 1.0;
+    md.focal_length_35mm = 35.0;
+    md.projection_type = 0;
+    md.k0 = 0.0;
+    md.k1 = 0.0;
+    md.x0 = 0.0;
+    md.save_image(&out, &out).unwrap();
+    println!("wrote {}", out);
+}
diff --git a/crates/gipfel-core/src/hill.rs b/crates/gipfel-core/src/hill.rs
new file mode 100644
index 0000000..a904948
--- /dev/null
+++ b/crates/gipfel-core/src/hill.rs
@@ -0,0 +1,224 @@
+//! Hill database: data structures, CSV loading, duplicate detection.
+//!
+//! Faithful port of gipfel's `Hill.cxx`.
+
+use std::fs::File;
+use std::io::{BufRead, BufReader};
+
+pub const FLAG_DUPLIC: u32 = 1 << 0;
+pub const FLAG_KNOWN: u32 = 1 << 1;
+pub const FLAG_HIDDEN: u32 = 1 << 2;
+pub const FLAG_VISIBLE: u32 = 1 << 3;
+pub const FLAG_TRACK_POINT: u32 = 1 << 4;
+pub const FLAG_EXPORT: u32 = 1 << 5;
+
+pub const PI: f64 = std::f64::consts::PI;
+pub const DEG2RAD: f64 = PI / 180.0;
+
+#[derive(Clone, Debug)]
+pub struct Hill {
+    pub name: String,
+    /// latitude, radians
+    pub phi: f64,
+    /// longitude, radians
+    pub lam: f64,
+    /// height, meters
+    pub height: f64,
+    /// azimuth from viewpoint, radians
+    pub alph: f64,
+    /// peak angle over horizon, radians
+    pub a_elev: f64,
+    /// distance from viewpoint, fraction of earth radius
+    pub dist: f64,
+    pub x: f64,
+    pub y: f64,
+    pub label_x: f64,
+    pub label_y: f64,
+    pub flags: u32,
+}
+
+impl Hill {
+    pub fn new(name: &str, phi: f64, lam: f64, height: f64) -> Self {
+        Hill {
+            name: name.to_string(),
+            phi,
+            lam,
+            height,
+            alph: 0.0,
+            a_elev: 0.0,
+            dist: 0.0,
+            x: 0.0,
+            y: 0.0,
+            label_x: 0.0,
+            label_y: 0.0,
+            flags: 0,
+        }
+    }
+
+    pub fn is_track_point(&self) -> bool {
+        self.flags & FLAG_TRACK_POINT != 0
+    }
+}
+
+/// `atof()` semantics: parse the longest numeric prefix, 0.0 on failure.
+fn atof(s: &str) -> f64 {
+    let s = s.trim();
+    if s.is_empty() {
+        return 0.0;
+    }
+    if let Ok(v) = s.parse::<f64>() {
+        return v;
+    }
+    for i in (1..s.len()).rev() {
+        if let Ok(v) = s[..i].parse::<f64>() {
+            return v;
+        }
+    }
+    0.0
+}
+
+/// Load a hill database. Two line formats are accepted:
+/// - `SHORT,NAME,DESC,LAT,LON,HEIGHT` (6 fields)
+/// - `LAT,LON,HEIGHT` (3 fields, track point)
+pub fn load(path: &str) -> std::io::Result<Vec<Hill>> {
+    let file = File::open(path)?;
+    let mut hills = Vec::new();
+
+    // the database is Latin-1; decode to Unicode code points
+    for line in BufReader::new(file).split(b'\n') {
+        let raw = line?;
+        let line: String = raw.iter().map(|b| *b as char).collect();
+        // strsep() semantics: split into at most 10 fields, keeping empties
+        let vals: Vec<&str> = line.split(',').take(10).collect();
+        let n = vals.len();
+
+        if n == 6 {
+            // standard format including name and description
+            let phi = atof(vals[3]) * DEG2RAD;
+            let lam = atof(vals[4]) * DEG2RAD;
+            let height = atof(vals[5]);
+            hills.push(Hill::new(vals[1], phi, lam, height));
+        } else if n == 3 {
+            // track point format
+            let phi = atof(vals[0]) * DEG2RAD;
+            let lam = atof(vals[1]) * DEG2RAD;
+            let height = atof(vals[2]);
+            let mut h = Hill::new("", phi, lam, height);
+            // NOTE: original does not set TRACK_POINT here
+            h.flags |= 0;
+            hills.push(h);
+        }
+    }
+
+    Ok(hills)
+}
+
+/// Mark hills that are (near-)exact duplicates of an earlier hill.
+/// Faithful port of `Hills::mark_duplicates`.
+///
+/// Unlike the C original (which sorts the pointer array in place), this
+/// must NOT reorder `hills`: in this port, hill indices are stable handles
+/// (track lists, known-hill sets, visible lists). Sorting an index
+/// permutation instead keeps the semantics identical while preserving
+/// insertion order.
+pub fn mark_duplicates(hills: &mut [Hill], dist: f64) {
+    // descending phi (mirrors the C qsort comparator order)
+    let mut order: Vec<usize> = (0..hills.len()).collect();
+    order.sort_by(|&a, &b| {
+        hills[b]
+            .phi
+            .partial_cmp(&hills[a].phi)
+            .unwrap_or(std::cmp::Ordering::Equal)
+    });
+
+    let len = hills.len();
+    for pos in 0..len {
+        let i = order[pos];
+        if hills[i].is_track_point() {
+            continue;
+        }
+
+        let mut j = pos + 1;
+        while j < len && (hills[order[j]].phi - hills[i].phi).abs() <= dist {
+            let jj = order[j];
+            if hills[jj].flags & FLAG_DUPLIC == 0
+                && (hills[jj].lam - hills[i].lam).abs() <= dist
+                && (hills[jj].height - hills[i].height).abs() <= 50.0
+            {
+                hills[jj].flags |= FLAG_DUPLIC;
+            }
+            j += 1;
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn atof_semantics() {
+        assert_eq!(atof("12.5"), 12.5);
+        assert_eq!(atof(""), 0.0);
+        assert_eq!(atof("abc"), 0.0);
+        assert_eq!(atof("12abc"), 12.0);
+        assert_eq!(atof(" -3.5x"), -3.5);
+    }
+
+    #[test]
+    fn load_formats() {
+        let dir = std::env::temp_dir().join("gipfel-core-test");
+        std::fs::create_dir_all(&dir).unwrap();
+        let path = dir.join("hills.csv");
+        std::fs::write(
+            &path,
+            "# comment line\n\
+             ZSP,Zugspitze,Germany,47.421,10.9855,2962\n\
+             FOO,Foo Peak,,47.5,11.0,1000\n\
+             47.6,11.1,1234\n\
+             bad,line\n",
+        )
+        .unwrap();
+
+        let hills = load(path.to_str().unwrap()).unwrap();
+        assert_eq!(hills.len(), 3);
+        assert_eq!(hills[0].name, "Zugspitze");
+        assert!((hills[0].phi - 47.421 * DEG2RAD).abs() < 1e-12);
+        assert_eq!(hills[0].height, 2962.0);
+        assert_eq!(hills[2].name, "");
+        let _ = std::fs::remove_file(&path);
+    }
+
+    #[test]
+    fn duplicates() {
+        let mut hills = vec![
+            Hill::new("A", 47.0, 11.0, 1000.0),
+            Hill::new("A2", 47.000001, 11.000001, 1010.0), // dup of A
+            Hill::new("B", 47.5, 11.5, 1000.0),
+            Hill::new("C", 47.0, 11.0, 1200.0), // >50m diff -> not dup
+        ];
+        mark_duplicates(&mut hills, 0.00001);
+        assert_eq!(hills.iter().filter(|h| h.flags & FLAG_DUPLIC != 0).count(), 1);
+    }
+
+    #[test]
+    fn mark_duplicates_preserves_order() {
+        // indices are stable handles in the port: mark_duplicates must not
+        // reorder the vec (regression: an in-place phi sort broke the app's
+        // track index range when add_hills ran it on the full hill list)
+        let mut hills = vec![
+            Hill::new("south", 46.0, 11.0, 1000.0),
+            Hill::new("north", 48.0, 11.0, 1000.0),
+            Hill::new("mid", 47.0, 11.0, 1000.0),
+            Hill::new("mid-dup", 47.000001, 11.0, 1010.0),
+        ];
+        mark_duplicates(&mut hills, 0.00001);
+        assert_eq!(hills[0].name, "south");
+        assert_eq!(hills[1].name, "north");
+        assert_eq!(hills[2].name, "mid");
+        assert_eq!(hills[3].name, "mid-dup");
+        // descending phi visits mid-dup first, so 'mid' is the one marked
+        assert_eq!(hills[2].flags & FLAG_DUPLIC, FLAG_DUPLIC);
+        assert_eq!(hills[3].flags & FLAG_DUPLIC, 0);
+    }
+}
diff --git a/crates/gipfel-core/src/lib.rs b/crates/gipfel-core/src/lib.rs
new file mode 100644
index 0000000..2e2fac9
--- /dev/null
+++ b/crates/gipfel-core/src/lib.rs
@@ -0,0 +1,23 @@
+//! gipfel-core: headless port of the gipfel computational core.
+//!
+//! Modules:
+//! - [`hill`]: hill database (CSV load, duplicate marking)
+//! - [`panorama`]: viewpoint, geodesy, refraction, occlusion, visibility
+//! - [`projection`]: rectilinear/cylindrical projections with symbolic
+//!   derivatives generated by maxima (see `codegen/`)
+//! - [`lsq`]: Levenberg-Marquardt fit of the view parameters (system GSL)
+//! - [`metadata`]: `gipfel:` JPEG comment + EXIF read/write
+
+pub mod hill;
+pub mod lsq;
+pub mod metadata;
+pub mod panorama;
+pub mod profile;
+pub mod projection;
+pub mod stitch;
+pub mod track;
+
+pub use hill::Hill;
+pub use metadata::ImageMetaData;
+pub use panorama::Panorama;
+pub use projection::{Projection, ProjectionType, ViewParams};
diff --git a/crates/gipfel-core/src/lsq.rs b/crates/gipfel-core/src/lsq.rs
new file mode 100644
index 0000000..6a9d6c7
--- /dev/null
+++ b/crates/gipfel-core/src/lsq.rs
@@ -0,0 +1,392 @@
+//! Non-linear least squares fit of the view parameters.
+//!
+//! Faithful port of `ProjectionLSQ.cxx`. The C++ original uses GSL's
+//! `gsl_multifit_fdfsolver_lmsder`; this is a self-contained Levenberg-
+//! Marquardt with the same scaled damping (`JᵀJ + λ·diag(JᵀJ)`) and the
+//! analytic Jacobian produced by the maxima-generated derivatives. For the
+//! well-conditioned gipfel fit (residuals reach zero at the optimum) it
+//! converges to the same minimum; the golden tests verify agreement with the
+//! C/GSL output.
+
+use crate::hill::Hill;
+use crate::projection::{Projection, ViewParams};
+
+const MAX_ITER: usize = 200;
+
+/// Number of fitted parameters at a given level.
+fn np_for_level(level: usize) -> usize {
+    match level {
+        0 => 2,
+        1 => 4,
+        _ => 7,
+    }
+}
+
+struct FitData<'a> {
+    proj: &'a dyn Projection,
+    hills: &'a [Hill],
+    known: &'a [usize],
+    level: usize,
+    /// parameters at solver start; entries above `level` keep these values
+    parms0: ViewParams,
+}
+
+impl FitData<'_> {
+    fn params_from_x(&self, x: &[f64]) -> ViewParams {
+        // C starts from the initial params and only overwrites the entries
+        // that are actually fitted at this level.
+        let mut p = self.parms0;
+        p.a_center = x[0];
+        p.a_pitch = x[1];
+        if self.level > 0 {
+            p.a_roll = x[2];
+            p.scale = x[3];
+        }
+        if self.level > 1 {
+            p.k0 = x[4];
+            p.k1 = x[5];
+            p.x0 = x[6];
+        }
+        p
+    }
+
+    fn args(&self, p: &ViewParams, i: usize) -> crate::projection::ProjArgs {
+        let m = &self.hills[self.known[i]];
+        crate::projection::ProjArgs::from_params(p, m.alph, m.a_elev)
+    }
+
+    /// Residual vector: 2 entries (x,y) per known hill, length `2*known.len()`.
+    fn residuals(&self, x: &[f64], f: &mut [f64]) {
+        let p = self.params_from_x(x);
+        for i in 0..self.known.len() {
+            let m = &self.hills[self.known[i]];
+            let a = self.args(&p, i);
+            f[2 * i] = self.proj.mac_x(a) - m.x;
+            f[2 * i + 1] = self.proj.mac_y(a) - m.y;
+        }
+    }
+
+    /// Jacobian, row-major `m`×`np` (`m = 2*known.len()`).
+    fn jacobian(&self, x: &[f64], j: &mut [f64]) {
+        let p = self.params_from_x(x);
+        let np = np_for_level(self.level);
+        let row = |r: usize, c: usize, v: f64, j: &mut [f64]| j[r * np + c] = v;
+        for i in 0..self.known.len() {
+            let a = self.args(&p, i);
+
+            row(2 * i, 0, self.proj.mac_x_dc_view(a), j);
+            row(2 * i + 1, 0, self.proj.mac_y_dc_view(a), j);
+            row(2 * i, 1, self.proj.mac_x_dc_pitch(a), j);
+            row(2 * i + 1, 1, self.proj.mac_y_dc_pitch(a), j);
+
+            if np > 2 {
+                row(2 * i, 2, self.proj.mac_x_dc_roll(a), j);
+                row(2 * i + 1, 2, self.proj.mac_y_dc_roll(a), j);
+                row(2 * i, 3, self.proj.mac_x_dscale(a), j);
+                row(2 * i + 1, 3, self.proj.mac_y_dscale(a), j);
+            }
+
+            if np > 4 {
+                row(2 * i, 4, self.proj.mac_x_dk0(a), j);
+                row(2 * i + 1, 4, self.proj.mac_y_dk0(a), j);
+                row(2 * i, 5, self.proj.mac_x_dk1(a), j);
+                row(2 * i + 1, 5, self.proj.mac_y_dk1(a), j);
+                row(2 * i, 6, self.proj.mac_x_dx0(a), j);
+                row(2 * i + 1, 6, self.proj.mac_y_dx0(a), j);
+            }
+        }
+    }
+}
+
+fn dot(a: &[f64], b: &[f64]) -> f64 {
+    a.iter().zip(b).map(|(x, y)| x * y).sum()
+}
+
+/// Solve `A·δ = b` in place (A row-major `n`×`n`, overwritten) by Gaussian
+/// elimination with partial pivoting. Returns `None` if the matrix is
+/// (numerically) singular.
+fn gauss_solve(a: &mut [f64], b: &[f64], n: usize) -> Option<Vec<f64>> {
+    let mut x = b.to_vec();
+    for col in 0..n {
+        // partial pivot
+        let mut piv = col;
+        let mut maxv = a[col * n + col].abs();
+        for r in (col + 1)..n {
+            let v = a[r * n + col].abs();
+            if v > maxv {
+                maxv = v;
+                piv = r;
+            }
+        }
+        if maxv < 1e-300 {
+            return None;
+        }
+        if piv != col {
+            for c in 0..n {
+                a.swap(col * n + c, piv * n + c);
+            }
+            x.swap(col, piv);
+        }
+        let diag = a[col * n + col];
+        for r in (col + 1)..n {
+            let factor = a[r * n + col] / diag;
+            a[r * n + col] = 0.0;
+            for c in (col + 1)..n {
+                a[r * n + c] -= factor * a[col * n + c];
+            }
+            x[r] -= factor * x[col];
+        }
+    }
+    for r in (0..n).rev() {
+        let mut s = x[r];
+        for c in (r + 1)..n {
+            s -= a[r * n + c] * x[c];
+        }
+        x[r] = s / a[r * n + r];
+    }
+    Some(x)
+}
+
+/// One Levenberg-Marquardt fit at a single parameter level. On return, `x`
+/// holds the fitted parameters (length `np`).
+fn lm(data: &FitData, x: &mut [f64], np: usize) {
+    let m = 2 * data.known.len();
+    let mut f = vec![0.0f64; m];
+    let mut fnew = vec![0.0f64; m];
+    let mut jac = vec![0.0f64; m * np];
+    let mut a = vec![0.0f64; np * np];
+    let mut ad = vec![0.0f64; np * np];
+    let mut g = vec![0.0f64; np];
+    let mut neg_g = vec![0.0f64; np];
+    let mut xnew = vec![0.0f64; np];
+
+    data.residuals(x, &mut f);
+    let mut cost = 0.5 * dot(&f, &f);
+    let mut lambda: f64 = 1e-3;
+
+    for _ in 0..MAX_ITER {
+        if !cost.is_finite() {
+            break;
+        }
+        data.jacobian(x, &mut jac);
+
+        // A = JᵀJ (symmetric), g = Jᵀf
+        for i in 0..np {
+            for k in i..np {
+                let mut s = 0.0;
+                for r in 0..m {
+                    s += jac[r * np + i] * jac[r * np + k];
+                }
+                a[i * np + k] = s;
+                a[k * np + i] = s;
+            }
+            let mut s = 0.0;
+            for r in 0..m {
+                s += jac[r * np + i] * f[r];
+            }
+            g[i] = s;
+        }
+
+        // Solve the damped normal equations (JᵀJ + λ·diag(JᵀJ))·δ = -Jᵀf,
+        // inflating λ until the system is non-singular.
+        let mut delta: Option<Vec<f64>> = None;
+        for _ in 0..50 {
+            ad.copy_from_slice(&a);
+            for i in 0..np {
+                let d = a[i * np + i];
+                ad[i * np + i] = if d > 1e-300 { d * (1.0 + lambda) } else { lambda.max(1e-300) };
+            }
+            for i in 0..np {
+                neg_g[i] = -g[i];
+            }
+            let mut ad_try = ad.clone();
+            if let Some(d) = gauss_solve(&mut ad_try, &neg_g, np) {
+                delta = Some(d);
+                break;
+            }
+            lambda *= 10.0;
+        }
+        let delta = match delta {
+            Some(d) => d,
+            None => break,
+        };
+
+        for i in 0..np {
+            xnew[i] = x[i] + delta[i];
+        }
+        data.residuals(&xnew, &mut fnew);
+        let cost_new = 0.5 * dot(&fnew, &fnew);
+
+        if cost_new < cost {
+            let improvement = cost - cost_new;
+            x.copy_from_slice(&xnew);
+            std::mem::swap(&mut f, &mut fnew);
+            cost = cost_new;
+            lambda = (lambda * 0.1).max(1e-12);
+            // Converged when the cost can no longer improve meaningfully.
+            if improvement <= 1e-18 * cost {
+                break;
+            }
+        } else {
+            lambda *= 10.0;
+            if lambda > 1e16 {
+                break;
+            }
+        }
+    }
+}
+
+/// Run the Levenberg-Marquardt fit for one parameter level and write the
+/// result back into `parms`. Mirrors `ProjectionLSQ::lsq`.
+fn lsq(proj: &dyn Projection, hills: &[Hill], known: &[usize], parms: &mut ViewParams, level: usize) {
+    let np = np_for_level(level);
+    let data = FitData {
+        proj,
+        hills,
+        known,
+        level,
+        parms0: *parms,
+    };
+
+    let mut x = vec![
+        parms.a_center,
+        parms.a_pitch,
+        parms.a_roll,
+        parms.scale,
+        parms.k0,
+        parms.k1,
+        parms.x0,
+    ];
+    x.truncate(np);
+
+    lm(&data, &mut x, np);
+
+    parms.a_center = x[0];
+    parms.a_pitch = x[1];
+    if level > 0 {
+        parms.a_roll = x[2];
+        parms.scale = x[3];
+    }
+    if level > 1 {
+        parms.k0 = x[4];
+        parms.k1 = x[5];
+        parms.x0 = x[6];
+    }
+}
+
+/// Fit the view parameters to the known hills.
+/// Mirrors `ProjectionLSQ::comp_params`. `known` indexes into `hills`.
+/// With `fit_distortion == false` the lens distortion parameters
+/// (k0/k1/x0) are held at their current values and only the viewpoint
+/// parameters are fitted, even for 4+ known hills.
+pub fn comp_params(
+    proj: &dyn Projection,
+    hills: &[Hill],
+    known: &[usize],
+    parms: &mut ViewParams,
+    fit_distortion: bool,
+) -> Result<(), String> {
+    let n = known.len();
+    if n < 1 {
+        return Err("Please position at least 1 hill".to_string());
+    }
+
+    let m1 = &hills[known[0]];
+    let mut new_parms = *parms;
+    new_parms.a_center = m1.alph;
+    new_parms.a_pitch = 0.0;
+
+    if n > 1 {
+        let m2 = &hills[known[1]];
+        let scale_tmp = (m2.x - m1.x).abs() / (m2.alph - m1.alph).abs();
+        if scale_tmp.is_nan() || scale_tmp < 50.0 {
+            return Err(format!("Could not determine initial scale value ({scale_tmp})"));
+        }
+        new_parms.scale = scale_tmp;
+        new_parms.a_roll = 0.0;
+    }
+
+    if n > 3 && fit_distortion {
+        new_parms.k0 = 0.0;
+        new_parms.k1 = 0.0;
+        new_parms.x0 = 0.0;
+        eprintln!("Performing calibration");
+    }
+
+    let level = if n == 1 { 0 } else { 1 };
+    lsq(proj, hills, known, &mut new_parms, level);
+    if n >= 4 && fit_distortion {
+        lsq(proj, hills, known, &mut new_parms, 2);
+    }
+
+    if new_parms.scale.is_nan() || new_parms.scale < 50.0 {
+        return Err("Could not determine reasonable view parameters".to_string());
+    }
+
+    *parms = new_parms;
+    Ok(())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::hill::Hill;
+    use crate::projection::{Projection, ProjectionType, ViewParams, projection_from_type};
+
+    /// Five observed hills projected with known parameters (incl. distortion).
+    fn fixture() -> (Box<dyn Projection>, Vec<Hill>, Vec<usize>, ViewParams) {
+        let proj = projection_from_type(ProjectionType::Rectilinear);
+        let truth = ViewParams {
+            a_center: 0.3,
+            scale: 2000.0,
+            a_pitch: 0.05,
+            a_roll: 0.0,
+            k0: -0.0012,
+            k1: 2e-5,
+            x0: 0.0,
+        };
+        let mut hills = Vec::new();
+        let mut known = Vec::new();
+        for i in 0..5 {
+            let alph = 0.15 + i as f64 * 0.07;
+            let elev = 0.02 + i as f64 * 0.015;
+            let (x, y) = proj.get_coordinates(alph, elev, &truth);
+            let mut h = Hill::new(&format!("h{i}"), 0.0, 0.0, 1000.0);
+            h.alph = alph;
+            h.a_elev = elev;
+            h.x = x;
+            h.y = y;
+            hills.push(h);
+            known.push(i);
+        }
+        (proj, hills, known, truth)
+    }
+
+    #[test]
+    fn locked_distortion_stays_untouched() {
+        let (proj, hills, known, truth) = fixture();
+        let mut p = truth;
+        p.a_center -= 0.01;
+        p.scale *= 1.02;
+        comp_params(&*proj, &hills, &known, &mut p, false).unwrap();
+        // 5 known hills, but the distortion params must not move
+        assert_eq!(p.k0, truth.k0);
+        assert_eq!(p.k1, truth.k1);
+        assert_eq!(p.x0, truth.x0);
+        // the viewpoint params did fit
+        assert!((p.a_center - truth.a_center).abs() < 1e-6);
+    }
+
+    #[test]
+    fn unlocked_distortion_fit_recovers() {
+        let (proj, hills, known, truth) = fixture();
+        let mut p = truth;
+        p.k0 = 0.0;
+        p.k1 = 0.0;
+        p.a_center -= 0.01;
+        p.scale *= 1.02;
+        comp_params(&*proj, &hills, &known, &mut p, true).unwrap();
+        assert!((p.k0 - truth.k0).abs() < 1e-5, "k0={}", p.k0);
+        assert!((p.k1 - truth.k1).abs() < 1e-7, "k1={}", p.k1);
+    }
+}
diff --git a/crates/gipfel-core/src/metadata.rs b/crates/gipfel-core/src/metadata.rs
new file mode 100644
index 0000000..e2a888f
--- /dev/null
+++ b/crates/gipfel-core/src/metadata.rs
@@ -0,0 +1,697 @@
+//! Image metadata: the `gipfel:` JPEG comment (read/write, byte-compatible
+//! with the C version) plus EXIF fallback (GPS position, focal length).
+//!
+//! Faithful port of gipfel's `ImageMetaData.cxx`.
+
+use std::fs::File;
+use std::io::{BufReader, Write};
+
+use exif::{In, Tag, Value};
+
+pub const NAN: f64 = f64::NAN;
+
+#[derive(Clone, Debug)]
+pub struct ImageMetaData {
+    pub manufacturer: Option<String>,
+    pub model: Option<String>,
+    pub longitude: f64,
+    pub latitude: f64,
+    pub height: f64,
+    pub direction: f64,
+    /// camera pitch, degrees (written as the "nick" key in the gipfel:
+    /// comment for backwards compatibility with the C original)
+    pub pitch: f64,
+    /// camera roll, degrees (written as the "tilt" key for compatibility)
+    pub roll: f64,
+    pub k0: f64,
+    pub k1: f64,
+    pub x0: f64,
+    pub focal_length: f64,
+    pub focal_length_35mm: f64,
+    pub projection_type: i32,
+}
+
+impl Default for ImageMetaData {
+    fn default() -> Self {
+        ImageMetaData {
+            manufacturer: None,
+            model: None,
+            longitude: NAN,
+            latitude: NAN,
+            height: NAN,
+            direction: NAN,
+            pitch: NAN,
+            roll: NAN,
+            k0: NAN,
+            k1: NAN,
+            x0: NAN,
+            focal_length: NAN,
+            focal_length_35mm: NAN,
+            projection_type: 0,
+        }
+    }
+}
+
+impl ImageMetaData {
+    pub fn new() -> Self {
+        Self::default()
+    }
+
+    /// Read JPEG comment first, then fill remaining gaps from EXIF.
+    pub fn load_image(&mut self, path: &str) {
+        *self = ImageMetaData::new();
+        if let Some(com) = read_jpeg_comment(path) {
+            self.parse_gipfel_comment(&com);
+        }
+        self.load_image_exif(path);
+    }
+
+    /// Write the gipfel comment into a copy of `in_img`, saved as `out_img`.
+    /// Writes atomically via a temp file, like the C original.
+    pub fn save_image(&self, in_img: &str, out_img: &str) -> std::io::Result<()> {
+        let comment = self.format_comment();
+        let tmp = format!("{}.tmp", out_img);
+        write_jpeg_comment(in_img, &tmp, comment.as_bytes())?;
+        std::fs::rename(&tmp, out_img)?;
+        Ok(())
+    }
+
+    pub fn format_comment(&self) -> String {
+        format!(
+            "gipfel: longitude {}, latitude {}, height {}, direction {}, nick {}, tilt {}, focal_length_35mm {}, projection type {}, k0 {}, k1 {}, x0 {}",
+            fmt_f64(self.longitude),
+            fmt_f64(self.latitude),
+            fmt_f64(self.height),
+            fmt_f64(self.direction),
+            fmt_f64(self.pitch),
+            fmt_f64(self.roll),
+            fmt_f64(self.focal_length_35mm),
+            self.projection_type,
+            fmt_f64(self.k0),
+            fmt_f64(self.k1),
+            fmt_f64(self.x0),
+        )
+    }
+
+    /// sscanf(gipfel_format_scan, ...) equivalent; returns true when at
+    /// least 8 fields were matched (the C requirement).
+    pub fn parse_gipfel_comment(&mut self, com: &[u8]) -> bool {
+        let s = String::from_utf8_lossy(com);
+        let mut sc = Scanner { s: s.as_ref() };
+
+        if !sc.lit("gipfel: longitude ") {
+            return false;
+        }
+        let Some(lo) = sc.f64() else { return false };
+        if !sc.lit(", latitude ") { return false; }
+        let la = sc.f64();
+        if la.is_none() { return false; }
+        if !sc.lit(", height ") { return false; }
+        let he = sc.f64();
+        if he.is_none() { return false; }
+        if !sc.lit(", direction ") { return false; }
+        let dir = sc.f64();
+        if dir.is_none() { return false; }
+        if !sc.lit(", nick ") { return false; }
+        let pit = sc.f64();
+        if pit.is_none() { return false; }
+        if !sc.lit(", tilt ") { return false; }
+        let rol = sc.f64();
+        if rol.is_none() { return false; }
+        if !sc.lit(", focal_length_35mm ") { return false; }
+        let fr = sc.f64();
+        if fr.is_none() { return false; }
+        if !sc.lit(", projection type ") { return false; }
+        let pt = sc.i32();
+
+        // 8 matches required: through projection type
+        let Some(pt) = pt else { return false };
+
+        self.longitude = lo;
+        self.latitude = la.unwrap();
+        self.height = he.unwrap();
+        self.direction = dir.unwrap();
+        self.pitch = pit.unwrap();
+        self.roll = rol.unwrap();
+        self.focal_length_35mm = fr.unwrap();
+        self.projection_type = pt;
+
+        // n >= 10: k0 and k1 were matched; x0 defaults to 0.0 when only 10
+        if let Some(k0) = sc.f64_after(", k0 ") {
+            if let Some(k1) = sc.f64_after(", k1 ") {
+                self.k0 = k0;
+                self.k1 = k1;
+                self.x0 = sc.f64_after(", x0 ").unwrap_or(0.0);
+            }
+        }
+
+        true
+    }
+
+    fn load_image_exif(&mut self, path: &str) {
+        let Ok(file) = File::open(path) else { return };
+        let mut reader = BufReader::new(file);
+        let Ok(exif) = exif::Reader::new().read_from_container(&mut reader) else {
+            return;
+        };
+
+        if self.manufacturer.is_none() {
+            if let Some(f) = exif.get_field(Tag::Make, In::PRIMARY) {
+                let s = field_string(f);
+                if !s.is_empty() {
+                    self.manufacturer = Some(s);
+                }
+            }
+        }
+        if self.model.is_none() {
+            if let Some(f) = exif.get_field(Tag::Model, In::PRIMARY) {
+                let s = field_string(f);
+                if !s.is_empty() {
+                    self.model = Some(s);
+                }
+            }
+        }
+
+        if self.focal_length.is_nan() {
+            exif_set_value(self, &exif, Tag::FocalLength, In::PRIMARY, |md, v| md.focal_length = v);
+        }
+        if self.focal_length_35mm.is_nan() {
+            exif_set_value(
+                self,
+                &exif,
+                Tag::FocalLengthIn35mmFilm,
+                In::PRIMARY,
+                |md, v| md.focal_length_35mm = v,
+            );
+        }
+        if self.longitude.is_nan() {
+            exif_set_coordinate(self, &exif, Tag::GPSLongitude, In::PRIMARY, |md, v| md.longitude = v);
+            if !self.longitude.is_nan() && exif_ref_is(&exif, Tag::GPSLongitudeRef, "W") {
+                self.longitude = -self.longitude;
+            }
+        }
+        if self.latitude.is_nan() {
+            exif_set_coordinate(self, &exif, Tag::GPSLatitude, In::PRIMARY, |md, v| md.latitude = v);
+            if !self.latitude.is_nan() && exif_ref_is(&exif, Tag::GPSLatitudeRef, "S") {
+                self.latitude = -self.latitude;
+            }
+        }
+        if self.height.is_nan() {
+            exif_set_value(self, &exif, Tag::GPSAltitude, In::PRIMARY, |md, v| md.height = v);
+            // GPSAltitudeRef: 0 = above sea level, 1 = below
+            if !self.height.is_nan() {
+                if let Some(f) = exif.get_field(Tag::GPSAltitudeRef, In::PRIMARY) {
+                    if let Value::Byte(b) = &f.value {
+                        if b.first() == Some(&1) {
+                            self.height = -self.height;
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+/// True if an EXIF ASCII reference tag (GPSLatitudeRef etc.) equals `val`.
+fn exif_ref_is(exif: &exif::Exif, tag: Tag, val: &str) -> bool {
+    exif.get_field(tag, In::PRIMARY)
+        .and_then(|f| match &f.value {
+            Value::Ascii(a) => a.first().map(|v| String::from_utf8_lossy(v).trim() == val),
+            _ => None,
+        })
+        .unwrap_or(false)
+}
+
+fn fmt_f64(v: f64) -> String {
+    if v.is_nan() {
+        // C printf("%f", nan) prints "nan"; keep it parseable by both
+        "nan".to_string()
+    } else {
+        format!("{:.6}", v)
+    }
+}
+
+fn field_string(f: &exif::Field) -> String {
+    let join = |v: &[String]| v.join(" ");
+    match &f.value {
+        Value::Ascii(a) => a
+            .first()
+            .map(|v| String::from_utf8_lossy(v).trim().to_string())
+            .unwrap_or_default(),
+        Value::Rational(r) => join(r.iter().map(|v| v.to_string()).collect::<Vec<_>>().as_slice()),
+        Value::Short(v) => join(v.iter().map(|v| v.to_string()).collect::<Vec<_>>().as_slice()),
+        Value::Long(v) => join(v.iter().map(|v| v.to_string()).collect::<Vec<_>>().as_slice()),
+        Value::SShort(v) => join(v.iter().map(|v| v.to_string()).collect::<Vec<_>>().as_slice()),
+        Value::SLong(v) => join(v.iter().map(|v| v.to_string()).collect::<Vec<_>>().as_slice()),
+        _ => String::new(),
+    }
+}
+
+/// exifSetValue: assign when the tag exists and its value is >= 0.
+fn exif_set_value<F>(md: &mut ImageMetaData, exif: &exif::Exif, tag: Tag, inp: In, set: F)
+where
+    F: FnOnce(&mut ImageMetaData, f64),
+{
+    if let Some(f) = exif.get_field(tag, inp) {
+        if let Some(v) = field_tofloat(f, 0) {
+            if v >= 0.0 {
+                set(md, v);
+            }
+        }
+    }
+}
+
+/// exifSetCoordinate: degrees + minutes/60 + seconds/3600, each component
+/// only used when >= 0 (faithful to the C original, incl. ignoring N/S refs).
+fn exif_set_coordinate<F>(md: &mut ImageMetaData, exif: &exif::Exif, tag: Tag, inp: In, set: F)
+where
+    F: FnOnce(&mut ImageMetaData, f64),
+{
+    let Some(f) = exif.get_field(tag, inp) else {
+        return;
+    };
+
+    let mut dest: Option<f64> = None;
+    for i in 0..3 {
+        if let Some(v) = field_tofloat(f, i) {
+            if v >= 0.0 {
+                dest = Some(match i {
+                    0 => v,
+                    1 => dest.unwrap_or(0.0) + v / 60.0,
+                    _ => dest.unwrap_or(0.0) + v / 3600.0,
+                });
+            }
+        }
+    }
+    if let Some(v) = dest {
+        set(md, v);
+    }
+}
+
+fn field_tofloat(f: &exif::Field, idx: usize) -> Option<f64> {
+    match &f.value {
+        Value::Rational(r) => r.get(idx).map(|v| v.num as f64 / v.denom as f64),
+        Value::Short(v) => v.get(idx).map(|v| *v as f64),
+        Value::Long(v) => v.get(idx).map(|v| *v as f64),
+        Value::SShort(v) => v.get(idx).map(|v| *v as f64),
+        Value::SLong(v) => v.get(idx).map(|v| *v as f64),
+        _ => None,
+    }
+}
+
+// ---------------------------------------------------------------------------
+// strtod-like scanner for the sscanf format
+// ---------------------------------------------------------------------------
+
+struct Scanner<'a> {
+    s: &'a str,
+}
+
+impl Scanner<'_> {
+    fn skip_ws(&mut self) {
+        self.s = self.s.trim_start();
+    }
+
+    fn lit(&mut self, lit: &str) -> bool {
+        if self.s.starts_with(lit) {
+            self.s = &self.s[lit.len()..];
+            true
+        } else {
+            false
+        }
+    }
+
+    /// %lf: skip whitespace, parse the longest numeric prefix.
+    fn f64(&mut self) -> Option<f64> {
+        self.skip_ws();
+        for i in (1..=self.s.len()).rev() {
+            if let Ok(v) = self.s[..i].parse::<f64>() {
+                self.s = &self.s[i..];
+                return Some(v);
+            }
+        }
+        None
+    }
+
+    fn f64_after(&mut self, lit: &str) -> Option<f64> {
+        if self.lit(lit) {
+            self.f64()
+        } else {
+            None
+        }
+    }
+
+    /// %d
+    fn i32(&mut self) -> Option<i32> {
+        self.skip_ws();
+        for i in (1..=self.s.len()).rev() {
+            if let Ok(v) = self.s[..i].parse::<i32>() {
+                self.s = &self.s[i..];
+                return Some(v);
+            }
+        }
+        None
+    }
+}
+
+// ---------------------------------------------------------------------------
+// JPEG COM segment reading / writing
+// ---------------------------------------------------------------------------
+
+fn marker_len_ok(marker: u8) -> bool {
+    // standalone markers without length field
+    marker == 0x01 || (0xD0..=0xD9).contains(&marker)
+}
+
+/// Return the payload of the first COM segment before SOS, if any.
+pub fn read_jpeg_comment(path: &str) -> Option<Vec<u8>> {
+    let bytes = std::fs::read(path).ok()?;
+    if bytes.len() < 2 || bytes[0] != 0xFF || bytes[1] != 0xD8 {
+        return None;
+    }
+
+    let mut i = 2;
+    while i + 1 < bytes.len() {
+        if bytes[i] != 0xFF {
+            return None;
+        }
+        let marker = bytes[i + 1];
+        if marker == 0xFF {
+            i += 1;
+            continue;
+        }
+        if marker == 0xDA || marker == 0xD9 {
+            return None;
+        }
+        if marker_len_ok(marker) {
+            i += 2;
+            continue;
+        }
+        if i + 4 > bytes.len() {
+            return None;
+        }
+        let len = (u16::from_be_bytes([bytes[i + 2], bytes[i + 3]]) as usize).checked_add(2)?;
+        if i + len > bytes.len() {
+            return None;
+        }
+        if marker == 0xFE {
+            return Some(bytes[i + 4..i + 2 + len - 2].to_vec());
+        }
+        i += len;
+    }
+    None
+}
+
+/// Copy `in_img` to `out_img`, dropping all COM segments and inserting a
+/// fresh COM segment right after SOI.
+pub fn write_jpeg_comment(in_img: &str, out_img: &str, comment: &[u8]) -> std::io::Result<()> {
+    let bytes = std::fs::read(in_img)?;
+    if bytes.len() < 2 || bytes[0] != 0xFF || bytes[1] != 0xD8 {
+        return Err(std::io::Error::new(std::io::ErrorKind::InvalidInput, "not a JPEG file"));
+    }
+    if comment.len() + 2 > 0xFFFF {
+        return Err(std::io::Error::new(std::io::ErrorKind::InvalidInput, "comment too long"));
+    }
+
+    let mut out: Vec<u8> = Vec::with_capacity(bytes.len() + comment.len() + 4);
+    out.extend_from_slice(&bytes[..2]); // SOI
+
+    // Copy the SOI-adjacent header region (APPn/COM/DPI/etc.), dropping COM
+    // segments, and insert our fresh COM after the last APPn segment:
+    // that is where exiv2 (and thus the C gipfel) expects it.
+    let mut i = 2;
+    let mut insert_at = out.len();
+    let mut copy_rest_from = bytes.len();
+    while i + 1 < bytes.len() {
+        if bytes[i] != 0xFF {
+            break;
+        }
+        let marker = bytes[i + 1];
+        if marker == 0xFF {
+            i += 1;
+            continue;
+        }
+        if marker == 0xDA || marker == 0xD9 {
+            copy_rest_from = i;
+            break;
+        }
+        if marker_len_ok(marker) {
+            out.push(0xFF);
+            out.push(marker);
+            i += 2;
+            continue;
+        }
+        if i + 4 > bytes.len() {
+            break;
+        }
+        let len = match u16::from_be_bytes([bytes[i + 2], bytes[i + 3]]) {
+            l => (l as usize).checked_add(2).unwrap_or(usize::MAX),
+        };
+        if len < 2 || i + len > bytes.len() {
+            break;
+        }
+        if marker != 0xFE {
+            out.extend_from_slice(&bytes[i..i + len]);
+            if (0xE0..=0xEF).contains(&marker) {
+                insert_at = out.len();
+            }
+        }
+        i += len;
+    }
+    copy_rest_from = copy_rest_from.min(i.max(2));
+
+    let mut com: Vec<u8> = Vec::with_capacity(comment.len() + 4);
+    com.push(0xFF);
+    com.push(0xFE);
+    com.extend_from_slice(&((comment.len() + 2) as u16).to_be_bytes());
+    com.extend_from_slice(comment);
+    out.splice(insert_at..insert_at, com);
+
+    out.extend_from_slice(&bytes[copy_rest_from..]);
+
+    let mut f = File::create(out_img)?;
+    f.write_all(&out)?;
+    f.flush()
+}
+
+/// Read raw image dimensions (width, height) of a JPEG.
+pub fn image_dimensions(path: &str) -> Option<(u32, u32)> {
+    image::image_dimensions(path).ok()
+}
+
+/// Apply image metadata to a panorama like `GipfelWidget::load_image` does.
+/// Returns false when the gipfel info is incomplete (C: `have_gipfel_info`).
+pub fn apply_image_metadata(
+    pan: &mut crate::panorama::Panorama,
+    md: &ImageMetaData,
+    img_w: u32,
+    img_h: u32,
+) -> bool {
+    use crate::projection::ProjectionType;
+
+    let mut have_info = true;
+
+    pan.set_view_long(md.longitude);
+    pan.set_view_lat(md.latitude);
+    pan.set_view_height(md.height);
+
+    pan.set_projection(match md.projection_type {
+        1 => ProjectionType::Cylindrical,
+        _ => ProjectionType::Rectilinear,
+    });
+
+    if md.direction.is_nan() {
+        pan.set_center_angle(0.0);
+        have_info = false;
+    } else {
+        pan.set_center_angle(md.direction);
+    }
+
+    if md.pitch.is_nan() {
+        pan.set_pitch_angle(0.0);
+        have_info = false;
+    } else {
+        pan.set_pitch_angle(md.pitch);
+    }
+
+    if md.roll.is_nan() {
+        pan.set_roll_angle(0.0);
+        have_info = false;
+    } else {
+        pan.set_roll_angle(md.roll);
+    }
+
+    // set_focal_length_35mm: scale = fl * max(w, h) / 35
+    let w = img_w.max(img_h);
+    if md.focal_length_35mm.is_nan() || md.focal_length_35mm == 0.0 {
+        pan.set_scale(35.0 * w as f64 / 35.0);
+        have_info = false;
+    } else {
+        pan.set_scale(md.focal_length_35mm * w as f64 / 35.0);
+    }
+
+    // distortion parameters, C ordering:
+    // 1. gipfel data in JPEG comment
+    // 2. matching distortion profile (camera + focal length)
+    // 3. 0.0, 0.0, 0.0
+    let (mut k0, mut k1, mut x0) = (md.k0, md.k1, md.x0);
+    if k0.is_nan() {
+        if let Some(p) = md
+            .manufacturer
+            .as_deref()
+            .zip(md.model.as_deref())
+            .and_then(|(mfr, model)| crate::profile::profile_name(mfr, model, md.focal_length))
+            .and_then(|name| crate::profile::load(&name).ok())
+        {
+            k0 = p.k0;
+            k1 = p.k1;
+            x0 = p.x0;
+        }
+    }
+    if k0.is_nan() {
+        k0 = 0.0;
+    }
+    if k1.is_nan() {
+        k1 = 0.0;
+    }
+    if x0.is_nan() {
+        x0 = 0.0;
+    }
+    pan.set_distortion_params(k0, k1, x0);
+
+    have_info
+}
+
+/// The profile name auto-matched from this image's camera and focal length
+/// (C `get_distortion_profile_name`). Suggested whenever the camera is
+/// known — [`apply_image_metadata`] only *applies* it when the JPEG
+/// comment carries no distortion params.
+pub fn distortion_profile_name(md: &ImageMetaData) -> Option<String> {
+    md.manufacturer
+        .as_deref()
+        .zip(md.model.as_deref())
+        .and_then(|(mfr, model)| crate::profile::profile_name(mfr, model, md.focal_length))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn jpeg_comment_write_read() {
+        let dir = std::env::temp_dir().join("gipfel-core-mdtest");
+        std::fs::create_dir_all(&dir).unwrap();
+        let img = dir.join("t.jpg");
+        let out = dir.join("o.jpg");
+
+        image::DynamicImage::new_rgb8(64, 48)
+            .grayscale()
+            .save(&img)
+            .unwrap();
+
+        let comment = b"gipfel: longitude 1.000000, latitude 2.000000, height 3.000000, direction 4.000000, nick 5.000000, tilt 6.000000, focal_length_35mm 7.000000, projection type 0, k0 0.000000, k1 0.000000, x0 0.000000";
+        write_jpeg_comment(img.to_str().unwrap(), out.to_str().unwrap(), comment).unwrap();
+
+        // no COM before, one after
+        assert!(read_jpeg_comment(img.to_str().unwrap()).is_none());
+        assert_eq!(read_jpeg_comment(out.to_str().unwrap()).unwrap(), comment);
+
+        // rewriting replaces the existing comment
+        let out2 = dir.join("o2.jpg");
+        write_jpeg_comment(out.to_str().unwrap(), out2.to_str().unwrap(), b"short").unwrap();
+        assert_eq!(read_jpeg_comment(out2.to_str().unwrap()).unwrap(), b"short");
+
+        // image still decodes
+        assert!(image::image_dimensions(out.to_str().unwrap()).is_ok());
+
+        let _ = std::fs::remove_dir_all(&dir);
+    }
+
+    #[test]
+    fn distortion_precedence_and_legacy_import() {
+        use crate::panorama::Panorama;
+        use crate::profile;
+
+        let dir = std::env::temp_dir().join(format!("gipfel-prec-{}", std::process::id()));
+        let _ = std::fs::remove_dir_all(&dir);
+        let data = dir.join("data");
+        let home = dir.join("home");
+        std::fs::create_dir_all(&data).unwrap();
+
+        let old = (std::env::var("HOME").ok(), std::env::var("XDG_DATA_HOME").ok());
+        std::env::set_var("HOME", &home);
+        std::env::set_var("XDG_DATA_HOME", &data);
+
+        // legacy FLTK prefs file (C version's location and format)
+        let legacy = home
+            .join(".fltk")
+            .join("Johannes.HofmannATgmx.de")
+            .join("gipfel")
+            .join("DistortionProfiles.prefs");
+        std::fs::create_dir_all(legacy.parent().unwrap()).unwrap();
+        std::fs::write(
+            &legacy,
+            "; FLTK preferences file format 1.0\n; vendor: Johannes.HofmannATgmx.de\n\
+             ; application: gipfel/DistortionProfiles\n\n[.]\n\n\n\
+             [./Canon_EOS_5D_50.00_mm]\n\nk0:-0.01\nk1:0.001\nx0:3\n",
+        )
+        .unwrap();
+
+        // first access imports into the XDG data dir
+        assert!(!profile::prefs_path().exists());
+        let all = profile::list();
+        assert!(all.iter().any(|(n, p)| n == "Canon_EOS_5D_50.00_mm" && p.k0 == -0.01));
+        assert!(profile::prefs_path().starts_with(&data));
+        assert!(profile::prefs_path().exists());
+
+        let mut md = ImageMetaData::new();
+        md.manufacturer = Some("Canon".into());
+        md.model = Some("EOS_5D".into());
+        md.focal_length = 50.0;
+        let mut pan = Panorama::new();
+
+        // 1. comment wins when present
+        md.k0 = -0.9;
+        md.k1 = 0.0;
+        md.x0 = 0.0;
+        apply_image_metadata(&mut pan, &md, 1000, 800);
+        assert!((pan.parms.k0 + 0.9).abs() < 1e-12, "comment should win");
+
+        // 2. profile used when comment absent
+        md.k0 = f64::NAN;
+        md.k1 = f64::NAN;
+        md.x0 = f64::NAN;
+        apply_image_metadata(&mut pan, &md, 1000, 800);
+        assert!(
+            (pan.parms.k0 + 0.01).abs() < 1e-12,
+            "profile should apply: k0={}",
+            pan.parms.k0
+        );
+        // name is suggested regardless of the comment (C save dialog)
+        assert_eq!(
+            distortion_profile_name(&md).as_deref(),
+            Some("Canon_EOS_5D_50.00_mm")
+        );
+        md.k0 = 0.0;
+        assert_eq!(
+            distortion_profile_name(&md).as_deref(),
+            Some("Canon_EOS_5D_50.00_mm")
+        );
+
+        // 3. zero when neither
+        md.k0 = f64::NAN;
+        md.manufacturer = None;
+        apply_image_metadata(&mut pan, &md, 1000, 800);
+        assert_eq!(pan.parms.k0, 0.0);
+
+        match old.0 {
+            Some(h) => std::env::set_var("HOME", h),
+            None => std::env::remove_var("HOME"),
+        }
+        match old.1 {
+            Some(h) => std::env::set_var("XDG_DATA_HOME", h),
+            None => std::env::remove_var("XDG_DATA_HOME"),
+        }
+        let _ = std::fs::remove_dir_all(&dir);
+    }
+}
diff --git a/crates/gipfel-core/src/panorama.rs b/crates/gipfel-core/src/panorama.rs
new file mode 100644
index 0000000..9e27251
--- /dev/null
+++ b/crates/gipfel-core/src/panorama.rs
@@ -0,0 +1,502 @@
+//! Panorama: viewpoint, hill visibility (occlusion, refraction) and the
+//! link between geographic coordinates and image coordinates.
+//!
+//! Faithful port of gipfel's `Panorama.cxx`.
+
+use crate::hill::{Hill, DEG2RAD, PI};
+use crate::projection::{projection_from_type, Projection, ProjectionType, ViewParams};
+
+const EARTH_RADIUS: f64 = 6371000.785;
+
+pub struct Panorama {
+    pub hills: Vec<Hill>,
+    /// indices into `hills`, sorted by alpha descending (rebuilt by update_angles)
+    pub alpha_order: Vec<usize>,
+    /// indices into `hills`, sorted by alpha descending
+    pub close: Vec<usize>,
+    /// indices into `hills`
+    pub visible: Vec<usize>,
+
+    pub view_name: Option<String>,
+    /// latitude, radians
+    pub view_phi: f64,
+    /// longitude, radians
+    pub view_lam: f64,
+    /// meters
+    pub view_height: f64,
+
+    pub height_dist_ratio: f64,
+    pub hide_value: f64,
+    /// Additional visibility limit: maximum angular distance (radians) from
+    /// the image center. The app sets this to the angle the projection maps
+    /// to the image edge (the lens' actual field of view, including fitted
+    /// distortion), so that directions outside the photo — which a distortion
+    /// model can fold back into the image — neither get markers nor compass
+    /// ticks. 0 = no additional limit (only the projection's own limit).
+    pub max_view_angle: f64,
+
+    pub parms: ViewParams,
+    pub projection_type: ProjectionType,
+    proj: Box<dyn Projection>,
+}
+
+impl Default for Panorama {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl Panorama {
+    pub fn new() -> Self {
+        Panorama {
+            hills: Vec::new(),
+            alpha_order: Vec::new(),
+            close: Vec::new(),
+            visible: Vec::new(),
+            view_name: None,
+            view_phi: 0.0,
+            view_lam: 0.0,
+            view_height: 0.0,
+            height_dist_ratio: 0.07,
+            hide_value: 1.2,
+            max_view_angle: 0.0,
+            parms: ViewParams::new(),
+            projection_type: ProjectionType::Rectilinear,
+            proj: projection_from_type(ProjectionType::Rectilinear),
+        }
+    }
+
+    pub fn load_data(&mut self, path: &str) -> std::io::Result<()> {
+        let mut hills = crate::hill::load(path)?;
+        crate::hill::mark_duplicates(&mut hills, 0.00001);
+        self.hills = hills;
+        self.update_angles();
+        Ok(())
+    }
+
+    pub fn add_hills(&mut self, extra: &[Hill]) {
+        self.hills.extend_from_slice(extra);
+        crate::hill::mark_duplicates(&mut self.hills, 0.00001);
+        self.update_angles();
+    }
+
+    pub fn remove_hills(&mut self, flags: u32) {
+        self.hills.retain(|m| m.flags & flags == 0);
+        // indices shifted; rebuild all derived lists
+        self.alpha_order.clear();
+        self.close.clear();
+        self.visible.clear();
+    }
+
+    pub fn set_viewpoint(&mut self, idx: Option<usize>) {
+        let Some(idx) = idx else { return };
+        let m = &self.hills[idx];
+        self.view_phi = m.phi;
+        self.view_lam = m.lam;
+        self.view_height = m.height;
+        self.view_name = Some(m.name.clone());
+        self.update_angles();
+    }
+
+    /// Fit view parameters to the known hills (indices into `self.hills`).
+    /// `fit_distortion` enables the 7-parameter fit (k0/k1/x0) for 4+ hills;
+    /// otherwise the distortion params stay locked at their current values.
+    pub fn comp_params(&mut self, known: &[usize], fit_distortion: bool) -> Result<(), String> {
+        let ret = crate::lsq::comp_params(&*self.proj, &self.hills, known, &mut self.parms, fit_distortion);
+        if ret.is_ok() {
+            self.update_visible_mountains(Some(known));
+        }
+        ret
+    }
+
+    /// Like [`Panorama::comp_params`] but the observed hills are an external
+    /// slice (used by the golden fit tests, mirroring the C harness).
+    pub fn comp_params_with(&mut self, hills: &[Hill], known: &[usize], fit_distortion: bool) -> Result<(), String> {
+        let ret = crate::lsq::comp_params(&*self.proj, hills, known, &mut self.parms, fit_distortion);
+        if ret.is_ok() {
+            self.update_visible_mountains(None);
+        }
+        ret
+    }
+
+    pub fn set_center_angle(&mut self, deg: f64) {
+        self.parms.a_center = deg * DEG2RAD;
+        self.update_visible_mountains(None);
+    }
+
+    pub fn set_pitch_angle(&mut self, deg: f64) {
+        self.parms.a_pitch = deg * DEG2RAD;
+        self.update_coordinates(None);
+    }
+
+    pub fn set_roll_angle(&mut self, deg: f64) {
+        self.parms.a_roll = deg * DEG2RAD;
+        self.update_coordinates(None);
+    }
+
+    pub fn set_scale(&mut self, s: f64) {
+        self.parms.scale = s;
+        self.update_coordinates(None);
+    }
+
+    pub fn set_distortion_params(&mut self, k0: f64, k1: f64, x0: f64) {
+        self.parms.k0 = k0;
+        self.parms.k1 = k1;
+        self.parms.x0 = x0;
+        self.update_coordinates(None);
+    }
+
+    pub fn set_height_dist_ratio(&mut self, r: f64) {
+        self.height_dist_ratio = r;
+        self.update_close_mountains();
+    }
+
+    pub fn set_view_long(&mut self, deg: f64) {
+        self.view_lam = deg * DEG2RAD;
+        self.update_angles();
+    }
+
+    pub fn set_view_lat(&mut self, deg: f64) {
+        self.view_phi = deg * DEG2RAD;
+        self.update_angles();
+    }
+
+    pub fn set_view_height(&mut self, v: f64) {
+        self.view_height = v;
+        self.update_angles();
+    }
+
+    pub fn set_projection(&mut self, t: ProjectionType) {
+        self.projection_type = t;
+        self.proj = projection_from_type(t);
+        self.update_angles();
+    }
+
+    pub fn get_center_angle(&self) -> f64 {
+        self.parms.a_center / DEG2RAD
+    }
+    pub fn get_pitch_angle(&self) -> f64 {
+        self.parms.a_pitch / DEG2RAD
+    }
+    pub fn get_roll_angle(&self) -> f64 {
+        self.parms.a_roll / DEG2RAD
+    }
+    pub fn get_scale(&self) -> f64 {
+        self.parms.scale
+    }
+    pub fn get_view_long(&self) -> f64 {
+        self.view_lam / DEG2RAD
+    }
+    pub fn get_view_lat(&self) -> f64 {
+        self.view_phi / DEG2RAD
+    }
+    pub fn get_view_height(&self) -> f64 {
+        self.view_height
+    }
+
+    /// Last hill whose name matches exactly (C keeps the last match).
+    pub fn get_pos(&self, name: &str) -> Option<usize> {
+        let mut ret = None;
+        for (i, m) in self.hills.iter().enumerate() {
+            if m.name == name {
+                ret = Some(i);
+            }
+        }
+        ret
+    }
+
+    pub fn update_angles(&mut self) {
+        for i in 0..self.hills.len() {
+            let (phi, lam) = (self.hills[i].phi, self.hills[i].lam);
+            let d = self.distance(phi, lam);
+            self.hills[i].dist = d;
+            if phi != self.view_phi || lam != self.view_lam {
+                let a = self.alpha(i);
+                self.hills[i].alph = a;
+            }
+        }
+
+        // keep `hills` in insertion order (indices stay stable) and maintain a
+        // separate alpha-descending order, mirroring the C pointer-array sort
+        let mut idx: Vec<usize> = (0..self.hills.len()).collect();
+        idx.sort_by(|&a, &b| {
+            self.hills[b]
+                .alph
+                .partial_cmp(&self.hills[a].alph)
+                .unwrap_or(std::cmp::Ordering::Equal)
+        });
+        self.alpha_order = idx;
+
+        self.update_close_mountains();
+    }
+
+    /// Mark hills that are hidden behind a nearer, taller neighbour.
+    /// Faithful port including the original's `fabs(m->alph - n->alph > pi/2)`
+    /// oddity (comparison inside fabs, i.e. skip when m.alph - n.alph > pi/2).
+    pub fn mark_hidden(&mut self, hills: &[usize]) {
+        for i in 0..hills.len() {
+            let m = hills[i];
+            self.hills[m].flags &= !crate::hill::FLAG_HIDDEN;
+
+            if self.hills[m].flags & crate::hill::FLAG_DUPLIC != 0 {
+                continue;
+            }
+
+            for j in 0..hills.len() {
+                let n = hills[j];
+                if self.hills[n].flags & crate::hill::FLAG_DUPLIC != 0
+                    || self.hills[n].flags & crate::hill::FLAG_TRACK_POINT != 0
+                {
+                    continue;
+                }
+                if n == m || self.hills[m].alph - self.hills[n].alph > PI / 2.0 {
+                    continue;
+                }
+                if self.hills[m].dist < self.hills[n].dist
+                    || self.hills[m].a_elev > self.hills[n].a_elev
+                {
+                    continue;
+                }
+
+                let h = (self.hills[n].a_elev - self.hills[m].a_elev)
+                    / (self.hills[m].alph - self.hills[n].alph).abs();
+                if h.is_infinite() || h > self.hide_value {
+                    self.hills[m].flags |= crate::hill::FLAG_HIDDEN;
+                }
+            }
+        }
+    }
+
+    pub fn update_close_mountains(&mut self) {
+        self.close.clear();
+
+        let order = self.alpha_order.clone();
+        for &i in &order {
+            let m = &self.hills[i];
+            if m.flags & crate::hill::FLAG_TRACK_POINT != 0
+                || ((m.phi != self.view_phi || m.lam != self.view_lam)
+                    && (m.height / (m.dist * EARTH_RADIUS) > self.height_dist_ratio))
+            {
+                let elev = self.elevation(i);
+                self.hills[i].a_elev = elev;
+                self.close.push(i);
+            }
+        }
+
+        let close = self.close.clone();
+        self.mark_hidden(&close);
+        self.update_visible_mountains(None);
+    }
+
+    pub fn update_visible_mountains(&mut self, excluded: Option<&[usize]>) {
+        self.visible.clear();
+
+        for k in 0..self.close.len() {
+            let i = self.close[k];
+            if self.is_visible(self.hills[i].alph) {
+                self.visible.push(i);
+                self.hills[i].flags |= crate::hill::FLAG_VISIBLE;
+            } else {
+                self.hills[i].flags &= !crate::hill::FLAG_VISIBLE;
+            }
+        }
+
+        let owned;
+        let excl: Option<&[usize]> = match excluded {
+            Some(e) => {
+                owned = e.to_vec();
+                Some(&owned)
+            }
+            None => None,
+        };
+        self.update_coordinates(excl);
+    }
+
+    pub fn update_coordinates(&mut self, excluded: Option<&[usize]>) {
+        for k in 0..self.visible.len() {
+            let i = self.visible[k];
+            if let Some(e) = excluded {
+                if e.contains(&i) {
+                    continue;
+                }
+            }
+            let (alph, elev) = (self.hills[i].alph, self.hills[i].a_elev);
+            let (x, y) = self.proj.get_coordinates(alph, elev, &self.parms);
+            self.hills[i].x = x;
+            self.hills[i].y = y;
+        }
+    }
+
+    /// angular distance (fraction of earth radius) viewpoint -> (phi, lam)
+    pub fn distance(&self, phi: f64, lam: f64) -> f64 {
+        let d_lam = self.view_lam - lam;
+
+        let t1 = (phi.cos() * d_lam.sin()).powi(2);
+        let t2 = (self.view_phi.cos() * phi.sin()
+            - self.view_phi.sin() * phi.cos() * d_lam.cos())
+        .powi(2);
+
+        (t1 + t2).sqrt().atan2(
+            self.view_phi.sin() * phi.sin() + self.view_phi.cos() * phi.cos() * d_lam.cos(),
+        )
+    }
+
+    pub fn alpha(&self, idx: usize) -> f64 {
+        let m = &self.hills[idx];
+        let sin_alph = (m.lam - self.view_lam).sin() * m.phi.cos() / m.dist.sin();
+        let cos_alph =
+            (m.phi.sin() - self.view_phi.sin() * m.dist.cos()) / (self.view_phi.cos() * m.dist.sin());
+
+        (sin_alph.atan2(cos_alph) + 2.0 * PI) % (2.0 * PI)
+    }
+
+    /// approximation of refraction effect as described by Tom Chester at
+    /// <http://tchester.org/sgm/analysis/peaks/refraction.html>
+    pub fn refraction(&self, idx: usize) -> f64 {
+        let (alpha, t0) = (6.5, 10.0);
+
+        let a = 2.9e-4 * (-self.view_height / 10000.0).exp() / (1.0 + 2.9 * t0 / 760.0);
+        let b = 2.9 * alpha / (760.0 * (1.0 + 2.9 * t0 / 760.0));
+        let c = a * (b - 1.0 / 10.0);
+
+        c * self.get_real_distance(idx) / (2000.0 * (1.0 + a))
+    }
+
+    pub fn elevation(&self, idx: usize) -> f64 {
+        let m = &self.hills[idx];
+        let theta = self.refraction(idx);
+
+        let b = m.height + self.get_earth_radius(m.phi);
+        let c = self.view_height + self.get_earth_radius(self.view_phi);
+
+        ((m.dist.cos() * b - c) / (m.dist.sin() * b)).atan() - theta
+    }
+
+    /// local distance to center of WGS84 ellipsoid
+    pub fn get_earth_radius(&self, phi: f64) -> f64 {
+        let a = 6378137.000f64;
+        let b = 6356752.315f64;
+        let ata = phi.tan();
+
+        a * (ata.powi(2) + 1.0).powf(1.0 / 2.0)
+            * b.abs()
+            * (b.powi(2) + a.powi(2) * ata.powi(2)).powf(-1.0 / 2.0)
+    }
+
+    /// real (line-of-sight) distance in meters
+    pub fn get_real_distance(&self, idx: usize) -> f64 {
+        let m = &self.hills[idx];
+        let a = self.view_height + self.get_earth_radius(self.view_phi);
+        let b = m.height + self.get_earth_radius(m.phi);
+        let gam = m.dist;
+
+        (a.powi(2) + b.powi(2) - 2.0 * a * b * gam.cos()).sqrt()
+    }
+
+    pub fn is_visible(&self, alph: f64) -> bool {
+        let mut center_dist = (alph - self.parms.a_center).abs();
+        if center_dist > PI {
+            center_dist = 2.0 * PI - center_dist;
+        }
+        if center_dist >= self.proj.get_view_angle() {
+            return false;
+        }
+        self.max_view_angle <= 0.0 || center_dist < self.max_view_angle
+    }
+
+    /// Forward-project (azimuth, elevation) to image-centered pixel coords.
+    /// `None` when the direction is not visible or the projection is
+    /// undefined there (mirrors C `Panorama::get_coordinates` returning 1,
+    /// plus the NaN guard of C `ScanImage::get_pixel_nearest`).
+    pub fn project_point(&self, alph: f64, elev: f64) -> Option<(f64, f64)> {
+        if !self.is_visible(alph) {
+            return None;
+        }
+        let (x, y) = self.proj.get_coordinates(alph, elev, &self.parms);
+        if x.is_nan() || y.is_nan() {
+            return None;
+        }
+        Some((x, y))
+    }
+
+    /// Parameter-only snapshot (no hill data) for resampling: same viewpoint
+    /// and projection parameters, but without the UI's lens-FOV visibility
+    /// limit (the stitch output range defines the FOV itself).
+    pub fn params_snapshot(&self) -> Panorama {
+        let mut p = Panorama::new();
+        p.projection_type = self.projection_type;
+        p.proj = projection_from_type(self.projection_type);
+        p.parms = self.parms;
+        p.view_phi = self.view_phi;
+        p.view_lam = self.view_lam;
+        p.view_height = self.view_height;
+        p
+    }
+
+    pub fn get_coordinates(&self, alph: f64, elev: f64) -> Option<(f64, f64)> {
+        if self.is_visible(alph) {
+            Some(self.proj.get_coordinates(alph, elev, &self.parms))
+        } else {
+            None
+        }
+    }
+
+    /// Export rows for the visible hills, mirroring the C `export_hills`
+    /// output loop: `name\theight\tx\ty\tdistance` with image-center origin.
+    /// When `require_export` is set, only hills with `FLAG_EXPORT` are listed.
+    pub fn export_hill_rows(&self, img_w: i64, img_h: i64, require_export: bool) -> Vec<String> {
+        let mut rows = Vec::new();
+        for &i in &self.visible {
+            let m = &self.hills[i];
+            if m.flags & crate::hill::FLAG_DUPLIC != 0
+                || m.flags & crate::hill::FLAG_HIDDEN != 0
+                || (require_export && m.flags & crate::hill::FLAG_EXPORT == 0)
+            {
+                continue;
+            }
+
+            let x = m.x.round_ties_even() as i64 + img_w / 2;
+            let y = m.y.round_ties_even() as i64 + img_h / 2;
+            if x < 0 || x > img_w || y < 0 || y > img_h {
+                continue;
+            }
+
+            rows.push(format!(
+                "{}\t{}\t{}\t{}\t{}",
+                m.name,
+                m.height.round_ties_even() as i64,
+                x,
+                y,
+                self.get_real_distance(i).round_ties_even() as i64
+            ));
+        }
+        rows
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn max_view_angle_limits_visibility() {
+        let mut pan = Panorama::new();
+        pan.set_center_angle(0.0);
+        let a30 = 30f64.to_radians();
+        let a70 = 70f64.to_radians();
+
+        // rectilinear default limit is +-60 deg
+        assert!(pan.is_visible(a30));
+        assert!(!pan.is_visible(a70));
+
+        // a tighter lens FOV (say +-35 deg) hides the 40 deg direction that
+        // the projection's own +-60 deg limit would allow
+        pan.max_view_angle = 35f64.to_radians();
+        assert!(pan.is_visible(a30));
+        assert!(!pan.is_visible(40f64.to_radians()));
+
+        // a wider limit never overrides the projection's own bound
+        pan.max_view_angle = 80f64.to_radians();
+        assert!(!pan.is_visible(a70));
+    }
+}
diff --git a/crates/gipfel-core/src/profile.rs b/crates/gipfel-core/src/profile.rs
new file mode 100644
index 0000000..e8bda68
--- /dev/null
+++ b/crates/gipfel-core/src/profile.rs
@@ -0,0 +1,362 @@
+//! Lens distortion profiles.
+//!
+//! Profiles live in the standard per-user data directory
+//! (`$XDG_DATA_HOME/gipfel/distortion-profiles.prefs`, i.e.
+//! `~/.local/share/gipfel/…`). The file format is the FLTK preferences
+//! format used by the C version (`[./name]` groups, `key:value` entries,
+//! `%g` values), so the legacy file
+//! `~/.fltk/Johannes.HofmannATgmx.de/gipfel/DistortionProfiles.prefs`
+//! imports automatically into the new location on first use.
+//!
+//! File layout (as written by FLTK):
+//! ```text
+//! ; FLTK preferences file format 1.0
+//! ; vendor: Johannes.HofmannATgmx.de
+//! ; application: gipfel/DistortionProfiles
+//!
+//! [.]
+//!
+//! [./Canon_EOS_5D_50.00_mm]
+//!
+//! k0:-0.00123
+//! k1:2.5e-05
+//! x0:0.5
+//! ```
+
+use std::path::{Path, PathBuf};
+
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Profile {
+    pub k0: f64,
+    pub k1: f64,
+    pub x0: f64,
+}
+
+#[derive(Debug)]
+pub enum Error {
+    NotFound,
+    Exists,
+    BadName,
+    Io(std::io::Error),
+}
+
+impl std::fmt::Display for Error {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            Error::NotFound => write!(f, "profile not found"),
+            Error::Exists => write!(f, "profile already exists"),
+            Error::BadName => write!(f, "invalid profile name"),
+            Error::Io(e) => write!(f, "{e}"),
+        }
+    }
+}
+
+/// The profile store in the standard data directory
+/// (`$XDG_DATA_HOME/gipfel/distortion-profiles.prefs`).
+pub fn prefs_path() -> PathBuf {
+    dirs::data_dir()
+        .unwrap_or_else(|| PathBuf::from(std::env::var("HOME").unwrap_or_default()))
+        .join("gipfel")
+        .join("distortion-profiles.prefs")
+}
+
+/// Legacy FLTK preferences file used by the C version.
+fn fltk_prefs_path() -> PathBuf {
+    let home = std::env::var("HOME").unwrap_or_default();
+    Path::new(&home)
+        .join(".fltk")
+        .join("Johannes.HofmannATgmx.de")
+        .join("gipfel")
+        .join("DistortionProfiles.prefs")
+}
+
+/// Path to the profile store; imports legacy FLTK profiles on first use.
+fn ensure_store() -> PathBuf {
+    let path = prefs_path();
+    if !path.exists() {
+        let legacy = fltk_prefs_path();
+        if legacy.exists() {
+            if let Some(dir) = path.parent() {
+                let _ = std::fs::create_dir_all(dir);
+            }
+            // same format, so a copy suffices
+            let _ = std::fs::copy(&legacy, &path);
+        }
+    }
+    path
+}
+
+/// C `get_distortion_profile_name`: `"%s_%s_%.2f_mm"`. `None` when the
+/// camera is unknown or the focal length is unusable (C would build a name
+/// containing "nan" that can never match a stored profile).
+pub fn profile_name(manufacturer: &str, model: &str, focal_mm: f64) -> Option<String> {
+    if manufacturer.is_empty() || model.is_empty() || !focal_mm.is_finite() || focal_mm <= 0.0 {
+        return None;
+    }
+    let name = format!("{manufacturer}_{model}_{focal_mm:.2}_mm");
+    sanitize_name(&name).ok().map(|_| name)
+}
+
+/// Reject names that would break the prefs file structure (FLTK group paths
+/// use '/' as separator; control chars/newlines corrupt the file).
+pub fn sanitize_name(name: &str) -> Result<(), Error> {
+    if name.trim().is_empty() || name.len() > 255 {
+        return Err(Error::BadName);
+    }
+    if name
+        .chars()
+        .any(|c| c == '/' || c == '\\' || c.is_control())
+    {
+        return Err(Error::BadName);
+    }
+    Ok(())
+}
+
+/// Parse a prefs file into ordered (name, profile) pairs. Tolerant:
+/// comments (`;`), the root group `[.]`, blank lines and unknown keys are
+/// ignored; values are parsed with Rust's `f64` (superset of C `%g`).
+fn parse(text: &str) -> Vec<(String, Profile)> {
+    let mut out: Vec<(String, Profile)> = Vec::new();
+    let mut cur: Option<usize> = None;
+    for line in text.lines() {
+        let line = line.trim_end();
+        if line.is_empty() || line.starts_with(';') {
+            continue;
+        }
+        if let Some(group) = line.strip_prefix('[').and_then(|s| s.strip_suffix(']')) {
+            cur = match group.strip_prefix("./") {
+                Some(name) => {
+                    out.push((
+                        name.to_string(),
+                        Profile {
+                            k0: f64::NAN,
+                            k1: f64::NAN,
+                            x0: f64::NAN,
+                        },
+                    ));
+                    Some(out.len() - 1)
+                }
+                None => None, // "[.]" root group
+            };
+            continue;
+        }
+        let Some(idx) = cur else { continue };
+        if let Some((key, val)) = line.split_once(':') {
+            let Ok(v) = val.trim().parse::<f64>() else {
+                continue;
+            };
+            let p = &mut out[idx].1;
+            match key.trim() {
+                "k0" => p.k0 = v,
+                "k1" => p.k1 = v,
+                "x0" => p.x0 = v,
+                _ => {}
+            }
+        }
+    }
+    out.retain(|(_, p)| p.k0.is_finite() && p.k1.is_finite() && p.x0.is_finite());
+    out
+}
+
+/// Serialize in FLTK's layout so the C version can read our file.
+fn serialize(profiles: &[(String, Profile)]) -> String {
+    let mut s = String::from(
+        "; FLTK preferences file format 1.0\n\
+         ; vendor: Johannes.HofmannATgmx.de\n\
+         ; application: gipfel/DistortionProfiles\n\
+         \n\
+         [.]\n\
+         \n",
+    );
+    for (name, p) in profiles {
+        s.push_str(&format!(
+            "\n[./{name}]\n\nk0:{}\nk1:{}\nx0:{}\n",
+            fmt_g(p.k0),
+            fmt_g(p.k1),
+            fmt_g(p.x0)
+        ));
+    }
+    s
+}
+
+/// Closest Rust equivalent of C `%g`: scientific notation for magnitudes
+/// outside [1e-4, 1e15), exponent zero-padded to two digits ("2.5e-05").
+/// Unlike C we keep the shortest round-trip mantissa (more precise; the
+/// FLTK reader uses strtod and reads either fine).
+fn fmt_g(v: f64) -> String {
+    if v != 0.0 && (v.abs() < 1e-4 || v.abs() >= 1e15) {
+        if let Some((mant, exp)) = format!("{v:e}").split_once('e') {
+            let neg = exp.starts_with('-');
+            let digits = exp.trim_start_matches(['-', '+']);
+            let digits = if digits.len() < 2 { format!("0{digits}") } else { digits.to_string() };
+            let sign = if neg { '-' } else { '+' };
+            return format!("{mant}e{sign}{digits}");
+        }
+    }
+    if v == v.trunc() && v.abs() < 1e16 {
+        format!("{}", v as i64)
+    } else {
+        format!("{v}")
+    }
+}
+
+fn read_all(path: &Path) -> Result<Vec<(String, Profile)>, Error> {
+    match std::fs::read_to_string(path) {
+        Ok(text) => Ok(parse(&text)),
+        Err(e) if e.kind() == std::io::ErrorKind::NotFound => Ok(Vec::new()),
+        Err(e) => Err(Error::Io(e)),
+    }
+}
+
+fn write_all(path: &Path, profiles: &[(String, Profile)]) -> Result<(), Error> {
+    if let Some(dir) = path.parent() {
+        std::fs::create_dir_all(dir).map_err(Error::Io)?;
+    }
+    // atomic: temp file + rename
+    let tmp = path.with_extension("prefs.tmp");
+    std::fs::write(&tmp, serialize(profiles)).map_err(Error::Io)?;
+    std::fs::rename(&tmp, path).map_err(Error::Io)?;
+    Ok(())
+}
+
+/// Load a profile by exact name from the default store.
+pub fn load(name: &str) -> Result<Profile, Error> {
+    load_from(&ensure_store(), name)
+}
+
+pub fn load_from(path: &Path, name: &str) -> Result<Profile, Error> {
+    read_all(path)?
+        .into_iter()
+        .find(|(n, _)| n == name)
+        .map(|(_, p)| p)
+        .ok_or(Error::NotFound)
+}
+
+/// Save a profile. With `force == false` an existing profile of that name is
+/// not overwritten (C `save_distortion_params(name, force)`).
+pub fn save(name: &str, profile: Profile, force: bool) -> Result<(), Error> {
+    save_to(&ensure_store(), name, profile, force)
+}
+
+pub fn save_to(path: &Path, name: &str, profile: Profile, force: bool) -> Result<(), Error> {
+    sanitize_name(name)?;
+    let mut all = read_all(path)?;
+    if let Some(pos) = all.iter().position(|(n, _)| n == name) {
+        if !force {
+            return Err(Error::Exists);
+        }
+        all[pos].1 = profile;
+    } else {
+        all.push((name.to_string(), profile));
+    }
+    write_all(path, &all)
+}
+
+/// All profiles in the default store, in file order.
+pub fn list() -> Vec<(String, Profile)> {
+    read_all(&ensure_store()).unwrap_or_default()
+}
+
+pub fn list_from(path: &Path) -> Vec<(String, Profile)> {
+    read_all(path).unwrap_or_default()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Exact bytes as written by the C version (verified against FLTK 1.4).
+    const FLTK_FIXTURE: &str = "; FLTK preferences file format 1.0\n\
+         ; vendor: Johannes.HofmannATgmx.de\n\
+         ; application: gipfel/DistortionProfiles\n\
+         \n\
+         [.]\n\
+         \n\
+         \n\
+         [./Canon_EOS_5D_50.00_mm]\n\
+         \n\
+         k0:-0.00123\n\
+         k1:2.5e-05\n\
+         x0:0.5\n\
+         \n\
+         [./Nikon D850_35.00_mm]\n\
+         \n\
+         k0:0.01\n\
+         k1:0\n\
+         x0:0\n";
+
+    #[test]
+    fn parses_fltk_layout() {
+        let all = parse(FLTK_FIXTURE);
+        assert_eq!(all.len(), 2);
+        assert_eq!(all[0].0, "Canon_EOS_5D_50.00_mm");
+        assert_eq!(all[0].1.k0, -0.00123);
+        assert_eq!(all[0].1.k1, 2.5e-5);
+        assert_eq!(all[0].1.x0, 0.5);
+        assert_eq!(all[1].0, "Nikon D850_35.00_mm");
+        assert_eq!(all[1].1.k0, 0.01);
+        assert_eq!(all[1].1.k1, 0.0);
+    }
+
+    #[test]
+    fn serialize_replicates_fltk_layout() {
+        let all = parse(FLTK_FIXTURE);
+        assert_eq!(serialize(&all), FLTK_FIXTURE);
+    }
+
+    #[test]
+    fn fmt_g_matches_c_percent_g() {
+        assert_eq!(fmt_g(2.5e-5), "2.5e-05");
+        assert_eq!(fmt_g(-0.00123), "-0.00123");
+        assert_eq!(fmt_g(0.0001), "0.0001");
+        assert_eq!(fmt_g(1e15), "1e+15");
+        assert_eq!(fmt_g(0.0), "0");
+        assert_eq!(fmt_g(0.01), "0.01");
+    }
+
+    #[test]
+    fn profile_name_matches_c() {
+        assert_eq!(
+            profile_name("Canon", "EOS_5D", 50.0).as_deref(),
+            Some("Canon_EOS_5D_50.00_mm")
+        );
+        assert_eq!(
+            profile_name("Nikon", "D850", 34.666).as_deref(),
+            Some("Nikon_D850_34.67_mm")
+        );
+        assert!(profile_name("", "D850", 35.0).is_none());
+        assert!(profile_name("Canon", "D850", f64::NAN).is_none());
+        assert!(profile_name("Can/on", "D850", 35.0).is_none());
+    }
+
+    #[test]
+    fn roundtrip_and_overwrite() {
+        let dir = std::env::temp_dir().join(format!("gipfel-prof-{}", std::process::id()));
+        let path = dir.join("DistortionProfiles.prefs");
+        let _ = std::fs::remove_dir_all(&dir);
+
+        let p = Profile { k0: -0.001, k1: 2e-5, x0: 0.0 };
+        save_to(&path, "Cam_A_24.00_mm", p, false).unwrap();
+        // exists -> error without force
+        assert!(matches!(
+            save_to(&path, "Cam_A_24.00_mm", p, false),
+            Err(Error::Exists)
+        ));
+        save_to(&path, "Cam_B_50.00_mm", Profile { k0: 0.1, k1: 0.0, x0: 1.0 }, false).unwrap();
+        // overwrite first, keep second
+        let p2 = Profile { k0: -0.002, k1: 2e-5, x0: 0.0 };
+        save_to(&path, "Cam_A_24.00_mm", p2, true).unwrap();
+
+        let all = list_from(&path);
+        assert_eq!(all.len(), 2);
+        assert_eq!(all[0].1, p2);
+        assert_eq!(load_from(&path, "Cam_B_50.00_mm").unwrap().k0, 0.1);
+        assert!(matches!(load_from(&path, "nope"), Err(Error::NotFound)));
+
+        // bad names rejected
+        assert!(save_to(&path, "a/b", p, true).is_err());
+        assert!(save_to(&path, "with\nnewline", p, true).is_err());
+
+        std::fs::remove_dir_all(&dir).unwrap();
+    }
+}
diff --git a/crates/gipfel-core/src/projection/cylindrical.rs b/crates/gipfel-core/src/projection/cylindrical.rs
new file mode 100644
index 0000000..9adf79c
--- /dev/null
+++ b/crates/gipfel-core/src/projection/cylindrical.rs
@@ -0,0 +1,86 @@
+//! Cylindrical projection.
+
+mod funcs {
+    #![allow(non_snake_case, unused_parens, unused_variables, dead_code, clippy::all)]
+
+    // math shims: the maxima-generated code uses C-style calls
+    #[inline(always)]
+    fn sin(x: f64) -> f64 {
+        x.sin()
+    }
+    #[inline(always)]
+    fn cos(x: f64) -> f64 {
+        x.cos()
+    }
+    #[inline(always)]
+    fn tan(x: f64) -> f64 {
+        x.tan()
+    }
+    #[inline(always)]
+    fn sqrt(x: f64) -> f64 {
+        x.sqrt()
+    }
+    #[inline(always)]
+    fn pow(x: f64, y: f64) -> f64 {
+        x.powf(y)
+    }
+    #[inline(always)]
+    fn exp(x: f64) -> f64 {
+        x.exp()
+    }
+    #[inline(always)]
+    fn fabs(x: f64) -> f64 {
+        x.abs()
+    }
+    #[inline(always)]
+    fn atan2(y: f64, x: f64) -> f64 {
+        y.atan2(x)
+    }
+    #[inline(always)]
+    fn sec(x: f64) -> f64 {
+        1.0 / x.cos()
+    }
+
+    include!("generated/cylindrical_funcs.rs");
+}
+
+use super::{ProjArgs, Projection};
+
+pub struct ProjectionCylindrical;
+
+impl Projection for ProjectionCylindrical {
+    fn get_view_angle(&self) -> f64 {
+        6.2831853 /* 360 deg */
+    }
+
+    fn mac_x(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dc_view(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dc_view(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dc_pitch(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dc_pitch(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dc_roll(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dc_roll(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dscale(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dscale(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dc_view(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dc_view(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dc_pitch(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dc_pitch(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dc_roll(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dc_roll(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dscale(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dscale(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+}
diff --git a/crates/gipfel-core/src/projection/generated/cylindrical_funcs.rs b/crates/gipfel-core/src/projection/generated/cylindrical_funcs.rs
new file mode 100644
index 0000000..220a7d8
--- /dev/null
+++ b/crates/gipfel-core/src/projection/generated/cylindrical_funcs.rs
@@ -0,0 +1,10 @@
+pub fn mac_x ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return (cos(c_roll)*(m_view-c_view)-(sin(c_roll)*(_2-_1))/(_1*_2+1.0))*scale ; } 
+pub fn mac_y ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return ((-sin(c_roll)*(m_view-c_view))-(cos(c_roll)*(_2-_1))/(_1*_2+1.0))*scale ; } 
+pub fn mac_x_dc_view ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { return -cos(c_roll)*scale ; } 
+pub fn mac_x_dc_pitch ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; _1 = pow(sec(c_pitch),2.0) ; _2 = sin(c_roll) ; _3 = tan(m_elev) ; _4 = tan(c_pitch) ; _5 = _4*_3+1.0 ; return ((_1*_2*_3*(_3-_4))/pow(_5,2.0)+(_1*_2)/_5)*scale ; } 
+pub fn mac_x_dc_roll ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return ((-sin(c_roll)*(m_view-c_view))-(cos(c_roll)*(_2-_1))/(_1*_2+1.0))*scale ; } 
+pub fn mac_x_dscale ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return cos(c_roll)*(m_view-c_view)-(sin(c_roll)*(_2-_1))/(_1*_2+1.0) ; } 
+pub fn mac_y_dc_view ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { return sin(c_roll)*scale ; } 
+pub fn mac_y_dc_pitch ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; _1 = pow(sec(c_pitch),2.0) ; _2 = cos(c_roll) ; _3 = tan(m_elev) ; _4 = tan(c_pitch) ; _5 = _4*_3+1.0 ; return ((_1*_2*_3*(_3-_4))/pow(_5,2.0)+(_1*_2)/_5)*scale ; } 
+pub fn mac_y_dc_roll ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return ((sin(c_roll)*(_2-_1))/(_1*_2+1.0)-cos(c_roll)*(m_view-c_view))*scale ; } 
+pub fn mac_y_dscale ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; _1 = tan(c_pitch) ; _2 = tan(m_elev) ; return (-sin(c_roll)*(m_view-c_view))-(cos(c_roll)*(_2-_1))/(_1*_2+1.0) ; } 
diff --git a/crates/gipfel-core/src/projection/generated/rectilinear_funcs.rs b/crates/gipfel-core/src/projection/generated/rectilinear_funcs.rs
new file mode 100644
index 0000000..7bd450d
--- /dev/null
+++ b/crates/gipfel-core/src/projection/generated/rectilinear_funcs.rs
@@ -0,0 +1,16 @@
+pub fn mac_x ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; _16 = pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0) ; return scale*_15*(k0*sqrt(_16)+k1*_16+1.0) ; } 
+pub fn mac_y ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; _16 = pow(x0+_1*_5*_14*_11+_13*_11*_12,2.0)+pow(_15,2.0) ; return _15*scale*(k0*sqrt(_16)+k1*_16+1.0) ; } 
+pub fn mac_x_dc_view ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; let mut _21 : f64 = 0.0; let mut _22 : f64 = 0.0; let mut _23 : f64 = 0.0; let mut _24 : f64 = 0.0; let mut _25 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = _4*_5*_10+_2*_3 ; _12 = 1.0/_11 ; _13 = _2*_5*_10-_4*_3 ; _14 = cos(c_roll) ; _15 = _6*_9-_8*_7 ; _16 = x0+_14*_5*_15*_12+_1*_12*_13 ; _17 = _14*_12*_13-_1*_5*_15*_12 ; _18 = 1.0/pow(_11,2.0) ; _19 = pow(_5,2.0) ; _20 = pow(_15,2.0) ; _21 = (-_8*_9)-_6*_7 ; _22 = _14*_5*_21*_12+_2*_1*_5*_15*_12-_4*_14*_19*_20*_18-_4*_1*_5*_15*_18*_13 ; _23 = 2.0*_22*_16+2.0*_17*((-_1*_5*_21*_12)+_2*_14*_5*_15*_12+_4*_1*_19*_20*_18-_4*_14*_5*_15*_18*_13) ; _24 = pow(_17,2.0)+pow(_16,2.0) ; _25 = sqrt(_24) ; return _22*scale*(k0*_25+k1*_24+1.0)+scale*_16*((k0*_23)/(2.0*_25)+k1*_23) ; } 
+pub fn mac_x_dc_pitch ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; let mut _21 : f64 = 0.0; let mut _22 : f64 = 0.0; let mut _23 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = _4*_5*_10+_2*_3 ; _12 = 1.0/_11 ; _13 = _2*_5*_10-_4*_3 ; _14 = cos(c_roll) ; _15 = _6*_9-_8*_7 ; _16 = x0+_14*_5*_15*_12+_1*_12*_13 ; _17 = 1.0/pow(_11,2.0) ; _18 = _4*_3-_2*_5*_10 ; _19 = _14*_12*_13-_1*_5*_15*_12 ; _20 = (-_1*_17*_18*_13)-_14*_5*_15*_17*_18+_1 ; _21 = 2.0*_20*_16+2.0*((-_14*_17*_18*_13)+_1*_5*_15*_17*_18+_14)*_19 ; _22 = pow(_19,2.0)+pow(_16,2.0) ; _23 = sqrt(_22) ; return _20*scale*(k0*_23+k1*_22+1.0)+scale*_16*((k0*_21)/(2.0*_23)+k1*_21) ; } 
+pub fn mac_x_dc_roll ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; _16 = _13*_11*_12-_1*_5*_14*_11 ; _17 = 2.0*_16*_15+2.0*((-_13*_5*_14*_11)-_1*_11*_12)*_16 ; _18 = pow(_16,2.0)+pow(_15,2.0) ; _19 = sqrt(_18) ; return _16*scale*(k0*_19+k1*_18+1.0)+scale*_15*((k0*_17)/(2.0*_19)+k1*_17) ; } 
+pub fn mac_x_dscale ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; _16 = pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0) ; return _15*(k0*sqrt(_16)+k1*_16+1.0) ; } 
+pub fn mac_x_dk0 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; return _15*sqrt(pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0))*scale ; } 
+pub fn mac_x_dk1 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; return _15*(pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0))*scale ; } 
+pub fn mac_x_dx0 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; _1 = sin(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = cos(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = x0+_13*_5*_14*_11+_1*_11*_12 ; _16 = pow(_15,2.0)+pow(_13*_11*_12-_1*_5*_14*_11,2.0) ; _17 = sqrt(_16) ; return scale*(k0*_17+k1*_16+1.0)+scale*_15*((k0*_15)/_17+2.0*k1*_15) ; } 
+pub fn mac_y_dc_view ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; let mut _21 : f64 = 0.0; let mut _22 : f64 = 0.0; let mut _23 : f64 = 0.0; let mut _24 : f64 = 0.0; let mut _25 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = _4*_5*_10+_2*_3 ; _12 = 1.0/_11 ; _13 = _2*_5*_10-_4*_3 ; _14 = sin(c_roll) ; _15 = _6*_9-_8*_7 ; _16 = _1*_12*_13-_14*_5*_15*_12 ; _17 = 1.0/pow(_11,2.0) ; _18 = pow(_5,2.0) ; _19 = pow(_15,2.0) ; _20 = (-_8*_9)-_6*_7 ; _21 = (-_14*_5*_20*_12)+_2*_1*_5*_15*_12+_4*_14*_18*_19*_17-_4*_1*_5*_15*_17*_13 ; _22 = x0+_1*_5*_15*_12+_14*_12*_13 ; _23 = 2.0*(_1*_5*_20*_12+_2*_14*_5*_15*_12-_4*_1*_18*_19*_17-_4*_14*_5*_15*_17*_13)*_22+2.0*_16*_21 ; _24 = pow(_22,2.0)+pow(_16,2.0) ; _25 = sqrt(_24) ; return _21*scale*(k0*_25+k1*_24+1.0)+_16*scale*((k0*_23)/(2.0*_25)+k1*_23) ; } 
+pub fn mac_y_dc_pitch ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; let mut _21 : f64 = 0.0; let mut _22 : f64 = 0.0; let mut _23 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = _4*_5*_10+_2*_3 ; _12 = 1.0/_11 ; _13 = _2*_5*_10-_4*_3 ; _14 = sin(c_roll) ; _15 = _6*_9-_8*_7 ; _16 = _1*_12*_13-_14*_5*_15*_12 ; _17 = 1.0/pow(_11,2.0) ; _18 = _4*_3-_2*_5*_10 ; _19 = (-_1*_17*_18*_13)+_14*_5*_15*_17*_18+_1 ; _20 = x0+_1*_5*_15*_12+_14*_12*_13 ; _21 = 2.0*((-_14*_17*_18*_13)-_1*_5*_15*_17*_18+_14)*_20+2.0*_19*_16 ; _22 = pow(_20,2.0)+pow(_16,2.0) ; _23 = sqrt(_22) ; return _19*scale*(k0*_23+k1*_22+1.0)+_16*scale*((k0*_21)/(2.0*_23)+k1*_21) ; } 
+pub fn mac_y_dc_roll ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; let mut _17 : f64 = 0.0; let mut _18 : f64 = 0.0; let mut _19 : f64 = 0.0; let mut _20 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; _16 = (-_1*_5*_14*_11)-_13*_11*_12 ; _17 = x0+_1*_5*_14*_11+_13*_11*_12 ; _18 = 2.0*_15*_17+2.0*_16*_15 ; _19 = pow(_17,2.0)+pow(_15,2.0) ; _20 = sqrt(_19) ; return _16*scale*(k0*_20+k1*_19+1.0)+_15*scale*((k0*_18)/(2.0*_20)+k1*_18) ; } 
+pub fn mac_y_dscale ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; _16 = pow(x0+_1*_5*_14*_11+_13*_11*_12,2.0)+pow(_15,2.0) ; return _15*(k0*sqrt(_16)+k1*_16+1.0) ; } 
+pub fn mac_y_dk0 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; return _15*scale*sqrt(pow(x0+_1*_5*_14*_11+_13*_11*_12,2.0)+pow(_15,2.0)) ; } 
+pub fn mac_y_dk1 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; return _15*scale*(pow(x0+_1*_5*_14*_11+_13*_11*_12,2.0)+pow(_15,2.0)) ; } 
+pub fn mac_y_dx0 ( c_view: f64, c_pitch: f64, c_roll: f64, scale: f64, k0: f64, k1: f64, x0: f64, m_view: f64, m_elev: f64 ) -> f64 { let mut _1 : f64 = 0.0; let mut _2 : f64 = 0.0; let mut _3 : f64 = 0.0; let mut _4 : f64 = 0.0; let mut _5 : f64 = 0.0; let mut _6 : f64 = 0.0; let mut _7 : f64 = 0.0; let mut _8 : f64 = 0.0; let mut _9 : f64 = 0.0; let mut _10 : f64 = 0.0; let mut _11 : f64 = 0.0; let mut _12 : f64 = 0.0; let mut _13 : f64 = 0.0; let mut _14 : f64 = 0.0; let mut _15 : f64 = 0.0; let mut _16 : f64 = 0.0; _1 = cos(c_roll) ; _2 = sin(c_pitch) ; _3 = sin(m_elev) ; _4 = cos(c_pitch) ; _5 = cos(m_elev) ; _6 = cos(c_view) ; _7 = cos(m_view) ; _8 = sin(c_view) ; _9 = sin(m_view) ; _10 = _8*_9+_6*_7 ; _11 = 1.0/(_4*_5*_10+_2*_3) ; _12 = _2*_5*_10-_4*_3 ; _13 = sin(c_roll) ; _14 = _6*_9-_8*_7 ; _15 = _1*_11*_12-_13*_5*_14*_11 ; _16 = x0+_1*_5*_14*_11+_13*_11*_12 ; return _15*((_16*k0)/sqrt(pow(_16,2.0)+pow(_15,2.0))+2.0*k1*_16)*scale ; } 
diff --git a/crates/gipfel-core/src/projection/mod.rs b/crates/gipfel-core/src/projection/mod.rs
new file mode 100644
index 0000000..65b74e5
--- /dev/null
+++ b/crates/gipfel-core/src/projection/mod.rs
@@ -0,0 +1,174 @@
+//! Projections: image coordinates as a function of azimuth/elevation angle,
+//! with symbolically pre-generated derivatives (see `codegen/`).
+//!
+//! Faithful port of gipfel's `ProjectionLSQ.{H,cxx}`, `ProjectionRectilinear.*`
+//! and `ProjectionCylindrical.*`.
+
+pub mod cylindrical;
+pub mod rectilinear;
+
+pub const NAN: f64 = f64::NAN;
+
+/// The 7 fitted view parameters (radians where angular).
+#[derive(Clone, Copy, Debug, Default, PartialEq)]
+pub struct ViewParams {
+    /// azimuth of image center, radians
+    pub a_center: f64,
+    /// focal scale (pixels)
+    pub scale: f64,
+    /// pitch angle of image center (camera nose up), radians
+    pub a_pitch: f64,
+    /// camera roll, radians
+    pub a_roll: f64,
+    /// radial distortion coefficient 0
+    pub k0: f64,
+    /// radial distortion coefficient 1
+    pub k1: f64,
+    /// distortion center x, pixels
+    pub x0: f64,
+}
+
+impl ViewParams {
+    pub fn new() -> Self {
+        ViewParams {
+            a_center: 0.0,
+            scale: 3500.0,
+            a_pitch: 0.0,
+            a_roll: 0.0,
+            k0: 0.0,
+            k1: 0.0,
+            x0: 0.0,
+        }
+    }
+}
+
+/// The 9 arguments every generated macro function takes.
+#[derive(Clone, Copy, Debug)]
+pub struct ProjArgs {
+    pub c_view: f64,
+    pub c_pitch: f64,
+    pub c_roll: f64,
+    pub scale: f64,
+    pub k0: f64,
+    pub k1: f64,
+    pub x0: f64,
+    pub m_view: f64,
+    pub m_elev: f64,
+}
+
+impl ProjArgs {
+    pub fn from_params(p: &ViewParams, m_view: f64, m_elev: f64) -> Self {
+        ProjArgs {
+            c_view: p.a_center,
+            c_pitch: p.a_pitch,
+            c_roll: p.a_roll,
+            scale: p.scale,
+            k0: p.k0,
+            k1: p.k1,
+            x0: p.x0,
+            m_view,
+            m_elev,
+        }
+    }
+}
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum ProjectionType {
+    Rectilinear = 0,
+    Cylindrical = 1,
+}
+
+/// A projection plus its symbolic partial derivatives w.r.t. the view
+/// parameters. Methods default to `NAN`, matching the C++ base class.
+pub trait Projection: Send + Sync {
+    /// half-width of the visible azimuth range, radians
+    fn get_view_angle(&self) -> f64;
+
+    fn mac_x(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_y(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_x_dc_view(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_x_dc_pitch(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_x_dc_roll(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_x_dscale(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_x_dk0(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_x_dk1(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_x_dx0(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_y_dc_view(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_y_dc_pitch(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_y_dc_roll(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_y_dscale(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_y_dk0(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_y_dk1(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+    fn mac_y_dx0(&self, a: ProjArgs) -> f64 {
+        let _ = a;
+        NAN
+    }
+
+    /// Project azimuth/elevation to image coordinates (relative to image center),
+    /// normalizing azimuth around the view center like the C++ original.
+    fn get_coordinates(&self, alph: f64, elev: f64, parms: &ViewParams) -> (f64, f64) {
+        let mut alph = alph;
+        let pi = crate::hill::PI;
+
+        if alph - parms.a_center > pi {
+            alph -= 2.0 * pi;
+        } else if alph - parms.a_center < -pi {
+            alph += 2.0 * pi;
+        }
+
+        let args = ProjArgs::from_params(parms, alph, elev);
+        (self.mac_x(args), self.mac_y(args))
+    }
+}
+
+pub fn projection_from_type(t: ProjectionType) -> Box<dyn Projection> {
+    match t {
+        ProjectionType::Rectilinear => Box::new(rectilinear::ProjectionRectilinear),
+        ProjectionType::Cylindrical => Box::new(cylindrical::ProjectionCylindrical),
+    }
+}
diff --git a/crates/gipfel-core/src/projection/rectilinear.rs b/crates/gipfel-core/src/projection/rectilinear.rs
new file mode 100644
index 0000000..a7037d6
--- /dev/null
+++ b/crates/gipfel-core/src/projection/rectilinear.rs
@@ -0,0 +1,104 @@
+//! Rectilinear (gnomonic) projection with lens distortion.
+
+mod funcs {
+    #![allow(non_snake_case, unused_parens, unused_variables, dead_code, clippy::all)]
+
+    // math shims: the maxima-generated code uses C-style calls
+    #[inline(always)]
+    fn sin(x: f64) -> f64 {
+        x.sin()
+    }
+    #[inline(always)]
+    fn cos(x: f64) -> f64 {
+        x.cos()
+    }
+    #[inline(always)]
+    fn tan(x: f64) -> f64 {
+        x.tan()
+    }
+    #[inline(always)]
+    fn sqrt(x: f64) -> f64 {
+        x.sqrt()
+    }
+    #[inline(always)]
+    fn pow(x: f64, y: f64) -> f64 {
+        x.powf(y)
+    }
+    #[inline(always)]
+    fn exp(x: f64) -> f64 {
+        x.exp()
+    }
+    #[inline(always)]
+    fn fabs(x: f64) -> f64 {
+        x.abs()
+    }
+    #[inline(always)]
+    fn atan2(y: f64, x: f64) -> f64 {
+        y.atan2(x)
+    }
+    #[inline(always)]
+    fn sec(x: f64) -> f64 {
+        1.0 / x.cos()
+    }
+
+    include!("generated/rectilinear_funcs.rs");
+}
+
+use super::{ProjArgs, Projection};
+
+pub struct ProjectionRectilinear;
+
+impl Projection for ProjectionRectilinear {
+    fn get_view_angle(&self) -> f64 {
+        1.0471976 /* 60 deg */
+    }
+
+    fn mac_x(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dc_view(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dc_view(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dc_pitch(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dc_pitch(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dc_roll(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dc_roll(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dscale(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dscale(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dk0(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dk0(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dk1(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dk1(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_x_dx0(&self, a: ProjArgs) -> f64 {
+        funcs::mac_x_dx0(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dc_view(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dc_view(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dc_pitch(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dc_pitch(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dc_roll(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dc_roll(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dscale(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dscale(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dk0(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dk0(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dk1(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dk1(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+    fn mac_y_dx0(&self, a: ProjArgs) -> f64 {
+        funcs::mac_y_dx0(a.c_view, a.c_pitch, a.c_roll, a.scale, a.k0, a.k1, a.x0, a.m_view, a.m_elev)
+    }
+}
diff --git a/crates/gipfel-core/src/stitch.rs b/crates/gipfel-core/src/stitch.rs
new file mode 100644
index 0000000..a311961
--- /dev/null
+++ b/crates/gipfel-core/src/stitch.rs
@@ -0,0 +1,323 @@
+//! Cylindrical panorama resampling — port of the C++ `Stitch.cxx` /
+//! `ScanImage.cxx` stitch mode.
+//!
+//! The output is a cylindrical (equirectangular-in-azimuth, tan-in-elevation)
+//! image spanning `view_start..view_end` degrees of azimuth. Each output
+//! pixel is mapped back through every source image's fitted projection; the
+//! first source that covers the direction contributes its (resampled) color.
+//!
+//! All sampling arithmetic reproduces the C implementation bit-for-bit:
+//! colors are handled on the 0..=65025 scale (`value * 255`), nearest uses
+//! `rint` (round-half-to-even), bicubic uses the same Keys cubic and only
+//! rounds after the final vertical interpolation.
+
+use crate::metadata::{ImageMetaData, apply_image_metadata};
+use crate::panorama::Panorama;
+use image::RgbImage;
+
+#[derive(Clone, Copy, PartialEq, Eq, Debug)]
+pub enum ScanMode {
+    Nearest,
+    Bicubic,
+}
+
+/// C MAX_VALUE: clamp bound on the 0..=65025 (value*255) scale.
+pub const MAX_VALUE: i64 = 65025;
+
+fn pixel(img: &RgbImage, x: i64, y: i64) -> Option<[i64; 3]> {
+    if x < 0 || y < 0 || x >= img.width() as i64 || y >= img.height() as i64 {
+        return None;
+    }
+    let p = img.get_pixel(x as u32, y as u32);
+    Some([p[0] as i64 * 255, p[1] as i64 * 255, p[2] as i64 * 255])
+}
+
+/// C ScanImage::get_pixel_nearest: NaN-guard, then round to nearest pixel.
+fn get_pixel_nearest(img: &RgbImage, x: f64, y: f64) -> Option<[i64; 3]> {
+    if x.is_nan() || y.is_nan() {
+        return None;
+    }
+    pixel(img, x.round_ties_even() as i64, y.round_ties_even() as i64)
+}
+
+/// C ScanImage interp_cubic (Keys cubic convolution kernel, evaluated in
+/// Bernstein-like power basis).
+#[inline]
+fn interp_cubic(x: f64, x2: f64, x3: f64, v: &[f64; 4]) -> f64 {
+    let a0 = v[3] - v[2] - v[0] + v[1];
+    let a1 = v[0] - v[1] - a0;
+    let a2 = v[2] - v[0];
+    let a3 = v[1];
+    a0 * x3 + a1 * x2 + a2 * x + a3
+}
+
+/// C ScanImage::get_pixel_bicubic. Any of the 16 taps outside the source
+/// makes the pixel unset (C returns 1), matching the "leave black" behavior.
+fn get_pixel_bicubic(img: &RgbImage, x: f64, y: f64) -> Option<[i64; 3]> {
+    if x.is_nan() || y.is_nan() {
+        return None;
+    }
+    let fl_x = x.floor();
+    let fl_y = y.floor();
+    let dx = x - fl_x;
+    let dx2 = dx * dx;
+    let dx3 = dx2 * dx;
+    let dy = y - fl_y;
+    let dy2 = dy * dy;
+    let dy3 = dy2 * dy;
+
+    let mut c1 = [[0f64; 4]; 3];
+    for iy in 0..4usize {
+        let mut c = [[0f64; 4]; 3];
+        for ix in 0..4usize {
+            let p = pixel(img, fl_x as i64 + ix as i64 - 1, fl_y as i64 + iy as i64 - 1)?;
+            for l in 0..3 {
+                c[l][ix] = p[l] as f64;
+            }
+        }
+        for l in 0..3 {
+            c1[l][iy] = interp_cubic(dx, dx2, dx3, &c[l]);
+        }
+    }
+    Some([
+        interp_cubic(dy, dy2, dy3, &c1[0]).round_ties_even() as i64,
+        interp_cubic(dy, dy2, dy3, &c1[1]).round_ties_even() as i64,
+        interp_cubic(dy, dy2, dy3, &c1[2]).round_ties_even() as i64,
+    ])
+}
+
+/// C ScanImage::get_pixel(mode, x, y). Values on the 0..=65025 scale;
+/// bicubic may overshoot (clamped by the caller, as in C).
+pub fn get_pixel(img: &RgbImage, mode: ScanMode, x: f64, y: f64) -> Option<[i64; 3]> {
+    match mode {
+        ScanMode::Nearest => get_pixel_nearest(img, x, y),
+        ScanMode::Bicubic => get_pixel_bicubic(img, x, y),
+    }
+}
+
+/// One input photo: its fitted projection plus its decoded pixels.
+pub struct StitchSource {
+    pub pan: Panorama,
+    pub img: RgbImage,
+}
+
+impl StitchSource {
+    /// Load a JPEG and apply its gipfel metadata (comment or EXIF), like
+    /// `GipfelWidget::load_image` does for the stitch mode.
+    pub fn load(path: &str) -> Result<Self, String> {
+        let img = image::open(path)
+            .map(|i| i.to_rgb8())
+            .map_err(|e| format!("{}: {e}", path))?;
+        let mut md = ImageMetaData::new();
+        md.load_image(path);
+        let mut pan = Panorama::new();
+        apply_image_metadata(&mut pan, &md, img.width(), img.height());
+        Ok(Self { pan, img })
+    }
+
+    /// C GipfelWidget::get_pixel: forward-project (azimuth, elevation)
+    /// through this image's projection and sample the source.
+    pub fn get_pixel(&self, mode: ScanMode, a_alph: f64, a_elev: f64) -> Option<[i64; 3]> {
+        let (px, py) = self.pan.project_point(a_alph, a_elev)?;
+        get_pixel(
+            &self.img,
+            mode,
+            px + self.img.width() as f64 / 2.0,
+            py + self.img.height() as f64 / 2.0,
+        )
+    }
+}
+
+/// The resampled cylindrical panorama. Pixel values are RGB triplets on the
+/// C 0..=65025 scale (clamped); `hit` marks pixels some source contributed
+/// to (unset pixels are black in the C output).
+pub struct StitchedImage {
+    pub w: u32,
+    pub h: u32,
+    pub pixels: Vec<[u16; 3]>,
+    pub hit: Vec<bool>,
+}
+
+/// Port of `Stitch::resample`: output column x maps linearly to azimuth,
+/// output row y maps to elevation via `atan((y_off - y) / radius)`.
+pub fn resample(
+    sources: &[StitchSource],
+    mode: ScanMode,
+    w: u32,
+    h: u32,
+    view_start_deg: f64,
+    view_end_deg: f64,
+) -> StitchedImage {
+    let view_start = view_start_deg.to_radians();
+    let view_end = view_end_deg.to_radians();
+    let mut pixels = vec![[0u16; 3]; (w as usize) * (h as usize)];
+    let mut hit = vec![false; (w as usize) * (h as usize)];
+    if !(view_end > view_start) || w == 0 || h == 0 {
+        return StitchedImage {
+            w,
+            h,
+            pixels,
+            hit,
+        };
+    }
+    let step_view = (view_end - view_start) / w as f64;
+    let y_off = (h / 2) as f64;
+    let radius = w as f64 / (view_end - view_start);
+
+    for y in 0..h {
+        let a_elev = ((y_off - y as f64) / radius).atan();
+        for x in 0..w {
+            let a_view = view_start + x as f64 * step_view;
+            for src in sources {
+                let Some(rgb) = src.get_pixel(mode, a_view, a_elev) else {
+                    continue;
+                };
+                let clamped = [
+                    rgb[0].clamp(0, MAX_VALUE) as u16,
+                    rgb[1].clamp(0, MAX_VALUE) as u16,
+                    rgb[2].clamp(0, MAX_VALUE) as u16,
+                ];
+                let i = (y as usize) * (w as usize) + (x as usize);
+                pixels[i] = clamped;
+                hit[i] = true;
+                break; // C: first contributing image wins
+            }
+        }
+    }
+
+    StitchedImage {
+        w,
+        h,
+        pixels,
+        hit,
+    }
+}
+
+/// Write an RGBA TIFF like the C `TIFFOutputImage`: 8 bit (value/255, alpha
+/// 255) or 16 bit (raw 0..=65025, alpha 65025), deflate-compressed.
+pub fn write_tiff(path: &str, img: &StitchedImage, bits16: bool) -> Result<(), String> {
+    use tiff::encoder::colortype::{RGBA16, RGBA8};
+    use tiff::encoder::{Compression, TiffEncoder};
+
+    let n = (img.w as usize) * (img.h as usize);
+    let file = std::fs::File::create(path).map_err(|e| format!("{}: {e}", path))?;
+    let mut enc = TiffEncoder::new(file)
+        .map_err(|e| format!("{}: {e}", path))?
+        .with_compression(Compression::Deflate(Default::default()));
+    if bits16 {
+        let mut buf: Vec<u16> = Vec::with_capacity(n * 4);
+        for p in &img.pixels {
+            buf.extend_from_slice(&[p[0], p[1], p[2], MAX_VALUE as u16]);
+        }
+        enc.write_image::<RGBA16>(img.w, img.h, &buf)
+            .map_err(|e| format!("{}: {e}", path))?;
+    } else {
+        let mut buf: Vec<u8> = Vec::with_capacity(n * 4);
+        for p in &img.pixels {
+            buf.extend_from_slice(&[(p[0] / 255) as u8, (p[1] / 255) as u8, (p[2] / 255) as u8, 255]);
+        }
+        enc.write_image::<RGBA8>(img.w, img.h, &buf)
+            .map_err(|e| format!("{}: {e}", path))?;
+    }
+    Ok(())
+}
+
+/// Write a JPEG (quality 90, like the C `JPEGOutputImage`).
+pub fn write_jpeg(path: &str, img: &StitchedImage, quality: u8) -> Result<(), String> {
+    use image::codecs::jpeg::JpegEncoder;
+    use image::ImageEncoder;
+
+    let rgb: Vec<u8> = img
+        .pixels
+        .iter()
+        .flat_map(|p| [(p[0] / 255) as u8, (p[1] / 255) as u8, (p[2] / 255) as u8])
+        .collect();
+    let enc = JpegEncoder::new_with_quality(
+        std::io::BufWriter::new(std::fs::File::create(path).map_err(|e| format!("{}: {e}", path))?),
+        quality,
+    );
+    enc.write_image(&rgb, img.w, img.h, image::ExtendedColorType::Rgb8)
+        .map_err(|e| format!("{}: {e}", path))?;
+    Ok(())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn solid(v: u8, w: u32, h: u32) -> RgbImage {
+        RgbImage::from_pixel(w, h, image::Rgb([v, v, v]))
+    }
+
+    #[test]
+    fn nearest_rounds_half_to_even() {
+        let img = solid(100, 4, 4);
+        // 1.5 -> 2 (ties to even), 2.5 -> 2
+        assert_eq!(get_pixel_nearest(&img, 1.5, 0.0), Some([100 * 255; 3]));
+        assert_eq!(get_pixel_nearest(&img, 2.5, 0.0), Some([100 * 255; 3]));
+        assert_eq!(get_pixel_nearest(&img, -0.4, 0.0), Some([100 * 255; 3]));
+        assert_eq!(get_pixel_nearest(&img, -0.6, 0.0), None);
+        assert_eq!(get_pixel_nearest(&img, f64::NAN, 0.0), None);
+    }
+
+    #[test]
+    fn bicubic_constant_and_bounds() {
+        let img = solid(100, 4, 4);
+        // constant image: cubic reproduces the constant exactly
+        assert_eq!(get_pixel_bicubic(&img, 1.3, 1.7), Some([100 * 255; 3]));
+        // any tap outside the 4x4 neighborhood -> unset
+        assert_eq!(get_pixel_bicubic(&img, 0.0, 2.0), None); // fl_x-1 = -1
+        assert_eq!(get_pixel_bicubic(&img, 2.0, 3.5), None); // fl_y+2 = 5
+        assert_eq!(get_pixel_bicubic(&img, 1.0, 1.0), Some([100 * 255; 3]));
+    }
+
+    #[test]
+    fn bicubic_matches_keys_cubic_on_ramp() {
+        // Keys cubic (a=-1/2): on a linear ramp the kernel yields
+        // a + b*(2 dx^3 - 3 dx^2 + 2 dx), not the linear value.
+        let mut img = RgbImage::new(8, 8);
+        for x in 0..8 {
+            for y in 0..8 {
+                img.put_pixel(x, y, image::Rgb([(x as u8) * 30, 0, 0]));
+            }
+        }
+        let got = get_pixel_bicubic(&img, 3.25, 4.0).unwrap();
+        let a = 3.0 * 30.0 * 255.0;
+        let b = 30.0 * 255.0;
+        let dx = 0.25f64;
+        let want = (a + b * (2.0 * dx * dx * dx - 3.0 * dx * dx + 2.0 * dx)).round_ties_even();
+        assert_eq!(got[0], want as i64);
+    }
+
+    #[test]
+    fn bicubic_exact_at_integer_positions() {
+        // at integer sample coordinates (dx = dy = 0) the cubic kernel
+        // reproduces the center tap exactly -> a 1:1 resample is lossless
+        let mut img = RgbImage::new(8, 8);
+        for x in 0..8 {
+            for y in 0..8 {
+                img.put_pixel(x, y, image::Rgb([(x * 7 + y * 3) as u8, (y * 5) as u8, 9]));
+            }
+        }
+        for y in 1..6 {
+            for x in 1..6 {
+                let p = img.get_pixel(x, y).0;
+                assert_eq!(
+                    get_pixel_bicubic(&img, x as f64, y as f64).unwrap(),
+                    [p[0] as i64 * 255, p[1] as i64 * 255, p[2] as i64 * 255]
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn resample_rejects_empty_range() {
+        let img = StitchSource {
+            pan: Panorama::new(),
+            img: solid(50, 8, 8),
+        };
+        let out = resample(&[img], ScanMode::Nearest, 4, 4, 10.0, 10.0);
+        assert_eq!(out.hit.len(), 16);
+        assert!(out.hit.iter().all(|h| !h));
+    }
+}
diff --git a/crates/gipfel-core/src/track.rs b/crates/gipfel-core/src/track.rs
new file mode 100644
index 0000000..5450aa8
--- /dev/null
+++ b/crates/gipfel-core/src/track.rs
@@ -0,0 +1,248 @@
+//! GPS track loading (GPX).
+//!
+//! Minimal streaming parser built directly on quick-xml: we only need the
+//! `trkpt` positions, and real-world files (Garmin etc.) are full of vendor
+//! `<extensions>` blocks that stricter parsers choke on. Waypoints and
+//! routes are ignored; only track points are read.
+//!
+//! Security posture for untrusted files:
+//! - quick-xml is used with `default-features = false` (no serde, no encoding
+//!   detection); it never fetches external DTDs/entities and expands only the
+//!   five predefined XML entities, so XXE and entity-expansion ("billion
+//!   laughs") attacks are structurally impossible; unknown entity references
+//!   make the attribute fail to parse and the point is skipped
+//! - input size is capped (MAX_GPX_BYTES) before anything is parsed
+//! - the parser only reads lat/lon attributes and <ele> text; all other
+//!   elements (including vendor extensions) are ignored without being
+//!   interpreted
+
+use crate::hill::{Hill, FLAG_TRACK_POINT};
+use quick_xml::events::Event;
+
+/// Reject absurdly large track files (a 64 MiB GPX is already ~500k points).
+pub const MAX_GPX_BYTES: u64 = 64 * 1024 * 1024;
+
+/// Load track points from a GPX file (all tracks/segments, in file order).
+///
+/// Consecutive duplicate positions are dropped; the result is decimated with
+/// a uniform stride to at most `max_points` (GPX logs easily reach 10k
+/// points, which would bloat the panorama's hill list).
+pub fn load_gpx(path: &str, max_points: usize) -> Result<Vec<Hill>, String> {
+    let size = std::fs::metadata(path).map_err(|e| format!("{path}: {e}"))?;
+    if size.len() > MAX_GPX_BYTES {
+        return Err(format!("{path}: file too large ({} bytes)", size.len()));
+    }
+    let data = std::fs::read_to_string(path).map_err(|e| format!("{path}: {e}"))?;
+    let mut reader = quick_xml::Reader::from_str(&data);
+    reader.config_mut().trim_text(true);
+
+    let mut pts: Vec<(f64, f64, f64)> = Vec::new(); // (lat, lon, ele) degrees
+    let mut depth: usize = 0;
+    let mut trkpt_depth: Option<usize> = None;
+    let mut in_ele = false;
+    let (mut lat, mut lon, mut ele) = (f64::NAN, f64::NAN, 0.0);
+
+    fn push(pts: &mut Vec<(f64, f64, f64)>, lat: f64, lon: f64, ele: f64) {
+        if lat.is_nan() || lon.is_nan() {
+            return;
+        }
+        if let Some(&(lata, lona, _)) = pts.last() {
+            if lata == lat && lona == lon {
+                return; // consecutive duplicate fix
+            }
+        }
+        pts.push((lat, lon, ele));
+    }
+
+    loop {
+        match reader.read_event().map_err(|e| format!("{path}: {e}"))? {
+            Event::Eof => break,
+            Event::Empty(e) => {
+                if e.local_name().as_ref() == b"trkpt" {
+                    let (la, lo) = read_latlon(&e);
+                    lat = la;
+                    lon = lo;
+                    ele = 0.0;
+                    push(&mut pts, lat, lon, ele);
+                }
+            }
+            Event::Start(e) => {
+                depth += 1;
+                match e.local_name().as_ref() {
+                    b"trkpt" => {
+                        trkpt_depth = Some(depth);
+                        ele = 0.0;
+                        let (la, lo) = read_latlon(&e);
+                        lat = la;
+                        lon = lo;
+                    }
+                    // accept <ele> only as a direct child of <trkpt>, so
+                    // elevation-like tags inside vendor extensions are ignored
+                    b"ele" => {
+                        in_ele = trkpt_depth == Some(depth - 1);
+                        if in_ele {
+                            ele = 0.0;
+                        }
+                    }
+                    _ => {}
+                }
+            }
+            Event::Text(t) => {
+                if in_ele {
+                    if let Ok(s) = t.decode() {
+                        if let Ok(v) = s.trim().parse() {
+                            ele = v;
+                        }
+                    }
+                }
+            }
+            Event::End(e) => {
+                match e.local_name().as_ref() {
+                    b"ele" => in_ele = false,
+                    b"trkpt" => {
+                        push(&mut pts, lat, lon, ele);
+                        trkpt_depth = None;
+                    }
+                    _ => {}
+                }
+                depth = depth.saturating_sub(1);
+            }
+            _ => {}
+        }
+    }
+
+    if pts.is_empty() {
+        return Err(format!("{path}: no track points"));
+    }
+
+    let stride = pts.len().div_ceil(max_points.max(1));
+    Ok(pts
+        .iter()
+        .step_by(stride)
+        .map(|&(lat, lon, ele)| {
+            let mut h = Hill::new("track", lat.to_radians(), lon.to_radians(), ele);
+            h.flags |= FLAG_TRACK_POINT;
+            h
+        })
+        .collect())
+}
+
+fn read_latlon(e: &quick_xml::events::BytesStart) -> (f64, f64) {
+    let (mut lat, mut lon) = (f64::NAN, f64::NAN);
+    for a in e.attributes().flatten() {
+        if let Ok(v) = a.normalized_value(quick_xml::XmlVersion::Implicit1_0) {
+            match a.key.local_name().as_ref() {
+                b"lat" => lat = v.trim().parse().unwrap_or(f64::NAN),
+                b"lon" => lon = v.trim().parse().unwrap_or(f64::NAN),
+                _ => {}
+            }
+        }
+    }
+    (lat, lon)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn tmp(name: &str, content: &str) -> std::path::PathBuf {
+        let dir = std::env::temp_dir().join(format!("gipfel-gpx-test-{name}"));
+        std::fs::create_dir_all(&dir).unwrap();
+        let p = dir.join("t.gpx");
+        std::fs::write(&p, content).unwrap();
+        p
+    }
+
+    #[test]
+    fn gpx_track_roundtrip() {
+        let gpx = r#"<?xml version="1.0" encoding="UTF-8"?>
+<gpx version="1.1" creator="test" xmlns="http://www.topografix.com/GPX/1/1">
+  <trk><name>tour</name><trkseg>
+    <trkpt lat="47.30" lon="11.10"><ele>1000</ele></trkpt>
+    <trkpt lat="47.31" lon="11.11"><ele>1100</ele></trkpt>
+    <trkpt lat="47.31" lon="11.11"><ele>1100</ele></trkpt>
+    <trkpt lat="47.32" lon="11.12"><ele>1200</ele></trkpt>
+  </trkseg></trk>
+</gpx>"#;
+        let p = tmp("rt", gpx);
+        let hills = load_gpx(p.to_str().unwrap(), 100).unwrap();
+        assert_eq!(hills.len(), 3); // duplicate dropped
+        assert!(hills.iter().all(|h| h.is_track_point()));
+        assert!((hills[1].phi - 47.31f64.to_radians()).abs() < 1e-12);
+        assert!((hills[2].height - 1200.0).abs() < 1e-9);
+    }
+
+    #[test]
+    fn gpx_vendor_extensions_ignored() {
+        // Garmin-style files are full of vendor extension blocks; they must
+        // not confuse the parser (and must not contribute elevation).
+        let gpx = r#"<?xml version="1.0" encoding="utf-8"?>
+<gpx creator="Garmin Desktop App" version="1.1" xmlns="http://www.topografix.com/GPX/1/1"
+     xmlns:gpxtpx="http://www.garmin.com/xmlschemas/TrackPointExtension/v1">
+  <metadata><extensions><gpxx:TrackPointExtension xmlns:gpxx="x" version="1.0"/></extensions></metadata>
+  <trk><name>T</name><extensions><gpxx2:x xmlns:gpxx2="y"><ele>9999</ele></gpxx2:x></extensions>
+    <trkseg>
+      <trkpt lat="47.30" lon="11.10">
+        <ele>1000</ele>
+        <extensions><gpxtpx:TrackPointExtension><gpxtpx:hr>120</gpxtpx:hr></gpxtpx:TrackPointExtension></extensions>
+      </trkpt>
+      <trkpt lat="47.31" lon="11.11"><ele>1100</ele></trkpt>
+    </trkseg>
+  </trk>
+</gpx>"#;
+        let p = tmp("ext", gpx);
+        let hills = load_gpx(p.to_str().unwrap(), 100).unwrap();
+        assert_eq!(hills.len(), 2);
+        assert!((hills[0].height - 1000.0).abs() < 1e-9); // not 9999
+    }
+
+    #[test]
+    fn gpx_decimation() {
+        let mut pts = String::new();
+        for i in 0..1000 {
+            pts.push_str(&format!(
+                "<trkpt lat=\"{}\" lon=\"{}\"><ele>1</ele></trkpt>",
+                47.0 + i as f64 * 1e-5,
+                11.0 + i as f64 * 1e-5
+            ));
+        }
+        let gpx = format!(
+            r#"<?xml version="1.0"?><gpx version="1.1" creator="t" xmlns="http://www.topografix.com/GPX/1/1"><trk><trkseg>{pts}</trkseg></trk></gpx>"#
+        );
+        let p = tmp("dec", &gpx);
+        let hills = load_gpx(p.to_str().unwrap(), 100).unwrap();
+        assert!(hills.len() <= 100);
+    }
+
+    #[test]
+    fn gpx_bom_is_tolerated() {
+        // Garmin Desktop App writes a UTF-8 BOM; it must not break parsing.
+        let gpx = format!(
+            "\u{feff}<?xml version=\"1.0\"?><gpx version=\"1.1\" creator=\"t\" xmlns=\"http://www.topografix.com/GPX/1/1\"><trk><trkseg><trkpt lat=\"47.30\" lon=\"11.10\"/></trkseg></trk></gpx>"
+        );
+        let p = tmp("bom", &gpx);
+        let hills = load_gpx(p.to_str().unwrap(), 100).unwrap();
+        assert_eq!(hills.len(), 1);
+    }
+
+    #[test]
+    fn gpx_entity_injection_is_inert() {
+        // XXE / entity-expansion attempts: quick-xml only expands the five
+        // predefined entities; unknown references make the attribute fail to
+        // parse, so the point is simply skipped (never expanded, never fetched).
+        let gpx = r#"<?xml version="1.0"?>
+<!DOCTYPE gpx [ <!ENTITY evil "47.30"> <!ENTITY sys SYSTEM "file:///etc/passwd"> ]>
+<gpx version="1.1" creator="t" xmlns="http://www.topografix.com/GPX/1/1">
+  <trk><trkseg>
+    <trkpt lat="&evil;" lon="11.10"><ele>1000</ele></trkpt>
+    <trkpt lat="47.31" lon="&sys;"><ele>1100</ele></trkpt>
+    <trkpt lat="47.32" lon="11.12"><ele>1200</ele></trkpt>
+  </trkseg></trk>
+</gpx>"#;
+        let p = tmp("ent", gpx);
+        let hills = load_gpx(p.to_str().unwrap(), 100).unwrap();
+        // only the clean point survives; the entity-bearing ones are skipped
+        assert_eq!(hills.len(), 1);
+        assert!((hills[0].phi - 47.32f64.to_radians()).abs() < 1e-12);
+    }
+}
diff --git a/crates/gipfel-core/tests/golden.rs b/crates/gipfel-core/tests/golden.rs
new file mode 100644
index 0000000..fb12e40
--- /dev/null
+++ b/crates/gipfel-core/tests/golden.rs
@@ -0,0 +1,511 @@
+//! Golden-value tests against the original C++ implementation.
+//!
+//! Goldens are generated by `tests/gen_golden.sh` (needs the built C++ tree).
+//! Run with `cargo test --test golden` (unit tests: `cargo test --lib`).
+
+use gipfel_core::hill;
+use gipfel_core::metadata::ImageMetaData;
+use gipfel_core::panorama::Panorama;
+use gipfel_core::projection::{
+    cylindrical::ProjectionCylindrical, rectilinear::ProjectionRectilinear, ProjArgs, Projection, ProjectionType,
+};
+
+fn golden(name: &str) -> Vec<String> {
+    let path = format!(
+        "{}/../../tests/golden/{}",
+        env!("CARGO_MANIFEST_DIR"),
+        name
+    );
+    std::fs::read_to_string(&path)
+        .unwrap_or_else(|_| panic!("missing golden file {} (run tests/gen_golden.sh)", path))
+        .lines()
+        .map(|l| l.to_string())
+        .collect()
+}
+
+fn data_path() -> Option<String> {
+    if let Ok(p) = std::env::var("GIPFEL_DATA") {
+        if std::path::Path::new(&p).exists() {
+            return Some(p);
+        }
+    }
+    let home = std::env::var("HOME").unwrap_or_default();
+    let p = format!("{}/gipfel/gipfel.dat", home);
+    std::path::Path::new(&p).exists().then_some(p)
+}
+
+/// relative-or-absolute closeness (libm differences between glibc and Rust)
+fn close(a: f64, b: f64) -> bool {
+    if a.is_nan() && b.is_nan() {
+        return true;
+    }
+    if a.is_infinite() || b.is_infinite() {
+        return a == b;
+    }
+    (a - b).abs() <= 1e-9f64.max(a.abs() * 1e-9)
+}
+
+#[test]
+fn generated_functions_match_c() {
+    let golden = golden("funcs.txt");
+    let rect = ProjectionRectilinear;
+    let cyl = ProjectionCylindrical;
+
+    assert_eq!(golden.len(), 200);
+
+    for (k, line) in golden.iter().enumerate() {
+        let kk = k as f64;
+        let c_view = -3.0 + 0.03 * kk;
+        let c_pitch = 0.10 * (0.11 * kk).sin();
+        let c_roll = 0.05 * (0.07 * kk).cos();
+        let scale = 2000.0 + 37.0 * (kk % 19.0);
+        let k0 = 0.001 * (0.31 * kk).sin();
+        let k1 = 0.0001 * (0.17 * kk).cos();
+        let x0 = 50.0 * (0.13 * kk).sin();
+        let m_view = c_view + 0.5 * (0.23 * kk).sin();
+        let m_elev = 0.05 * (0.19 * kk).cos();
+
+        let a = ProjArgs { c_view, c_pitch, c_roll, scale, k0, k1, x0, m_view, m_elev };
+
+        let mine: Vec<f64> = vec![
+            rect.mac_x(a),
+            rect.mac_y(a),
+            rect.mac_x_dc_view(a),
+            rect.mac_x_dc_pitch(a),
+            rect.mac_x_dc_roll(a),
+            rect.mac_x_dscale(a),
+            rect.mac_x_dk0(a),
+            rect.mac_x_dk1(a),
+            rect.mac_x_dx0(a),
+            rect.mac_y_dc_view(a),
+            rect.mac_y_dc_pitch(a),
+            rect.mac_y_dc_roll(a),
+            rect.mac_y_dscale(a),
+            rect.mac_y_dk0(a),
+            rect.mac_y_dk1(a),
+            rect.mac_y_dx0(a),
+            cyl.mac_x(a),
+            cyl.mac_y(a),
+            cyl.mac_x_dc_view(a),
+            cyl.mac_x_dc_pitch(a),
+            cyl.mac_x_dc_roll(a),
+            cyl.mac_x_dscale(a),
+            cyl.mac_y_dc_view(a),
+            cyl.mac_y_dc_pitch(a),
+            cyl.mac_y_dc_roll(a),
+            cyl.mac_y_dscale(a),
+        ];
+
+        let theirs: Vec<f64> = line.split_whitespace().map(|t| t.parse().unwrap()).collect();
+        assert_eq!(mine.len(), theirs.len(), "line {}", k);
+
+        for (i, (m, t)) in mine.iter().zip(theirs.iter()).enumerate() {
+            assert!(close(*m, *t), "k={} fn#{}: rust {} != c {}", k, i, m, t);
+        }
+    }
+}
+
+/// derivatives of the generated code vs. central finite differences
+#[test]
+fn derivatives_match_finite_differences() {
+    let rect = ProjectionRectilinear;
+    let base = ProjArgs {
+        c_view: 0.3,
+        c_pitch: 0.05,
+        c_roll: 0.02,
+        scale: 3000.0,
+        k0: 0.001,
+        k1: 0.0002,
+        x0: 12.0,
+        m_view: 0.5,
+        m_elev: 0.04,
+    };
+
+    let h = 1e-6f64;
+
+    // (arg index, analytic derivative fn)
+    let dx: [(usize, &dyn Fn(ProjArgs) -> f64); 7] = [
+        (0, &|a| rect.mac_x_dc_view(a)),
+        (1, &|a| rect.mac_x_dc_pitch(a)),
+        (2, &|a| rect.mac_x_dc_roll(a)),
+        (3, &|a| rect.mac_x_dscale(a)),
+        (4, &|a| rect.mac_x_dk0(a)),
+        (5, &|a| rect.mac_x_dk1(a)),
+        (6, &|a| rect.mac_x_dx0(a)),
+    ];
+    let dy: [(usize, &dyn Fn(ProjArgs) -> f64); 7] = [
+        (0, &|a| rect.mac_y_dc_view(a)),
+        (1, &|a| rect.mac_y_dc_pitch(a)),
+        (2, &|a| rect.mac_y_dc_roll(a)),
+        (3, &|a| rect.mac_y_dscale(a)),
+        (4, &|a| rect.mac_y_dk0(a)),
+        (5, &|a| rect.mac_y_dk1(a)),
+        (6, &|a| rect.mac_y_dx0(a)),
+    ];
+
+    let set = |a: ProjArgs, idx: usize, v: f64| -> ProjArgs {
+        let mut b = a;
+        match idx {
+            0 => b.c_view = v,
+            1 => b.c_pitch = v,
+            2 => b.c_roll = v,
+            3 => b.scale = v,
+            4 => b.k0 = v,
+            5 => b.k1 = v,
+            _ => b.x0 = v,
+        }
+        b
+    };
+
+    for (idx, f) in dx {
+        let plus = rect.mac_x(set(base, idx, base_val(&base, idx) + h));
+        let minus = rect.mac_x(set(base, idx, base_val(&base, idx) - h));
+        let fd = (plus - minus) / (2.0 * h);
+        let an = f(base);
+        assert!(
+            (an - fd).abs() <= 1e-4f64.max(an.abs() * 1e-5),
+            "dx deriv arg {}: analytic {} finite-diff {}",
+            idx,
+            an,
+            fd
+        );
+    }
+    for (idx, f) in dy {
+        let plus = rect.mac_y(set(base, idx, base_val(&base, idx) + h));
+        let minus = rect.mac_y(set(base, idx, base_val(&base, idx) - h));
+        let fd = (plus - minus) / (2.0 * h);
+        let an = f(base);
+        assert!(
+            (an - fd).abs() <= 1e-4f64.max(an.abs() * 1e-5),
+            "dy deriv arg {}: analytic {} finite-diff {}",
+            idx,
+            an,
+            fd
+        );
+    }
+}
+
+fn base_val(a: &ProjArgs, idx: usize) -> f64 {
+    match idx {
+        0 => a.c_view,
+        1 => a.c_pitch,
+        2 => a.c_roll,
+        3 => a.scale,
+        4 => a.k0,
+        5 => a.k1,
+        _ => a.x0,
+    }
+}
+
+/// full fit (geodesy + projection + GSL LM) against the C++ fitdump
+fn fit_golden(projtype: ProjectionType, file: &str) {
+    let Some(data) = data_path() else {
+        eprintln!("gipfel.dat not found, skipping");
+        return;
+    };
+    let golden = golden(file);
+
+    let mut pan = Panorama::new();
+    pan.load_data(&data).unwrap();
+    pan.set_projection(projtype);
+
+    pan.set_view_lat(47.36);
+    pan.set_view_long(11.15);
+    pan.set_view_height(2100.0);
+
+    // ground truth
+    pan.set_center_angle(20.0);
+    pan.set_pitch_angle(0.0);
+    pan.set_roll_angle(1.0);
+    pan.set_scale(3000.0);
+
+    // snapshot first 6 visible hills as observations
+    let mut obs: Vec<hill::Hill> = pan.visible[..6].iter().map(|&i| pan.hills[i].clone()).collect();
+    let n = obs.len();
+    let known: Vec<usize> = (0..n).collect();
+    // the C harness copies a_elev explicitly (C copy ctor omits it)
+    for (k, &i) in pan.visible[..6].iter().enumerate() {
+        obs[k].a_elev = pan.hills[i].a_elev;
+    }
+
+    // perturb initial guess
+    pan.set_center_angle(pan.get_center_angle() + 0.5);
+    pan.set_pitch_angle(0.3);
+    pan.set_roll_angle(0.5);
+    pan.set_scale(2800.0);
+
+    let ret = pan.comp_params_with(&obs, &known, true);
+
+    let g_ret: i32 = golden[0]["ret ".len()..].parse().unwrap();
+    assert_eq!(ret.is_ok(), g_ret == 0);
+
+    let g: Vec<f64> = golden[1..8]
+        .iter()
+        .map(|l| l.split_whitespace().nth(1).unwrap().parse().unwrap())
+        .collect();
+
+    let mine = [
+        pan.parms.a_center,
+        pan.parms.a_pitch,
+        pan.parms.a_roll,
+        pan.parms.scale,
+        pan.parms.k0,
+        pan.parms.k1,
+        pan.parms.x0,
+    ];
+    for (m, t) in mine.iter().zip(g.iter()) {
+        // solver results: allow small drift, exact zeros must stay ~zero
+        assert!(
+            (m - t).abs() <= 1e-6f64.max(m.abs() * 1e-6),
+            "fit param {} vs {}",
+            m,
+            t
+        );
+    }
+
+    // final coordinates
+    let proj = gipfel_core::projection::projection_from_type(projtype);
+    for line in &golden[8..] {
+        let parts: Vec<&str> = line.split_whitespace().collect();
+        if parts.len() != 4 {
+            continue;
+        }
+        let i: usize = parts[1].parse().unwrap();
+        let tx: f64 = parts[2].parse().unwrap();
+        let ty: f64 = parts[3].parse().unwrap();
+        let (x, y) = proj.get_coordinates(obs[i].alph, obs[i].a_elev, &pan.parms);
+        assert!(close(x, tx), "hill {} x: {} vs {}", i, x, tx);
+        assert!(close(y, ty), "hill {} y: {} vs {}", i, y, ty);
+    }
+}
+
+#[test]
+fn fit_rectilinear_matches_c() {
+    fit_golden(ProjectionType::Rectilinear, "fit_rect.txt");
+}
+
+#[test]
+fn fit_cylindrical_matches_c() {
+    fit_golden(ProjectionType::Cylindrical, "fit_cyl.txt");
+}
+
+#[test]
+fn real_database_loads() {
+    let Some(data) = data_path() else {
+        eprintln!("gipfel.dat not found, skipping");
+        return;
+    };
+    let hills = hill::load(&data).unwrap();
+    assert!(hills.len() > 21000, "loaded {}", hills.len());
+    assert!(hills.iter().any(|h| h.name == "ZUGSPITZE"));
+}
+
+#[test]
+fn viewpoint_geometry() {
+    let Some(data) = data_path() else {
+        eprintln!("gipfel.dat not found, skipping");
+        return;
+    };
+    let mut pan = Panorama::new();
+    pan.load_data(&data).unwrap();
+    let idx = pan.get_pos("ZUGSPITZE").expect("ZUGSPITZE in db");
+    pan.set_viewpoint(Some(idx));
+
+    // NOTE: set_viewpoint re-sorts hills by azimuth; look the name up again
+    let idx = pan.get_pos("ZUGSPITZE").unwrap();
+    let m = &pan.hills[idx];
+    assert!((m.height - 2963.0).abs() < 1.0);
+    // distance to self is zero
+    assert!(pan.distance(m.phi, m.lam).abs() < 1e-9);
+    // a far hill: Mont Blanc (~45.83 N, 6.86 E)
+    let mb = pan.get_pos("MONT BLANC").expect("MONT BLANC in db");
+    let d = pan.distance(pan.hills[mb].phi, pan.hills[mb].lam);
+    let km = d * 6371000.785 / 1000.0;
+    assert!(km > 340.0 && km < 380.0, "Mont Blanc distance {} km", km);
+}
+
+#[test]
+fn export_hills_matches_c() {
+    let Some(data) = data_path() else {
+        eprintln!("gipfel.dat not found, skipping");
+        return;
+    };
+    let jpg = format!(
+        "{}/../../tests/golden/test.jpg",
+        env!("CARGO_MANIFEST_DIR")
+    );
+    let golden = golden("export.txt");
+
+    let (w, h) = gipfel_core::metadata::image_dimensions(&jpg).expect("test.jpg readable");
+    let mut md = ImageMetaData::new();
+    md.load_image(&jpg);
+
+    let mut pan = Panorama::new();
+    assert!(gipfel_core::metadata::apply_image_metadata(&mut pan, &md, w, h));
+    pan.load_data(&data).unwrap();
+    pan.set_height_dist_ratio(0.07);
+
+    let mut out = vec![
+        "#".to_string(),
+        "# name\theight\tx\ty\tdistance\tflags".to_string(),
+        "#".to_string(),
+    ];
+    out.extend(pan.export_hill_rows(w as i64, h as i64, false));
+
+    assert_eq!(out.len(), golden.len());
+    for (a, b) in out.iter().zip(golden.iter()) {
+        assert_eq!(a, b);
+    }
+}
+
+#[test]
+fn metadata_comment_roundtrip_and_c_format() {
+    let mut md = ImageMetaData::new();
+    md.longitude = 11.1533;
+    md.latitude = 47.3601;
+    md.height = 2100.5;
+    md.direction = 20.0;
+    md.pitch = 0.3;
+    md.roll = 1.0;
+    md.focal_length_35mm = 35.0;
+    md.projection_type = 1;
+    md.k0 = 0.001;
+    md.k1 = -0.0002;
+    md.x0 = 12.0;
+
+    let s = md.format_comment();
+    // exact C printf("%f") compatibility
+    assert_eq!(
+        s,
+        "gipfel: longitude 11.153300, latitude 47.360100, height 2100.500000, direction 20.000000, nick 0.300000, tilt 1.000000, focal_length_35mm 35.000000, projection type 1, k0 0.001000, k1 -0.000200, x0 12.000000"
+    );
+
+    let mut md2 = ImageMetaData::new();
+    assert!(md2.parse_gipfel_comment(s.as_bytes()));
+    assert!((md2.longitude - 11.1533).abs() < 1e-9);
+    assert!((md2.latitude - 47.3601).abs() < 1e-9);
+    assert!((md2.height - 2100.5).abs() < 1e-9);
+    assert!((md2.direction - 20.0).abs() < 1e-9);
+    assert_eq!(md2.projection_type, 1);
+    assert!((md2.k0 - 0.001).abs() < 1e-9);
+    assert!((md2.k1 + 0.0002).abs() < 1e-9);
+    assert!((md2.x0 - 12.0).abs() < 1e-9);
+
+    // short form (8 fields, pre-distortion files): k0/k1/x0 stay NaN
+    let short = b"gipfel: longitude 1.000000, latitude 2.000000, height 3.000000, direction 4.000000, nick 5.000000, tilt 6.000000, focal_length_35mm 7.000000, projection type 0";
+    let mut md3 = ImageMetaData::new();
+    assert!(md3.parse_gipfel_comment(short));
+    assert!(md3.k0.is_nan());
+    assert_eq!(md3.projection_type, 0);
+
+    // trailing garbage after x0 must not break parsing (sscanf semantics)
+    let mut md4 = ImageMetaData::new();
+    let mut s4 = s.clone();
+    s4.push_str("  some other comment junk");
+    assert!(md4.parse_gipfel_comment(s4.as_bytes()));
+    assert!((md4.x0 - 12.0).abs() < 1e-9);
+}
+
+/// The least-squares fit must run on a small stack: the LM solver keeps its
+/// working buffers (residuals, Jacobian, normal equations) on the heap, so a
+/// 128 KiB thread stack is plenty. Guards against regressions that move those
+/// buffers back onto the stack.
+#[test]
+fn fit_survives_small_stack() {
+    let Some(data) = data_path() else {
+        eprintln!("gipfel.dat not found, skipping");
+        return;
+    };
+    std::thread::Builder::new()
+        .stack_size(128 * 1024)
+        .spawn(move || {
+            let mut pan = Panorama::new();
+            pan.load_data(&data).unwrap();
+            pan.set_view_lat(47.36);
+            pan.set_view_long(11.15);
+            pan.set_view_height(2100.0);
+            pan.set_center_angle(20.0);
+            pan.set_pitch_angle(0.0);
+            pan.set_roll_angle(1.0);
+            pan.set_scale(3000.0);
+            let known: Vec<usize> = pan.visible[..4].to_vec();
+            // perturb so the solver actually iterates
+            pan.set_center_angle(20.5);
+            pan.set_scale(2800.0);
+            pan.comp_params(&known, true).unwrap();
+        })
+        .unwrap()
+        .join()
+        .unwrap();
+}
+
+/// Stitch resample against the C Stitch::resample + ScanImage pipeline.
+fn stitch_golden(
+    name: &str,
+    proj: ProjectionType,
+    center: f64,
+    pitch: f64,
+    roll: f64,
+    scale: f64,
+    k0: f64,
+    k1: f64,
+    bicubic: bool,
+) {
+    use gipfel_core::stitch::{self, ScanMode};
+
+    let ppm = format!(
+        "{}/../../tests/golden/stitch.ppm",
+        env!("CARGO_MANIFEST_DIR")
+    );
+    let img = image::open(&ppm)
+        .unwrap_or_else(|e| panic!("missing {}: {e}", ppm))
+        .to_rgb8();
+
+    let mut pan = Panorama::new();
+    pan.set_projection(proj);
+    pan.set_center_angle(center);
+    pan.set_pitch_angle(pitch);
+    pan.set_roll_angle(roll);
+    pan.set_scale(scale);
+    pan.set_distortion_params(k0, k1, 0.0);
+
+    let src = stitch::StitchSource { pan, img };
+    let mode = if bicubic { ScanMode::Bicubic } else { ScanMode::Nearest };
+    let out = stitch::resample(&[src], mode, 64, 48, 10.0, 30.0);
+
+    let mut checked = 0;
+    for line in golden(name) {
+        let f: Vec<&str> = line.split_whitespace().collect();
+        let (x, y, hit): (usize, usize, i32) = (f[0].parse().unwrap(), f[1].parse().unwrap(), f[2].parse().unwrap());
+        let want = [
+            f[3].parse::<u16>().unwrap(),
+            f[4].parse::<u16>().unwrap(),
+            f[5].parse::<u16>().unwrap(),
+        ];
+        let i = y * 64 + x;
+        assert_eq!(out.hit[i], hit == 1, "{name} pixel ({x},{y}) hit flag");
+        assert_eq!(out.pixels[i], want, "{name} pixel ({x},{y})");
+        checked += 1;
+    }
+    assert!(checked > 1000);
+}
+
+#[test]
+fn stitch_rect_nearest_matches_c() {
+    stitch_golden("stitch_rect_near.txt", ProjectionType::Rectilinear, 20.0, 2.0, 0.0, 200.0, 0.0, 0.0, false);
+}
+
+#[test]
+fn stitch_rect_bicubic_matches_c() {
+    stitch_golden("stitch_rect_bic.txt", ProjectionType::Rectilinear, 20.0, 2.0, 0.0, 200.0, 0.0, 0.0, true);
+}
+
+#[test]
+fn stitch_cyl_bicubic_matches_c() {
+    stitch_golden("stitch_cyl_bic.txt", ProjectionType::Cylindrical, 20.0, 2.0, 0.0, 200.0, 0.0, 0.0, true);
+}
+
+#[test]
+fn stitch_distorted_bicubic_matches_c() {
+    stitch_golden("stitch_dist_bic.txt", ProjectionType::Rectilinear, 20.0, 2.0, 3.0, 200.0, -0.3, 0.1, true);
+}