Harpoons, landing points, step sizes

Documentation/HarpoonAnchoring.md

@@ -0,0 +1,51 @@
+# Harpoons and terrain tethers
+
+`HarpoonAnchoring.swift` extends the compliant landing mechanics with explicit
+terrain anchoring and mass-conserving projectile launch. This is an engineering
+model; supply measured material properties only when they are actually available.
+
+`MechanicalBody.launching` separates a small rigid component whose centre is
+specified in the original body's axes. It subtracts component inertia and the
+parallel-axis contributions, moves the retained centre of mass, and carries
+rotational point velocity into both new bodies. Equal/opposite firing impulses
+supply separation energy and produce recoil/mount torque. No mass is created or
+counted twice. Component inertia is in the same body axes. Invalid resulting
+inertia throws. Use `retainedCentreBodyMetres` to adjust every retained contact
+and attachment coordinate. The projectile starts with the parent's orientation
+and angular velocity; an impulse through its own centre adds no projectile spin.
+
+`HarpoonAnchor` declares the surface `Frame` identity. Firing is one-shot and
+can explicitly misfire. Incoming impacts test incidence and normal kinetic
+energy, using a constant penetration-resistance force:
+
+    estimated depth = m * normalSpeed² / (2 * resistance)
+
+Attachment requires depth between the declared minimum and maximum supported
+burial depths and an adequate incidence cosine. Too-hard, too-shallow, grazing
+and overpenetrating impacts are rejected. This is not a soil mechanics, ice
+fracture or fabric/rock interaction solver. Holding strength is a separate input;
+it is not inferred from a visual texture, illumination or the penetration depth.
+`surfacePointMetres` records the cable exit point on the terrain in surface axes;
+penetration depth is recorded separately, not added as rendered surface relief.
+The application supplies rejected-projectile rebound/contact behavior.
+
+`TerrainTether` evaluates the existing cable tension law against a prescribed
+moving endpoint. Transform the surface-local cable exit point and its velocity
+into the inertial integration frame at every stage. Apply the returned force at
+the lander's body-local mount and record the equal/opposite reaction on the
+prescribed celestial body. Do not integrate a fictitious anchor mass, aggregate
+its inertia with the lander, or reparent the lander's physical motion to it.
+The full celestial ephemeris is not perturbed by these ground reactions.
+
+`commitLoad` checks breaking and pullout at accepted substep boundaries. Cable
+breakage takes precedence if both thresholds are exceeded. A broken or cut cable
+leaves a successfully embedded anchor in place. Pullout instead releases the
+projectile to independent dynamics while its cable can remain connected. Cutting
+or changing an attachment state applies no hidden impulse or velocity reset to
+the lander. The application transfers projectile stopping momentum into the
+prescribed terrain, rather than removing that momentum from the lander twice.
+
+Winch payout/retraction is commanded mechanical work. Motor energy/force limits,
+capstan friction, subsurface stress distribution, temperature-dependent strength,
+foot screws and hold-down thrusters are outside this initial model. Event timing
+is resolved to the chosen mechanical substep, not continuously certified.

Documentation/LandingMechanics.md

