Composite rigid bodies and atmosphere improvements

.github/workflows/swift.yml

@@ -20,7 +20,7 @@ jobs:
           - os: macos-latest
             build-system: swiftbuild
           # Use the previous backend while investigating opaque Linux/Windows failures.
-          - os: ubuntu-24.04
+          - os: ubuntu-26.04
             build-system: native
           - os: windows-latest
             build-system: native

Documentation/MassGeometry.md

@@ -0,0 +1,25 @@
+# Rigid mass components
+
+`MassComponent` is a value tree, independent of scene nodes and ephemeris frames.
+Each component defines its own mass, centre and intrinsic inertia, a rigid pose
+relative to its parent, and child components. Structural nodes may have no mass.
+Identifiers are unique throughout an assembly. Intrinsic inertia excludes children.
+Cylinder helpers use the local X axis as the cylinder axis.
+
+`massGeometry()` rotates each intrinsic tensor and applies the parallel-axis
+theorem about the assembled centre. The result is expressed in the root parent's
+axes and origin, including the root's own placement.
+
+`detaching(id:)` returns retained and detached value trees atomically. The detached
+root carries the composed parent transforms, so both results remain in the same
+structural coordinate system. Recompute both geometries, then call
+`offsetKinematics(from:angularVelocity:)` for each result. Rotate these offsets
+into the simulation reference frame and add the original assembly position and
+velocity. Both bodies initially inherit the original orientation and angular
+velocity; their subsequent dynamics use their respective inertia tensors.
+
+Separation springs and motors require explicit impulses with matching recoil.
+This API does not integrate flexible structures, slosh, articulated joints,
+propellant momentum transport, or time-dependent internal motion. Fuel depletion
+can rebuild the value tree at a step boundary; instantaneous reference-point
+kinematics alone is not a complete variable-mass dynamics model.

README.md

@@ -107,3 +107,19 @@ optional incompressible rotor/body CFD and a compressible Sod reference problem.
 Liquid properties and vehicle coefficients in these examples are illustrative.
 See the FluidDynamics and Hydrodynamics DocC guides for coverage, numerical cost,
 pressure-reference rules and the remaining real-gas/free-surface fidelity work.
+
+
+Single finite-burn position targeting is available through
+`FiniteBurnPositionTargeting.refine`. It keeps ignition, thrust and exhaust speed
+fixed, correcting duration and inertial direction through numerical endpoint
+sensitivities and bounded backtracking. The result reports position convergence,
+velocity residual, propellant, evaluation count and rejection count. Velocity is
+**not** a targeted constraint; callers must separately verify it along with
+clearance, attitude and reserve. Every trial uses `ScheduledThrustPropagator` and
+caller-supplied force models; no gravity or physical origin is assumed implicitly.
+
+`FiniteBurnStateTargeting.refine` corrects two fixed-ignition finite burns against
+both terminal position and velocity using six duration-vector controls, central
+differences, pivoted elimination and bounded backtracking/evaluation count.
+It returns the best evaluated schedule, including on nonconvergence: callers must
+check both residuals and independently validate clearance, resources and operations.

Sources/Atmospheres/EarthEntryAtmosphere.swift