@@ -0,0 +1,56 @@
+# Compliant landing mechanics
+
+`CelMek/LandingMechanics.swift` provides independent rigid bodies, cable edges,
+spherical contact lobes/feet, local terrain planes and adiabatic gas cushions.
+This is a compliant engineering model, not a rigid constraint solver or deformable
+fabric model. All quantities use SI units.
+
+A `MechanicalSystem` declares one pseudo-inertial `Frame` and an `Epoch`. Body
+positions and velocities are in that reference, orientations map body axes into
+it, angular velocities and inertia are in body axes, and mount coordinates are
+relative to each body's centre of mass. The connection graph is separate from
+frame parentage. Multiple cables between the same pair of bodies are allowed.
+Rigid assemblies can still use `MassComponent`; detach rigid components before
+creating independent bodies and do not count their masses twice.
+
+Cable extension and attachment-point relative speed produce tension. Slack cables
+exert no force, even when the endpoints are moving apart. Tension is clamped to
+zero rather than allowing compression. Payout subtracts from extension rate and
+advances rest length; retrieval can do positive mechanical work. Each load acts
+at both mounts with equal and opposite reference forces and body-axis torques.
+Cutting changes connectivity only. Breaking commits at accepted substep boundaries;
+`cableLoads` is pure and never changes connectivity. A broken cable contributes no
+load after the break boundary. Event time uncertainty is up to one substep.
+
+Contacts use a plane supplied by the terrain provider, normal spring support,
+compression-only normal damping and speed-regularized Coulomb friction. Forces
+act at the lobe's surface point and include angular point velocity. There is no
+adhesion. Optional rolling resistance applies a dissipative torque proportional
+to normal load and a resistance length, opposing angular motion tangent to the
+plane. Its default is zero; a rigid sphere can roll without friction stopping it.
+Regularized friction approaches zero with slip speed; this is not an
+exact static-friction constraint. Supply local surface velocity explicitly;
+for a rotating surface, evaluate a sufficiently local patch per contact. No
+terrain queries depend on visual LOD.
+
+Gas pressure follows `P V^gamma = gasInventory`, reduced by ambient pressure.
+Pressure support ramps with penetration over half a lobe radius, approximating
+the growing contact footprint rather than applying full-area preload at first
+touch. This footprint approximation is distinct from the gas compression volume.
+Volume has a lower bound; normal contact stiffness supplies the bottom-out guard.
+`vent(retaining:)` removes gas inventory monotonically. Creating a cushion at its
+initial absolute pressure is an explicit inflation initialization, not a modeled
+inflator or gas flow network. Geometry changes, gas vent rates and deployment
+controllers belong to the application. Fabric stress, leaks, rupture, temperature
+transport and inter-lobe plumbing are not included.
+
+`step` uses explicit midpoint integration and normalizes quaternions. External
+loads are evaluated at the beginning and midpoint and must be pure functions of
+sample epoch/state. Gravity, aerodynamic loads and controls belong there; do not
+also apply the internal cable reactions. A failed step leaves the system unchanged.
+Steps must be positive and at most 0.01 s, but that cap is **not** a stability
+guarantee for arbitrary stiffness. Choose a smaller step relative to the highest
+contact/cable natural frequency and check convergence. Open Orbit uses 0.002 s
+for its declared low-mass landing demonstrations. No continuous collision search,
+inter-body collision, rigid joints, line/cable terrain collision or long-duration
+tether orbital integrator is provided.

Sources/CelMek/FiniteBurnStateTargeting.swift

@@ -4,7 +4,9 @@ import Math
 
 /// Two fixed-ignition inertial burns corrected against a full terminal state.
 /// Six duration-vector components are solved with finite differences and