@@ -0,0 +1,40 @@
+// SPDX-License-Identifier: Apache-2.0
+// Copyright 2026 Mattias Holm
+import Foundation
+
+/// ISA lower layers joined continuously to a sparse US Standard Atmosphere 1976
+/// upper table. Log-linear pressure/density, linear kinetic temperature. Fixed
+/// nominal air, not weather/solar activity, chemistry or continuum certification.
+public struct EarthEntryAtmosphere: AtmosphericProfile {
+    public let altitudeRangeMetres = -500.0...300000.0
+    public let altitudeDatum = "Geometric metres above the simulation mean sea level"
+    public let source = "https://www.ngdc.noaa.gov/stp/space-weather/online-publications/miscellaneous/us-standard-atmosphere-1976/us-standard-atmosphere_st76-1562_noaa.pdf"
+    public let assumptions = "ISA to its explicit upper limit; USSA1976 Table I pp.68–70 at 90–300 km. Sparse log interpolation; match ISA at join. Upper sound speed uses gamma=1.4 and p/rho; Sutherland viscosity is a nominal transport approximation, not resolved composition."
+    public init() {}
+    public func properties(at altitudeMetres: Double) throws -> AtmosphericProperties {
+        guard altitudeMetres.isFinite, altitudeRangeMetres.contains(altitudeMetres) else { throw AtmosphericModelError.altitudeOutsideCoverage }
+        let isa = ISAAtmosphere(), join = isa.altitudeRangeMetres.upperBound
+        if altitudeMetres <= join { return try isa.properties(at:altitudeMetres) }
+        let base = try isa.properties(at:join)
+        // Geometric km, kinetic K, pressure Pa (publication mb ×100), rho kg/m³.
+        let rows: [(Double,Double,Double,Double)] = [
+            (join,base.temperatureKelvin,base.pressurePascals,base.densityKilogramsPerCubicMetre),
+            (90000,186.87,0.18359,3.416e-6), (95000,188.42,0.075966,1.393e-6),
+            (100000,195.08,0.032011,5.604e-7), (105000,208.84,0.014477,2.325e-7),
+            (110000,240,0.0071042,9.708e-8), (115000,300,0.0040075,4.061e-8),
+            (120000,360,0.0025382,2.222e-8), (125000,417.23,0.0017340,1.291e-8),
+            (130000,469.27,0.0012505,8.152e-9),
+            (140000,559.63,0.00072028,3.831e-9), (150000,634.39,0.00045422,2.076e-9),
+            (160000,696.29,0.00030395,1.233e-9), (170000,747.57,0.00021210,7.815e-10),
+            (180000,790.07,0.00015271,5.194e-10), (190000,825.31,0.00011266,3.581e-10),
+            (200000,854.56,0.000084736,2.541e-10), (220000,899.01,0.000050149,1.367e-10),
+            (230000,915.78,0.000039276,1.029e-10), (240000,929.73,0.000031059,7.858e-11),
+            (250000,941.33,0.000024767,6.073e-11), (300000,976.00,0.0000087704,1.916e-11)]
+        let i = min(rows.count-2,rows.lastIndex { $0.0 <= altitudeMetres }!)
+        let a = rows[i], b = rows[i+1], u = (altitudeMetres-a.0)/(b.0-a.0)
+        let t = a.1+(b.1-a.1)*u, p = exp(log(a.2)+(log(b.2)-log(a.2))*u)
+        let rho = exp(log(a.3)+(log(b.3)-log(a.3))*u)
+        return AtmosphericProperties(pressurePascals:p,temperatureKelvin:t,densityKilogramsPerCubicMetre:rho,
+            speedOfSoundMetresPerSecond:sqrt(1.4*p/rho),dynamicViscosityPascalSeconds:1.458e-6*pow(t,1.5)/(t+110.4))
+    }
+}

Sources/CelMek/FiniteBurnPositionTargeting.swift