-/// pivoted elimination. Clearance and operational constraints remain caller checks.
+/// pivoted elimination. `relativeControlStep` scales the duration-vector finite
+/// differences (with a 1 ms floor); a smaller value can resolve long-burn
+/// terminal encounters. Existing callers retain the 1e-4 default. Clearance and operational constraints remain caller checks.
 public struct FiniteBurnStateTargetingReport {
     public let best: FiniteTransferEvaluation
     public let evaluations: Int
@@ -16,19 +18,22 @@ public enum FiniteBurnStateTargeting {
         positionToleranceMetres: Double = 10, velocityToleranceMetresPerSecond: Double = 0.01,
         minimumRemainingPropellantKilograms: Double = 0, maximumEvaluations: Int = 160,
         numericalOptions: NumericalPropagationOptions = .preset(for:.high),
-        exhaustVelocityForDirection: ((SIMD3<Double>) throws -> Double)? = nil) throws -> FiniteTransferEvaluation {
+        exhaustVelocityForDirection: ((SIMD3<Double>) throws -> Double)? = nil,
+        relativeControlStep: Double = 1e-4) throws -> FiniteTransferEvaluation {
         try refineWithReport(initialState:initialState,initialSchedule:initialSchedule,target:target,context:context,
             backgroundModels:backgroundModels,positionToleranceMetres:positionToleranceMetres,
             velocityToleranceMetresPerSecond:velocityToleranceMetresPerSecond,minimumRemainingPropellantKilograms:minimumRemainingPropellantKilograms,
-            maximumEvaluations:maximumEvaluations,numericalOptions:numericalOptions,exhaustVelocityForDirection:exhaustVelocityForDirection).best
+            maximumEvaluations:maximumEvaluations,numericalOptions:numericalOptions,exhaustVelocityForDirection:exhaustVelocityForDirection,relativeControlStep:relativeControlStep).best
     }
     public static func refineWithReport(initialState: SpacecraftState, initialSchedule: ScheduledThrust,
         target: CartesianState, context: ForceEvaluationContext, backgroundModels: [any ForceModel],
         positionToleranceMetres: Double = 10, velocityToleranceMetresPerSecond: Double = 0.01,
         minimumRemainingPropellantKilograms: Double = 0, maximumEvaluations: Int = 160,
         numericalOptions: NumericalPropagationOptions = .preset(for:.high),
-        exhaustVelocityForDirection: ((SIMD3<Double>) throws -> Double)? = nil) throws -> FiniteBurnStateTargetingReport {
-        guard initialSchedule.burns.count == 2, (14...10_000).contains(maximumEvaluations),
+        exhaustVelocityForDirection: ((SIMD3<Double>) throws -> Double)? = nil,
+        relativeControlStep: Double = 1e-4) throws -> FiniteBurnStateTargetingReport {
+        guard relativeControlStep.isFinite,relativeControlStep>0,relativeControlStep<=0.01,
+            initialSchedule.burns.count == 2, (14...10_000).contains(maximumEvaluations),
             initialSchedule.frame === initialState.frame,
             positionToleranceMetres.isFinite, velocityToleranceMetresPerSecond.isFinite,
             minimumRemainingPropellantKilograms.isFinite,
@@ -78,7 +83,7 @@ public enum FiniteBurnStateTargeting {
             var matrix = Array(repeating:Array(repeating:0.0,count:7),count:6)
             let r = residual(best)
             for j in 0..<6 {
-                let h = max(0.001,abs(controls[j])*1e-4)
+                let h = max(0.001,abs(controls[j])*relativeControlStep)
                 var a = controls, b = controls; a[j] += h; b[j] -= h
                 let ea = try evaluate(a), eb = try evaluate(b)
                 guard ea != nil || eb != nil else { throw FiniteTransferPlanningError.invalidInput }

Sources/CelMek/HarpoonAnchoring.swift

@@ -0,0 +1,105 @@
+// SPDX-License-Identifier: Apache-2.0
+// Copyright 2026 Mattias Holm
+import Foundation
+import Math
+
+public enum HarpoonPhase: String, Codable { case stowed, flying, anchored, rejected, pulledOut, misfire }
+public enum AnchorCableState: String, Codable { case connected, cut, broken }
+/// Explicit engineering assumptions, not an inference of subsurface strength from
+/// the visible terrain. Normal kinetic energy is spent against constant resistance.
+public struct HarpoonSubstrate {
+    public let penetrationResistanceNewtons, minimumDepthMetres, maximumDepthMetres, holdingForceNewtons, minimumIncidenceCosine: Double
+    public init(penetrationResistanceNewtons: Double, minimumDepthMetres: Double, maximumDepthMetres: Double,
+        holdingForceNewtons: Double, minimumIncidenceCosine: Double = 0.5) throws {
+        guard [penetrationResistanceNewtons,minimumDepthMetres,maximumDepthMetres,holdingForceNewtons].allSatisfy({$0.isFinite && $0 > 0}),
+              maximumDepthMetres >= minimumDepthMetres, minimumIncidenceCosine.isFinite, (0...1).contains(minimumIncidenceCosine) else { throw LandingMechanicsError.invalidParameter }
+        self.penetrationResistanceNewtons=penetrationResistanceNewtons; self.minimumDepthMetres=minimumDepthMetres
+        self.maximumDepthMetres=maximumDepthMetres; self.holdingForceNewtons=holdingForceNewtons; self.minimumIncidenceCosine=minimumIncidenceCosine
+    }
+}
+/// Terrain anchor state is in the declared rotating surface frame. Its tether
+/// couples to an independent inertial body; it never reparents or grounds it.
+public struct HarpoonAnchor {
+    public let id:String
+    public let surfaceFrame:Frame
+    public private(set) var phase:HarpoonPhase = .stowed
+    public private(set) var cableState:AnchorCableState = .connected
+    public private(set) var surfacePointMetres:SIMD3<Double>?
+    public private(set) var penetrationMetres=0.0
+    public private(set) var holdingForceNewtons=0.0
+    public init(id:String,surfaceFrame:Frame) throws {
+        guard !id.isEmpty else { throw LandingMechanicsError.invalidParameter }
+        self.id=id; self.surfaceFrame=surfaceFrame
+    }
+    public mutating func fire(succeeded:Bool = true) throws {
+        guard phase == .stowed else { throw LandingMechanicsError.invalidParameter }
+        phase=succeeded ? .flying : .misfire
+    }
+    /// Returns true only for a captured, incoming impact within the supported
+    /// burial layer. Rejected impacts remain dynamic and can bounce away.
+    @discardableResult public mutating func impact(pointMetres:SIMD3<Double>,in frame:Frame,
+        massKilograms:Double,normalSpeedMetresPerSecond:Double,incidenceCosine:Double,substrate:HarpoonSubstrate) throws -> Bool {
+        guard phase == .flying,frame === surfaceFrame,attitudeFinite(pointMetres),massKilograms.isFinite,massKilograms > 0,
+              normalSpeedMetresPerSecond.isFinite,normalSpeedMetresPerSecond > 0,incidenceCosine.isFinite,(0...1).contains(incidenceCosine) else { throw LandingMechanicsError.invalidParameter }
+        let depth=0.5*massKilograms*normalSpeedMetresPerSecond*normalSpeedMetresPerSecond/substrate.penetrationResistanceNewtons
+        guard depth.isFinite else { throw LandingMechanicsError.numericalRange }
+        penetrationMetres=depth
+        if incidenceCosine < substrate.minimumIncidenceCosine || depth < substrate.minimumDepthMetres || depth > substrate.maximumDepthMetres {
+            phase = .rejected; return false
+        }
+        surfacePointMetres=pointMetres; holdingForceNewtons=substrate.holdingForceNewtons;phase = .anchored;return true
+    }
+    /// Commit only at accepted integration boundaries. Break precedence is
+    /// explicit when both thresholds are crossed. No change to vehicle velocity.
+    public mutating func commitLoad(tensionNewtons:Double,breakingTensionNewtons:Double) throws {
+        guard tensionNewtons.isFinite,tensionNewtons >= 0,breakingTensionNewtons.isFinite,breakingTensionNewtons > 0 else { throw LandingMechanicsError.invalidParameter }
+        guard cableState == .connected else { return }
+        if tensionNewtons > breakingTensionNewtons { cableState = .broken }
+        else if phase == .anchored && tensionNewtons > holdingForceNewtons { phase = .pulledOut }
+    }
+    public mutating func cut() { cableState = .cut }
+}
+
+/// A prescribed terrain endpoint uses the same tension law as two-body cables.
+/// The equal and opposite force is a reaction on the prescribed celestial body,
+/// not a force to apply a second time to the lander or an integrated dummy mass.
+public struct TerrainTether {
+    public private(set) var cable:MechanicalTether
+    public init(id:String,body:String,mountMetres:SIMD3<Double>,lengthMetres:Double,
+        stiffnessNewtonsPerMetre:Double,dampingNewtonSecondsPerMetre:Double,breakingTensionNewtons:Double) throws {
+        cable=try MechanicalTether(id:id,firstBody:body,secondBody:"terrain:\(id)",firstMountMetres:mountMetres,
+            lengthMetres:lengthMetres,stiffnessNewtonsPerMetre:stiffnessNewtonsPerMetre,
+            dampingNewtonSecondsPerMetre:dampingNewtonSecondsPerMetre,breakingTensionNewtons:breakingTensionNewtons)
+    }
+    public mutating func setLength(metres:Double) throws { try cable.setLength(metres:metres) }
+    public func load(body:MechanicalBody,pointReferenceMetres:SIMD3<Double>,velocityReferenceMetresPerSecond:SIMD3<Double>,payoutMetresPerSecond:Double=0) throws -> TetherLoad {
+        let endpoint=try MechanicalBody(id:cable.secondBody,massKilograms:1,inertia:RigidBodyInertia(xx:1,yy:1,zz:1),positionMetres:pointReferenceMetres,velocityMetresPerSecond:velocityReferenceMetresPerSecond)
+        return try cable.load(first:body,second:endpoint,payoutMetresPerSecond:payoutMetresPerSecond)
+    }
+}
+
+public struct MechanicalLaunch {
+    public let retained,projectile:MechanicalBody
+    public let retainedCentreBodyMetres:SIMD3<Double>
+}
+extension MechanicalBody {
+    /// Separate a small rigid component, then apply an internal firing impulse
+    /// at the component centre. Mass, linear and angular momentum are conserved;
+    /// launch energy is supplied by the firing mechanism. Inertia is explicitly
+    /// removed and translated to the remaining body's new centre of mass.
+    public func launching(id projectileID:String,massKilograms m:Double,inertia projectileInertia:RigidBodyInertia,
+        centreBodyMetres r:SIMD3<Double>,impulseBodyNewtonSeconds j:SIMD3<Double>) throws -> MechanicalLaunch {
+        guard m.isFinite,m > 0,m < massKilograms,projectileID != id,attitudeFinite(r),attitudeFinite(j) else { throw LandingMechanicsError.invalidParameter }
+        let remaining=massKilograms-m,s = -r*(m/remaining)
+        func parallel(_ v:SIMD3<Double>,_ mass:Double,_ i:Int,_ k:Int)->Double { mass*((i == k ? (v*v).sum() : 0)-v[i]*v[k]) }
+        func entry(_ i:Int,_ k:Int)->Double { inertia.matrix[i,k]-projectileInertia.matrix[i,k]-parallel(r,m,i,k)-parallel(s,remaining,i,k) }
+        let retainedInertia=try RigidBodyInertia(xx:entry(0,0),yy:entry(1,1),zz:entry(2,2),xy:entry(0,1),xz:entry(0,2),yz:entry(1,2))
+        let impulse=bodyToReference.matrix*j
+        let retained=try MechanicalBody(id:id,massKilograms:remaining,inertia:retainedInertia,positionMetres:point(s),
+            velocityMetresPerSecond:pointVelocity(s)-impulse/remaining,bodyToReference:bodyToReference,
+            angularVelocityBody:angularVelocityBody+retainedInertia.solve(attitudeCross(r-s,-j)))
+        let projectile=try MechanicalBody(id:projectileID,massKilograms:m,inertia:projectileInertia,positionMetres:point(r),
+            velocityMetresPerSecond:pointVelocity(r)+impulse/m,bodyToReference:bodyToReference,angularVelocityBody:angularVelocityBody)
+        return MechanicalLaunch(retained:retained,projectile:projectile,retainedCentreBodyMetres:s)
+    }
+}