@@ -0,0 +1,100 @@
+// SPDX-License-Identifier: Apache-2.0
+import Foundation
+import Math
+
+public enum FiniteBurnPositionTargetingStatus: Equatable {
+    case positionConverged, evaluationLimitReached, stepLimitReached, singularSensitivity
+}
+public struct FiniteBurnPositionTargetingResult {
+    public let status: FiniteBurnPositionTargetingStatus
+    public let evaluations, rejectedTrials: Int
+    public let initial, best: FiniteTransferEvaluation
+}
+
+/// Differential correction for one fixed-ignition, constant inertial burn.
+/// Adjusts duration and direction using a numerically differentiated position
+/// sensitivity and bounded backtracking. Velocity is diagnostic, NOT targeted.
+/// Callers must independently verify clearance, attitude, resources and velocity.
+public enum FiniteBurnPositionTargeting {
+    public static func refine(initialState: SpacecraftState, initialBurn: FiniteBurn,
+        target: CartesianState, context: ForceEvaluationContext, backgroundModels: [any ForceModel] = [],
+        positionToleranceMetres: Double = 1, minimumRemainingPropellantKilograms: Double = 0,
+        maximumEvaluations: Int = 120, burnWindow: ClosedRange<Epoch>? = nil,
+        numericalOptions: NumericalPropagationOptions = .preset(for:.high)) throws -> FiniteBurnPositionTargetingResult {
+        guard target.epoch > initialState.epoch, initialBurn.start >= initialState.epoch,
+            initialBurn.end <= target.epoch, initialBurn.thrust.thrustNewtons > 0,
+            positionToleranceMetres.isFinite, positionToleranceMetres > 0,
+            minimumRemainingPropellantKilograms.isFinite, minimumRemainingPropellantKilograms >= 0,
+            (1...10_000).contains(maximumEvaluations) else { throw FiniteTransferPlanningError.invalidInput }
+        guard target.frame === initialState.frame, initialBurn.thrust.frame === initialState.frame,
+            context.frame === initialState.frame else { throw SpacecraftDynamicsError.frameMismatch }
+        let thrust = initialBurn.thrust
+        var evaluations = 0, rejected = 0
+        func norm(_ v: SIMD3<Double>) -> Double { hypot(hypot(v.x,v.y),v.z) }
+        func cross(_ a: SIMD3<Double>, _ b: SIMD3<Double>) -> SIMD3<Double> {
+            SIMD3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x)
+        }
+        func evaluate(_ control: SIMD3<Double>) throws -> FiniteTransferEvaluation? {
+            let seconds = norm(control)
+            guard seconds.isFinite, seconds > 1e-8 else { rejected += 1; return nil }
+            let end = try initialBurn.start.advanced(by:Duration(seconds:seconds))
+            let used = seconds*thrust.propellantConsumptionKilogramsPerSecond
+            guard end <= target.epoch, initialState.mass.propellantMassKilograms-used > 0,
+                initialState.mass.propellantMassKilograms-used >= minimumRemainingPropellantKilograms,
+                burnWindow.map({ $0.contains(initialBurn.start) && $0.contains(end) }) ?? true else { rejected += 1; return nil }
+            let burn = try FiniteBurn(start:initialBurn.start,end:end,thrust:FiniteThrust(frame:thrust.frame,
+                thrustNewtons:thrust.thrustNewtons,effectiveExhaustVelocityMetresPerSecond:thrust.effectiveExhaustVelocityMetresPerSecond,direction:control))
+            let schedule = try ScheduledThrust(frame:initialState.frame,burns:[burn])
+            evaluations += 1
+            do {
+                let state = try ScheduledThrustPropagator(initialState:initialState,schedule:schedule,
+                    backgroundModels:backgroundModels,context:context,options:numericalOptions).spacecraftState(at:target.epoch)
+                return FiniteTransferEvaluation(schedule:schedule,arrivalState:state,
+                    positionErrorMetres:norm(state.cartesian.position-target.position),
+                    velocityErrorMetresPerSecond:norm(state.cartesian.velocity-target.velocity),
+                    propellantUsedKilograms:initialState.mass.propellantMassKilograms-state.mass.propellantMassKilograms)
+            } catch is NumericalPropagationError { rejected += 1; return nil }
+        }
+        var control = thrust.direction * (try initialBurn.end.duration(since:initialBurn.start).seconds)
+        guard let initial = try evaluate(control) else { throw FiniteTransferPlanningError.invalidInput }
+        var best = initial, status = FiniteBurnPositionTargetingStatus.evaluationLimitReached
+        while best.positionErrorMetres > positionToleranceMetres {
+            guard evaluations+7 <= maximumEvaluations else { break }
+            let step = max(1e-4,norm(control)*1e-4)
+            var columns: [SIMD3<Double>] = []
+            for axis in 0..<3 {
+                var plus = control, minus = control
+                plus[axis] += step; minus[axis] -= step
+                let a = try evaluate(plus), b = try evaluate(minus)
+                if let a, let b { columns.append((a.arrivalState.cartesian.position-b.arrivalState.cartesian.position)/(2*step)) }
+                else if let a { columns.append((a.arrivalState.cartesian.position-best.arrivalState.cartesian.position)/step) }
+                else if let b { columns.append((best.arrivalState.cartesian.position-b.arrivalState.cartesian.position)/step) }
+                else { break }
+            }
+            guard columns.count == 3 else { status = .stepLimitReached; break }
+            let lengths = columns.map(norm)
+            guard lengths.allSatisfy({$0.isFinite && $0 > 1e-15}) else { status = .singularSensitivity; break }
+            let a = columns[0]/lengths[0], b = columns[1]/lengths[1], c = columns[2]/lengths[2]
+            let determinant = (a*cross(b,c)).sum()
+            guard determinant.isFinite, abs(determinant) > 1e-10 else { status = .singularSensitivity; break }
+            let residual = target.position-best.arrivalState.cartesian.position
+            var correction = SIMD3((residual*cross(b,c)).sum()/lengths[0],
+                (a*cross(residual,c)).sum()/lengths[1],(a*cross(b,residual)).sum()/lengths[2])/determinant
+            let length = norm(correction)
+            guard length.isFinite else { status = .singularSensitivity; break }
+            correction *= min(1,max(0.01,norm(control)*0.2)/max(length,1e-30))
+            var accepted = false, fraction = 1.0
+            for _ in 0..<10 where evaluations < maximumEvaluations {
+                let candidate = control+correction*fraction
+                if let trial = try evaluate(candidate), trial.positionErrorMetres < best.positionErrorMetres {
+                    control = candidate; best = trial; accepted = true; break
+                }
+                fraction *= 0.5
+            }
+            guard accepted else { status = evaluations >= maximumEvaluations ? .evaluationLimitReached : .stepLimitReached; break }
+        }
+        if best.positionErrorMetres <= positionToleranceMetres { status = .positionConverged }
+        return FiniteBurnPositionTargetingResult(status:status,evaluations:evaluations,rejectedTrials:rejected,initial:initial,best:best)
+    }
+}
+