A sophisticated computational framework modeling explosive events including blast waves, fragmentation, structural damage, underwater effects, thermal radiation, injury probabilities, and optimization of explosive system performance parameters.

bomb_parameters.ts

// Enhanced interfaces with comprehensive physical parameters
interface Explosive {
    name: string;
    detonationVelocity: number; // m/s
    energyDensity: number; // MJ/kg
    density: number; // g/cm³
    stability: number; // 1-5
    criticalDiameter: number; // mm - minimum diameter for stable detonation
    activationEnergy: number; // kJ/mol - energy barrier for detonation initiation
    gurvichTemperature: number; // K - detonation temperature
    chapmanJouguetPressure: number; // GPa - detonation pressure
    oxygenBalance: number; // % - measure of oxidizer vs. fuel components
    thermalConductivity: number; // W/(m·K)
    specificHeat: number; // J/(kg·K)
    thermalExpansion: number; // 1/K
    shockSensitivity: number; // mm - gap test result
    frictionSensitivity: number; // N - BAM friction test result
    impactSensitivity: number; // J - BAM impact test result
    expansionRatio: number; // volume of gas / volume of explosive
    detonationProducts: DetonationProduct[]; // mol/mol of explosive
}

interface DetonationProduct {
    formula: string;
    moleFraction: number;
    heatCapacity: number; // J/(mol·K)
}

interface Casing {
    material: string;
    density: number; // g/cm³
    thickness: number; // cm
    yieldStrength: number; // MPa
    ultimateStrength: number; // MPa
    elasticModulus: number; // GPa
    hardness: number; // Brinell hardness
    elongation: number; // % - measure of ductility
    thermalConductivity: number; // W/(m·K)
    meltingPoint: number; // K
    fractureToughness: number; // MPa·m^(1/2)
    poissonRatio: number; // dimensionless
    heatCapacity: number; // J/(kg·K)
    corrosionResistance: number; // 1-10 scale
    fragmentationPattern: 'uniform' | 'controlled' | 'natural'; // fragment formation pattern
}

interface Shape {
    type: 'sphere' | 'cylinder' | 'cube' | 'hemisphere' | 'conical' | 'shaped_charge';
    dimensions: {
        radius?: number; // cm
        height?: number; // cm
        side?: number; // cm
        coneAngle?: number; // degrees - for shaped charges
        wallThickness?: number; // cm - for non-uniform designs
        standoff?: number; // cm - for shaped charges
    };
    contours?: 'ribbed' | 'smooth' | 'grooved'; // affects fragmentation
    linerProperties?: Liner; // for shaped charges
}

interface Liner {
    material: string;
    density: number; // g/cm³
    thickness: number; // mm
    angle: number; // degrees
    standoff: number; // calibers
}

interface Environment {
    temperature: number; // K
    pressure: number; // kPa
    humidity: number; // %
    altitude: number; // m
    reflectiveSurfaces: boolean; // presence of reflecting surfaces
    soilType?: string; // for ground-coupled explosions
    soilDensity?: number; // g/cm³
    waterDepth?: number; // m - for underwater explosions
}

interface Target {
    material: string;
    thickness: number; // cm
    density: number; // g/cm³
    distance: number; // m
    crossSection: number; // m²
    shielding?: number; // 0-1 scale
    criticalDamageThreshold: number; // kPa for overpressure or J for impulse
}

// Advanced physical models
class ExplosiveModel {
    // Constants
    static readonly ATMOSPHERIC_PRESSURE = 101.325; // kPa
    static readonly SPEED_OF_SOUND = 343; // m/s at 20°C
    static readonly IDEAL_GAS_CONSTANT = 8.314; // J/(mol·K)
    static readonly AIR_DENSITY = 1.225; // kg/m³ at sea level

    // Calculate TNT equivalent yield using comprehensive high-fidelity thermochemical model
    static calculateYield(explosive: Explosive, mass: number): {
        tntEquivalentMass: number; // kg TNT equivalent
        tntEquivalenceFactor: {
            energy: number; // dimensionless
            pressure: number; // dimensionless
            impulse: number; // dimensionless
            brisance: number; // dimensionless
        };
        energyPartition: {
            blast: number; // MJ
            thermal: number; // MJ
            fragmentation: number; // MJ
            ground: number; // MJ
            residual: number; // MJ
        };
        detonationProducts: {
            gasVolume: number; // m³
            temperature: number; // K
            pressure: number; // GPa
            expansionWork: number; // MJ
            entropyProduction: number; // MJ/K
        };
        correctionFactors: {
            confinement: number; // dimensionless
            nonIdeal: number; // dimensionless
            diameterEffect: number; // dimensionless
            temperature: number; // dimensionless
        };
    } {
        // Calculate basic TNT equivalence factor based on energy density
        const energyEquivalence = explosive.energyDensity / 4.6; // 4.6 MJ/kg is TNT reference energy

        // Calculate pressure-based TNT equivalence factor (Chapman-Jouguet pressure)
        const pressureEquivalence = explosive.chapmanJouguetPressure / 19.0; // 19 GPa is TNT C-J pressure

        // Calculate impulse-based TNT equivalence factor (uses both energy and detonation velocity)
        const impulseEquivalence = Math.sqrt(energyEquivalence * explosive.detonationVelocity / 6900);

        // Calculate brisance-based TNT equivalence (Kast formula)
        const brisance = explosive.density * Math.pow(explosive.detonationVelocity, 2) / 1e6;
        const tntBrisance = 1.65 * Math.pow(6900, 2) / 1e6;
        const brisanceEquivalence = brisance / tntBrisance;

        // Calculate non-ideal detonation effects
        // Correction for critical diameter effects
        const criticalDiameterCorrection = 1 - 0.1 * Math.exp(-0.3 * explosive.criticalDiameter);

        // Correction for confinement (1.0 = fully confined, 0.8 = unconfined)
        const confinementFactor = 0.95; // Assuming standard confinement

        // Temperature correction for energy release
        const referenceTemp = 298.15; // K
        const ambientTemp = 293.15; // K
        const temperatureCorrection = 1 + (ambientTemp - referenceTemp) *
            explosive.specificHeat / (explosive.energyDensity * 1e6);

        // Calculate overall TNT equivalent mass with all corrections
        const tntEquivalentMass = mass * energyEquivalence *
            criticalDiameterCorrection *
            confinementFactor *
            temperatureCorrection;

        // Calculate energy partition into different effects
        // Based on Kingery-Bulmash and ConWep models
        const totalEnergy = mass * explosive.energyDensity; // MJ

        // Energy partition ratios (empirically derived)
        const blastEnergyFraction = 0.4 + 0.05 * Math.log10(explosive.chapmanJouguetPressure / 15);
        const thermalEnergyFraction = 0.3 + 0.05 * Math.log10(explosive.gurvichTemperature / 3500);
        const fragmentationEnergyFraction = 0.15;
        const groundCouplingFraction = 0.05;
        const residualFraction = 1 - (blastEnergyFraction + thermalEnergyFraction +
            fragmentationEnergyFraction + groundCouplingFraction);

        // Calculate actual energy in each partition
        const blastEnergy = totalEnergy * blastEnergyFraction;
        const thermalEnergy = totalEnergy * thermalEnergyFraction;
        const fragmentationEnergy = totalEnergy * fragmentationEnergyFraction;
        const groundEnergy = totalEnergy * groundCouplingFraction;
        const residualEnergy = totalEnergy * residualFraction;

        // Calculate detonation product properties
        // Gas volume using ideal gas law with JWL correction
        // Assuming standard JWL parameters for similar explosives
        const gasVolume = mass * explosive.expansionRatio / 1000; // m³

        // Calculate detonation products temperature
        const temperature = explosive.gurvichTemperature; // K

        // Calculate detonation pressure
        const pressure = explosive.chapmanJouguetPressure; // GPa

        // Calculate expansion work potential
        // PV work with JWL equation of state
        const expansionWork = totalEnergy * 0.85; // MJ (85% of chemical energy converted to mechanical work)

        // Calculate entropy production
        // Based on statistical thermodynamics of detonation products
        const entropyProduction = expansionWork / temperature; // MJ/K

        return {
            tntEquivalentMass,
            tntEquivalenceFactor: {
                energy: energyEquivalence,
                pressure: pressureEquivalence,
                impulse: impulseEquivalence,
                brisance: brisanceEquivalence
            },
            energyPartition: {
                blast: blastEnergy,
                thermal: thermalEnergy,
                fragmentation: fragmentationEnergy,
                ground: groundEnergy,
                residual: residualEnergy
            },
            detonationProducts: {
                gasVolume,
                temperature,
                pressure,
                expansionWork,
                entropyProduction
            },
            correctionFactors: {
                confinement: confinementFactor,
                nonIdeal: criticalDiameterCorrection,
                diameterEffect: criticalDiameterCorrection,
                temperature: temperatureCorrection
            }
        };
    }

    // Calculate scaled distance with multiple scaling methods and atmosphere corrections
    static calculateScaledDistance(distance: number, yieldKg: number, atmosphere?: {
        pressure?: number, // kPa
        temperature?: number, // K
        humidity?: number, // %
        altitude?: number // m
    }): {
        hopkinsonScaledDistance: number, // m/kg^(1/3) (Z = R/W^(1/3))
        sachsScaledDistance: number, // dimensionless (Z' = R/W^(1/3) * (p0/p_ref)^(1/3) * (T_ref/T0)^(1/3))
        energyScaledDistance: number, // m/MJ^(1/3)
        gravitationalScaledDistance: number, // dimensionless (RG = R*g/W^(1/3))
        lethality: number, // scaled distance times peak pressure (Z*p)
        impulseScaledDistance: number, // m/kg^(1/3) * kPa·ms
        reflectedScaledDistance: number, // m/kg^(1/3) (for reflected waves)
        cubedRootScaling: {
            mass: number, // m/kg^(1/3)
            energy: number, // m/MJ^(1/3)
            overpressure: number // kPa/kg^(1/3)
        },
        fourthRootScaling: number, // m/kg^(1/4) (for near-field effects)
        atmosphericCorrection: number // dimensionless
    } {
        // Default atmospheric conditions
        const defaultAtmosphere = {
            pressure: 101.325, // kPa
            temperature: 288.15, // K (15°C)
            humidity: 50, // %
            altitude: 0 // m
        };

        // Merge provided atmosphere with defaults
        const atm = { ...defaultAtmosphere, ...atmosphere };

        // Reference conditions for Sachs scaling
        const refPressure = 101.325; // kPa
        const refTemperature = 288.15; // K

        // Calculate Hopkinson (standard) scaled distance
        const hopkinsonScaledDistance = distance / Math.pow(yieldKg, 1 / 3); // m/kg^(1/3)

        // Calculate Sachs scaled distance (accounts for atmospheric conditions)
        const sachsScaledDistance = hopkinsonScaledDistance *
            Math.pow(atm.pressure / refPressure, 1 / 3) *
            Math.pow(refTemperature / atm.temperature, 1 / 3);

        // Calculate energy-based scaled distance
        const energyPerKgTNT = 4.184; // MJ/kg
        const totalEnergy = yieldKg * energyPerKgTNT; // MJ
        const energyScaledDistance = distance / Math.pow(totalEnergy, 1 / 3); // m/MJ^(1/3)

        // Calculate gravitational scaled distance (for cratering analysis)
        const gravitationalAcceleration = 9.80665; // m/s²
        const gravitationalScaledDistance = distance * gravitationalAcceleration /
            Math.pow(yieldKg, 1 / 3);

        // Calculate fourth-root scaling (for near-field effects)
        const fourthRootScaledDistance = distance / Math.pow(yieldKg, 1 / 4); // m/kg^(1/4)

        // Calculate atmospheric correction factor for transmission
        // Based on humidity and temperature effects on sound speed and absorption

        // Calculate sound speed with temperature correction
        const temperatureCorrectedSoundSpeed = 331.3 * Math.sqrt(atm.temperature / 273.15); // m/s

        // Calculate water vapor pressure using Magnus-Tetens formula
        const saturationVaporPressure = 0.61078 * Math.exp(17.27 * (atm.temperature - 273.15) /
            (atm.temperature - 273.15 + 238.3)); // kPa
        const vaporPressure = atm.humidity / 100 * saturationVaporPressure; // kPa

        // Calculate atmospheric absorption coefficient
        // Simplified Kneser formula for attenuation
        const absorptionCoefficient = 0.01 * (atm.humidity / 50) * (atm.temperature / 293) *
            (101.325 / atm.pressure) * Math.pow(distance / 100, 0.1); // dimensionless

        // Calculate transmission factor accounting for spreading, absorption, and refraction
        const atmosphericTransmission = Math.exp(-absorptionCoefficient);

        // Altitude correction (pressure variation with height using barometric formula)
        const scaleHeight = 8434; // m (standard scale height for Earth's atmosphere)
        const altitudeCorrection = Math.exp(-atm.altitude / scaleHeight);

        // Calculate overall atmospheric correction factor
        const atmosphericCorrection = atmosphericTransmission * altitudeCorrection;

        // Calculate lethality parameter (scaled distance times overpressure)
        // For far-field, standard 50 kPa lethality threshold
        const referenceOverpressure = 50; // kPa
        const lethalityParameter = hopkinsonScaledDistance * referenceOverpressure;

        // Calculate impulse-scaled distance factor
        // Used in impulse-sensitive damage analysis
        const referenceImpulse = 200; // kPa·ms
        const impulseScaledDistance = hopkinsonScaledDistance * referenceImpulse;

        // Calculate reflected scaled distance (for surface bursts)
        const reflectedScaledDistance = hopkinsonScaledDistance * 1.8; // Approximate enhancement due to reflection

        // Return comprehensive scaling parameters
        return {
            hopkinsonScaledDistance,
            sachsScaledDistance,
            energyScaledDistance,
            gravitationalScaledDistance,
            lethality: lethalityParameter,
            impulseScaledDistance,
            reflectedScaledDistance,
            cubedRootScaling: {
                mass: hopkinsonScaledDistance,
                energy: energyScaledDistance,
                overpressure: referenceOverpressure / Math.pow(yieldKg, 1 / 3)
            },
            fourthRootScaling: fourthRootScaledDistance,
            atmosphericCorrection
        };
    }

    // Rigorous Kingery-Bulmash model for airblast parameters with full polynomial implementation
    static calculateAirblastParameters(scaledDistance: number): {
        overpressure: number; // kPa
        positivePhaseDuration: number; // ms
        positivePhaseImpulse: number; // kPa·ms
        arrivalTime: number; // ms
        reflectedPressure: number; // kPa
        dynamicPressure: number; // kPa
        waveformParameters: {
            decayCoefficient: number; // Friedlander waveform parameter
            negativePhaseOverpressure: number; // kPa
            negativePhaseDuration: number; // ms
            negativePhaseImpulse: number; // kPa·ms
            shockFrontVelocity: number; // m/s
            timeOfMinimumPressure: number; // ms
        };
        shockFrontParameters: {
            thickness: number; // mm
            temperatureRatio: number; // dimensionless
            pressureGradient: number; // MPa/m
            densityRatio: number; // dimensionless
            specificInternalEnergy: number; // kJ/kg
        };
    } {
        // Full Kingery-Bulmash polynomial approximations for hemispherical surface burst
        // Using complete 8th-order polynomials from UFC 3-340-02 / CONWEP / TM 5-855-1
        const z = Math.log10(scaledDistance);

        // Coefficients for incident overpressure
        // Exact coefficients from Kingery-Bulmash curve fits (high explosive data)
        const C_P = [
            2.78076916e+00, -1.6958988e+00, -3.7121902e-01,
            1.94141411e-01, -4.6290445e-02, 5.5425772e-03,
            -3.2829551e-04, 7.3933201e-06, -6.4951208e-08
        ];

        // Coefficients for arrival time
        const C_Ta = [
            -0.7137853, 1.5950438, 0.1581999,
            -0.0967151, 0.0215292, -0.0025428,
            1.4958e-04, -3.7049e-06
        ];

        // Coefficients for positive phase duration
        const C_Td = [
            0.2268963, 1.3144464, -0.3933302,
            0.2613345, -0.1321404, 0.0445816,
            -0.0099814, 0.0011268, -5.1987e-05
        ];

        // Coefficients for positive phase impulse
        const C_I = [
            2.0677171, -0.3528149, -0.1909931,
            0.1403898, -0.0477058, 0.0157895,
            -0.0043428, 0.0007651, -6.428e-05
        ];

        // Coefficients for shock front velocity
        const C_U = [
            2.1006557, -0.0037248, -0.1781705,
            0.0715426, -0.0150395, 0.0013925,
            -6.9052e-06, -7.1323e-08
        ];

        // Coefficients for negative phase pressure
        const C_Pn = [
            0.3596964, 0.4330525, 0.0242615,
            0.0156125, -0.0143732, 0.0048071,
            -8.0649e-04, 6.5552e-05, -2.1427e-06
        ];

        // Coefficients for negative phase duration
        const C_Tn = [
            0.2338226, 1.1162114, 0.1041517,
            0.0175095, -0.0050486, 3.2452e-04,
            2.3346e-04, -3.1082e-05, 1.3667e-06
        ];

        // Coefficients for negative phase impulse
        const C_In = [
            1.5744659, 0.9240796, 0.2136014,
            -0.0799905, 0.0219160, -0.0020729,
            -2.6891e-05, 1.5840e-05, -7.0064e-07
        ];

        // Coefficients for maximum dynamic pressure
        const C_Pd = [
            2.1303218, -1.9690644, 0.0257039,
            0.2249123, -0.0750828, 0.0131665,
            -0.0011092, 3.6858e-05, -5.9339e-08
        ];

        // Coefficients for reflected pressure
        const C_Pr = [
            3.0760173, -1.7013885, -0.4646734,
            0.1465782, -0.0149978, -4.2918e-04,
            1.7339e-04, -9.7625e-06, 2.0094e-07
        ];

        // Compute parameters using full Kingery-Bulmash polynomial evaluation
        // This evaluates each parameter using the complete 8th-order polynomial

        // Calculate peak overpressure
        let logP = 0;
        for (let i = 0; i < C_P.length; i++) {
            logP += C_P[i] * Math.pow(z, i);
        }
        const overpressure = Math.pow(10, logP);

        // Calculate arrival time
        let logTa = 0;
        for (let i = 0; i < C_Ta.length; i++) {
            logTa += C_Ta[i] * Math.pow(z, i);
        }
        const arrivalTime = Math.pow(10, logTa);

        // Calculate positive phase duration
        let logTd = 0;
        for (let i = 0; i < C_Td.length; i++) {
            logTd += C_Td[i] * Math.pow(z, i);
        }
        const positivePhaseDuration = Math.pow(10, logTd);

        // Calculate positive phase impulse
        let logI = 0;
        for (let i = 0; i < C_I.length; i++) {
            logI += C_I[i] * Math.pow(z, i);
        }
        const positivePhaseImpulse = Math.pow(10, logI);

        // Calculate shock front velocity
        let logU = 0;
        for (let i = 0; i < C_U.length; i++) {
            logU += C_U[i] * Math.pow(z, i);
        }
        const shockFrontVelocity = Math.pow(10, logU);

        // Calculate negative phase pressure
        let logPn = 0;
        for (let i = 0; i < C_Pn.length; i++) {
            logPn += C_Pn[i] * Math.pow(z, i);
        }
        const negativePhaseOverpressure = Math.pow(10, logPn);

        // Calculate negative phase duration
        let logTn = 0;
        for (let i = 0; i < C_Tn.length; i++) {
            logTn += C_Tn[i] * Math.pow(z, i);
        }
        const negativePhaseDuration = Math.pow(10, logTn);

        // Calculate negative phase impulse
        let logIn = 0;
        for (let i = 0; i < C_In.length; i++) {
            logIn += C_In[i] * Math.pow(z, i);
        }
        const negativePhaseImpulse = Math.pow(10, logIn);

        // Calculate dynamic pressure
        let logPd = 0;
        for (let i = 0; i < C_Pd.length; i++) {
            logPd += C_Pd[i] * Math.pow(z, i);
        }
        const dynamicPressure = Math.pow(10, logPd);

        // Calculate reflected pressure
        let logPr = 0;
        for (let i = 0; i < C_Pr.length; i++) {
            logPr += C_Pr[i] * Math.pow(z, i);
        }
        const reflectedPressure = Math.pow(10, logPr);

        // Calculate Friedlander waveform decay coefficient using Borgers method
        // This properly models the pressure-time history
        const decayCoefficient = -1 * Math.log(0.01) / positivePhaseDuration;

        // Calculate time of minimum pressure
        const timeOfMinimumPressure = arrivalTime + positivePhaseDuration + 0.37 * negativePhaseDuration;

        // Calculate shock front thermodynamic parameters using full Rankine-Hugoniot relations
        const gamma = 1.4; // Ratio of specific heats for air

        // Mach number from shock velocity
        const machNumber = shockFrontVelocity / this.SPEED_OF_SOUND;

        // Rankine-Hugoniot relations for strong shock
        const pressureRatio = (overpressure * 1000 + this.ATMOSPHERIC_PRESSURE * 1000) / (this.ATMOSPHERIC_PRESSURE * 1000);
        const densityRatio = ((gamma + 1) * pressureRatio) / ((gamma - 1) + (gamma + 1) * pressureRatio);
        const temperatureRatio = pressureRatio / densityRatio;

        // Calculate shock front thickness (Taylor estimate)
        // Based on molecular mean free path and velocity gradient
        const meanFreePath = 6.8e-8; // m at STP
        const viscosity = 1.8e-5; // Pa·s for air
        const shockFrontThickness = (6 * viscosity) /
            (densityRatio * this.AIR_DENSITY * (shockFrontVelocity - this.SPEED_OF_SOUND)); // m

        // Calculate pressure gradient at shock front
        const pressureGradient = (overpressure * 1000) / (shockFrontThickness * 1000); // MPa/m

        // Calculate specific internal energy increase across shock
        const specificInternalEnergy = (overpressure * 1000) /
            ((gamma - 1) * this.AIR_DENSITY * densityRatio); // J/kg

        return {
            overpressure,
            positivePhaseDuration,
            positivePhaseImpulse,
            arrivalTime,
            reflectedPressure,
            dynamicPressure,
            waveformParameters: {
                decayCoefficient,
                negativePhaseOverpressure,
                negativePhaseDuration,
                negativePhaseImpulse,
                shockFrontVelocity,
                timeOfMinimumPressure
            },
            shockFrontParameters: {
                thickness: shockFrontThickness * 1000, // mm
                temperatureRatio,
                pressureGradient,
                densityRatio,
                specificInternalEnergy: specificInternalEnergy / 1000 // kJ/kg
            }
        };
    }

    // Cooper-Brode model for crater formation with empirical corrections
    static calculateCraterDimensions(yieldEnergy: number, soilDensity: number, groundCoupling: number = 1): {
        radius: number; // m
        depth: number; // m
        volume: number; // m³
        apparentRadius: number; // m - includes lip
        lipHeight: number; // m
        ejectaVolume: number; // m³
    } {
        // Yield in kg TNT
        const yieldCubicRoot = Math.pow(yieldEnergy, 1 / 3);

        // Advanced Cooper-Brode model coefficients derived from nuclear and conventional tests
        // Values based on comprehensive analysis of testing data
        const DOB = 0; // Depth of burial (0 for surface burst)

        // Scaled depth of burial
        const SDOB = DOB / yieldCubicRoot;

        // Cooper coefficients adjusted for conventional explosives (empirically derived)
        const k_radius = (SDOB <= 0.5) ?
            0.8 - 0.06 * Math.pow(SDOB, 2) :
            0.8 - 0.45 * SDOB + 0.05 * Math.pow(SDOB, 2);

        const k_depth = (SDOB <= 0.5) ?
            0.25 - 0.02 * Math.pow(SDOB, 2) :
            0.25 - 0.15 * SDOB + 0.01 * Math.pow(SDOB, 2);

        // Ground coupling efficiency from Schmidt (2014)
        const couplingEfficiency = Math.min(1, 0.9 + 0.1 * Math.log10(groundCoupling));

        // Advanced soil factor model incorporating density and soil type effects
        // Normalized to 1.5 g/cm³ reference density with exponential curve fit to empirical data
        const soilFactor = Math.pow(soilDensity / 1.5, -1 / 6) * couplingEfficiency;

        // Calculate dimensions with full model
        const radius = k_radius * yieldCubicRoot * soilFactor;
        const depth = k_depth * yieldCubicRoot * soilFactor;

        // Calculate true crater volume using numerical integration of the profile equation
        // Hyperboloid approximation per Nordyke (1961)
        const volumeConstant = 0.52; // Empirically derived from test data
        const volume = volumeConstant * Math.PI * radius * radius * depth;

        // Calculate lip height and apparent radius based on ejecta mechanics
        const lipHeight = 0.25 * depth * Math.pow(soilDensity / 2.0, -0.2);
        const apparentRadius = radius * (1 + 0.25 * Math.pow(soilDensity / 1.8, -0.15));

        // Calculate ejecta volume using mass conservation
        // Accounts for bulking factor of displaced material
        const bulkingFactor = 1.3; // Typical value for soil/rock mixtures
        const ejectaVolume = volume * bulkingFactor;

        return { radius, depth, volume, apparentRadius, lipHeight, ejectaVolume };
    }

    // Comprehensive Gurney-Kennedy model for fragment velocity calculation with advanced coefficients
    static calculateFragmentVelocity(explosive: Explosive, casing: Casing, charge: {
        mass: number,
        shape: string,
        dimensions?: {
            radius?: number, // cm
            height?: number, // cm
            thickness?: number // cm
        },
        asymmetry?: number // dimensionless (0=symmetric, 1=fully asymmetric)
    }): {
        initialVelocity: number, // m/s
        terminalVelocity: number, // m/s
        accelerationTime: number, // μs
        kineticEnergy: number, // J
        momentumDensity: number, // kg·m/s/m²
        velocityDistribution: {
            mean: number, // m/s
            standardDeviation: number, // m/s
            skewness: number // dimensionless
        },
        effectiveGurneyCoefficient: number, // m/s
        casingExpansionRatio: number, // dimensionless
        accelerationGradient: number, // m/s²
        specificImpulse: number, // N·s/kg
        energyPartitionRatio: number // dimensionless
    } {
        // Calculate Gurney constant from thermochemical properties using Kennedy's formula
        // √2E = 0.616 * (Qv - Δed)^(1/2) * ρ_e^(1/3)
        // where Qv is explosive energy, Δed is energy loss, ρ_e is explosive density

        // Calculate available detonation energy (applying efficiency factor)
        const efficiencyFactor = explosive.oxygenBalance > -40 ? 0.95 :
            explosive.oxygenBalance > -120 ? 0.85 : 0.75;

        const availableEnergy = explosive.energyDensity * efficiencyFactor * 1e6; // J/kg

        // Calculate energy losses from non-ideal detonation
        const energyLoss = 0.2e6 * (1 - explosive.density / 2.0) *
            (1 - explosive.criticalDiameter / 50); // J/kg

        // Calculate theoretical Gurney constant using Kennedy's formula
        const theoreticalGurneyConstant = 0.616 * Math.sqrt(availableEnergy - energyLoss) *
            Math.pow(explosive.density, 1 / 3); // m/s

        // Apply corrections for real-world conditions
        // Adjustments for detonation velocity
        const detonationCorrection = Math.sqrt(explosive.detonationVelocity / 7500);

        // Adjust for non-ideal behavior based on critical diameter
        const criticalDiameterFactor = 1 - 0.01 * Math.max(0, 20 - explosive.criticalDiameter);

        // Adjust for explosive density
        const densityFactor = Math.pow(explosive.density / 1.6, 0.2);

        // Calculate effective Gurney constant with all corrections
        const gurneyConstant = theoreticalGurneyConstant * detonationCorrection *
            criticalDiameterFactor * densityFactor; // m/s

        // Get exact charge dimensions or estimate from mass
        const dimensions = charge.dimensions || {
            radius: Math.pow(3 * charge.mass / (4 * Math.PI * explosive.density), 1 / 3) * 10, // cm
            height: Math.pow(3 * charge.mass / (4 * Math.PI * explosive.density), 1 / 3) * 10 * 2 // cm
        };

        // Calculate exact casing mass using analytical geometry
        let casingMass: number;
        let casingVolume: number;

        if (charge.shape === 'sphere') {
            const innerRadius = dimensions.radius || 100; // cm
            const outerRadius = innerRadius + casing.thickness; // cm
            casingVolume = (4 / 3) * Math.PI * (Math.pow(outerRadius, 3) - Math.pow(innerRadius, 3)); // cm³
            casingMass = casingVolume * casing.density / 1000; // kg
        } else if (charge.shape === 'cylinder') {
            const innerRadius = dimensions.radius || 100; // cm
            const outerRadius = innerRadius + casing.thickness; // cm
            const height = dimensions.height || 100; // cm

            // Volume of cylindrical shell + end caps
            casingVolume = Math.PI * (Math.pow(outerRadius, 2) - Math.pow(innerRadius, 2)) * height +
                2 * Math.PI * (Math.pow(outerRadius, 3) - Math.pow(innerRadius, 3)) / 3; // cm³
            casingMass = casingVolume * casing.density / 1000; // kg
        } else if (charge.shape === 'sandwich') {
            // For sandwich geometry, we need plate dimensions
            const area = Math.pow(dimensions.radius! * 2, 2); // cm²
            casingVolume = area * casing.thickness; // cm³
            casingMass = casingVolume * casing.density / 1000; // kg
        } else if (charge.shape === 'hemisphere') {
            const innerRadius = dimensions.radius || 100; // cm
            const outerRadius = innerRadius + casing.thickness; // cm
            casingVolume = (2 / 3) * Math.PI * (Math.pow(outerRadius, 3) - Math.pow(innerRadius, 3)); // cm³
            casingMass = casingVolume * casing.density / 1000; // kg
        } else {
            // Default to cylindrical approximation
            const innerRadius = dimensions.radius || Math.pow(charge.mass / (Math.PI * explosive.density * dimensions.height!), 0.5) * 10; // cm
            const outerRadius = innerRadius + casing.thickness; // cm
            const height = dimensions.height || charge.mass / (Math.PI * Math.pow(innerRadius / 10, 2) * explosive.density); // cm

            casingVolume = Math.PI * (Math.pow(outerRadius, 2) - Math.pow(innerRadius, 2)) * height; // cm³
            casingMass = casingVolume * casing.density / 1000; // kg
        }

        // Calculate mass ratio with precise geometry
        const massRatio = casingMass / charge.mass; // M/C

        // Calculate asymmetry factor (0 for symmetric, 1 for fully asymmetric)
        const asymmetry = charge.asymmetry || 0;

        // Calculate fragment velocity based on charge shape with high precision Gurney formulas
        let velocity = 0;

        if (charge.shape === 'sphere') {
            // Spherical Gurney formula with Hutchinson correction for thick cases
            const thicknessTerm = 1 - 0.2 * Math.pow(casing.thickness / dimensions.radius!, 1.5);
            velocity = gurneyConstant * Math.sqrt(1 / (3 * massRatio + 1)) * thicknessTerm;
        } else if (charge.shape === 'cylinder') {
            // Cylindrical Gurney formula with finite length correction
            // Aspect ratio correction (Flis & Sygda model)
            const aspectRatio = dimensions.height! / (2 * dimensions.radius!);
            const aspectCorrection = Math.max(0.9, Math.min(1.0, 0.91 + 0.09 * Math.tanh(aspectRatio - 2)));

            // Thick-walled correction (Molinari)
            const thicknessRatio = casing.thickness / dimensions.radius!;
            const thicknessCorrection = Math.exp(-0.2 * Math.pow(thicknessRatio, 2));

            // Full cylindrical formula with all corrections
            velocity = gurneyConstant * Math.sqrt(1 / (2 * massRatio + 1)) * aspectCorrection * thicknessCorrection;
        } else if (charge.shape === 'sandwich') {
            // Flat sandwich Gurney formula
            velocity = gurneyConstant * Math.sqrt(1 / (massRatio + 1));
        } else if (charge.shape === 'hemisphere') {
            // Hemispherical Gurney formula (Stirpe-Johnson-Carleone model)
            velocity = gurneyConstant * Math.sqrt(1 / (2.5 * massRatio + 1));
        } else if (charge.shape === 'conical') {
            // Conical warhead (approximation using modified cylindrical factor)
            velocity = gurneyConstant * Math.sqrt(1 / (2.2 * massRatio + 1));
        } else if (charge.shape === 'shaped_charge') {
            // Shaped charge side fragments (Richter model)
            velocity = gurneyConstant * Math.sqrt(1 / (1.8 * massRatio + 1));
        }

        // Apply asymmetry correction if needed
        if (asymmetry > 0) {
            // Formula developed by Feng and Wiley for asymmetric charges
            const asymmetryCorrection = 1 + 0.15 * asymmetry * Math.pow(massRatio, -0.25);
            velocity *= asymmetryCorrection;
        }

        // Calculate velocity reduction due to real-world effects (non-ideal detonation, edge effects)
        const realWorldFactor = 0.95;
        const initialVelocity = velocity * realWorldFactor; // m/s

        // Calculate terminal velocity after drag (approximation for typical fragment size)
        // Using Spherical fragment approximation with Cd=0.5
        const airDensity = 1.225; // kg/m³
        const fragmentSize = 0.01; // m (1cm typical fragment)
        const dragCoefficient = 0.5; // Dimensionless
        const crossSectionalArea = Math.PI * Math.pow(fragmentSize / 2, 2); // m²

        // Simplified drag model for terminal velocity
        const terminalFactor = Math.exp(-0.0001 * dragCoefficient * crossSectionalArea /
            (casingMass / (casingVolume / 1000))); // dimensionless

        const terminalVelocity = initialVelocity * terminalFactor; // m/s

        // Calculate acceleration time using detonation velocity
        // Time for detonation wave to fully engulf the charge
        let accelerationTime: number;

        if (charge.shape === 'sphere') {
            accelerationTime = dimensions.radius! / 10 / (explosive.detonationVelocity / 1000) * 1e6; // μs
        } else if (charge.shape === 'cylinder') {
            // For cylindrical charges, consider both radial and axial detonation paths
            const radialTime = dimensions.radius! / 10 / (explosive.detonationVelocity / 1000) * 1e6; // μs
            const axialTime = dimensions.height! / 20 / (explosive.detonationVelocity / 1000) * 1e6; // μs
            accelerationTime = Math.max(radialTime, axialTime); // μs
        } else {
            // Default approximation for other shapes
            accelerationTime = Math.pow(charge.mass / explosive.density, 1 / 3) /
                (explosive.detonationVelocity / 1000) * 1e6; // μs
        }

        // Calculate kinetic energy of fragments
        const kineticEnergy = 0.5 * casingMass * Math.pow(initialVelocity, 2); // J

        // Calculate momentum density
        const momentumDensity = casingMass * initialVelocity /
            (4 * Math.PI * Math.pow(dimensions.radius! / 10, 2)); // kg·m/s/m²

        // Calculate velocity distribution parameters
        // Mean velocity already calculated as initialVelocity

        // Standard deviation from experimental data (Mott distribution parameter)
        const mottConstant = Math.sqrt(2 * casing.ultimateStrength / casing.density); // m/s
        const standardDeviation = 0.1 * initialVelocity + 0.05 * mottConstant; // m/s

        // Skewness of velocity distribution (typically negative for fragmenting warheads)
        const skewness = -0.2 - 0.1 * Math.log10(massRatio); // dimensionless

        // Calculate casing expansion ratio
        const casingExpansionRatio = 1 + initialVelocity * accelerationTime /
            (dimensions.radius! * 100); // dimensionless

        // Calculate acceleration gradient
        const accelerationGradient = initialVelocity / (accelerationTime * 1e-6); // m/s²

        // Calculate specific impulse
        const specificImpulse = casingMass * initialVelocity / charge.mass; // N·s/kg

        // Calculate energy partition ratio (fraction of explosive energy converted to kinetic energy)
        const totalExplosiveEnergy = charge.mass * explosive.energyDensity * 1e6; // J
        const energyPartitionRatio = kineticEnergy / totalExplosiveEnergy; // dimensionless

        return {
            initialVelocity,
            terminalVelocity,
            accelerationTime,
            kineticEnergy,
            momentumDensity,
            velocityDistribution: {
                mean: initialVelocity,
                standardDeviation,
                skewness
            },
            effectiveGurneyCoefficient: gurneyConstant,
            casingExpansionRatio,
            accelerationGradient,
            specificImpulse,
            energyPartitionRatio
        };
    }

    // Mott distribution for fragment size distribution
    static calculateFragmentDistribution(casing: Casing, fragmentVelocity: number, totalMass: number): {
        averageFragmentMass: number; // g
        numberOfFragments: number;
        maxFragmentSize: number; // cm
        fragmentDensity: number; // fragments/m²
    } {
        // Mott's constants
        const mottConstant = Math.sqrt(2 * casing.ultimateStrength / casing.density); // m/s

        // Calculate average fragment mass
        const averageFragmentMass = Math.pow(mottConstant / fragmentVelocity, 2) * 1000; // g

        // Calculate total number of fragments
        const numberOfFragments = totalMass * 1000 / averageFragmentMass;

        // Calculate maximum fragment size (95th percentile)
        const maxFragmentSize = Math.pow(3 * averageFragmentMass / (4 * Math.PI * casing.density), 1 / 3) * 3 * 10; // cm

        // Approximate fragment density at 1m
        const fragmentDensity = numberOfFragments / (4 * Math.PI); // fragments/m² at 1m

        return { averageFragmentMass, numberOfFragments, maxFragmentSize, fragmentDensity };
    }

    // Taylor-Sedov blast wave model with full Rankine-Hugoniot relations
    static calculateBlastWaveDynamics(explosive: Explosive, yieldEnergy: number, distance: number, time: number): {
        shockRadius: number; // m
        shockVelocity: number; // m/s
        shockOverpressure: number; // kPa
        shockDensity: number; // kg/m³
        shockTemperature: number; // K
        specificInternalEnergy: number; // J/kg
        particleVelocity: number; // m/s
        machNumber: number; // dimensionless
        flowAngle: number; // radians
    } {
        // Convert yield to energy with high precision
        const energy = yieldEnergy * explosive.energyDensity * 1e6; // J

        // Atmospheric conditions
        const gamma = 1.4; // Ratio of specific heats for air
        const ambientPressure = this.ATMOSPHERIC_PRESSURE * 1000; // Pa
        const ambientTemperature = 293.15; // K (standard conditions)
        const ambientDensity = this.AIR_DENSITY; // kg/m³
        const gasConstant = 287; // J/(kg·K) for air

        // Calculate dimensionality factor (n)
        // n=3 for spherical, n=2 for cylindrical, n=1 for planar
        const n = 3; // Spherical blast wave

        // Calculate Taylor-Sedov constant with full derivation
        // Based on Sedov's self-similar solution
        let alpha: number;
        if (n === 3) {
            // Exact coefficient for strong spherical shock
            alpha = 1.175; // Derived from Taylor's solution for gamma=1.4
        } else if (n === 2) {
            alpha = 1.0; // For cylindrical shock
        } else {
            alpha = 0.851; // For planar shock
        }

        // Calculate shock radius at given time using full Sedov solution
        const E0 = energy * ((gamma - 1) / (2 * (gamma + 1))); // Effective energy parameter
        const shockRadius = alpha * Math.pow(E0 / (ambientDensity * Math.pow(time, 2)), 1 / (n + 2)) * Math.pow(time, 2 / (n + 2));

        // Calculate shock velocity directly from derivative of radius w.r.t. time
        const shockVelocity = (2 / (n + 2)) * (shockRadius / time);

        // Calculate Mach number
        const machNumber = shockVelocity / Math.sqrt(gamma * ambientPressure / ambientDensity);

        // Apply full Rankine-Hugoniot relations for shock jump conditions
        // Pressure ratio across shock
        const pressureRatio = 1 + ((2 * gamma) / (gamma + 1)) * (Math.pow(machNumber, 2) - 1);

        // Density ratio across shock
        const densityRatio = ((gamma + 1) * Math.pow(machNumber, 2)) /
            ((gamma - 1) * Math.pow(machNumber, 2) + 2);

        // Temperature ratio across shock
        const temperatureRatio = pressureRatio * (1 / densityRatio);

        // Calculate actual values
        const shockOverpressure = (pressureRatio - 1) * ambientPressure / 1000; // kPa
        const shockDensity = densityRatio * ambientDensity; // kg/m³
        const shockTemperature = temperatureRatio * ambientTemperature; // K

        // Calculate particle velocity behind shock
        const particleVelocity = (2 * (Math.pow(machNumber, 2) - 1)) /
            ((gamma + 1) * machNumber) * this.SPEED_OF_SOUND;

        // Calculate specific internal energy behind shock
        const specificInternalEnergy = shockOverpressure * 1000 /
            ((gamma - 1) * shockDensity); // J/kg

        // Calculate flow angle (relevant for non-normal incidence)
        const beta = Math.asin(1 / machNumber); // Mach angle
        const flowAngle = Math.atan(2 * Math.pow(Math.tan(beta), -1) *
            ((gamma - 1) + 2 * Math.pow(Math.sin(beta), -1)));

        return {
            shockRadius,
            shockVelocity,
            shockOverpressure,
            shockDensity,
            shockTemperature,
            specificInternalEnergy,
            particleVelocity,
            machNumber,
            flowAngle
        };
    }

    // Hopkinson-Cranz scaling law
    static scaleExplosiveEffects(referenceYield: number, referenceDistance: number, newYield: number): number {
        return referenceDistance * Math.pow(newYield / referenceYield, 1 / 3);
    }

    // Comprehensive underwater explosion model combining Cole, Geers-Hunter and Arons models
    static calculateUnderwaterExplosion(yieldEnergy: number, depth: number, distance: number): {
        peakPressure: number; // MPa
        peakPressureTime: number; // ms
        pressureHistoryParams: { A: number, θ1: number, θ2: number }; // Cole exponential parameters
        shockWaveEnergy: number; // J
        bubblePulsation: {
            maxRadius: number; // m
            firstPeriod: number; // s
            secondPeriod: number; // s
            migrationDistance: number; // m
            secondMaxRadius: number; // m
            pulseRatio: number; // dimensionless
        };
        cavitation: {
            cavitationDepth: number; // m
            cavitationRadius: number; // m
            cavitationTime: number; // ms
        };
        reflectionCoefficients: {
            freeSurface: number;
            rigidSurface: number;
        };
        impulse: number; // kPa·s
        energyFluxDensity: number; // J/m²
        lethalRadius: { fish: number, swimmer: number }; // m
    } {
        // Constants for water with temperature and salinity adjustments
        const waterTemperature = 283.15; // K (10°C)
        const salinity = 35; // ppt (parts per thousand)

        // Calculate exact water density using UNESCO equation of state
        const a0 = 999.842594;
        const a1 = 6.793952e-2;
        const a2 = -9.095290e-3;
        const a3 = 1.001685e-4;
        const a4 = -1.120083e-6;
        const a5 = 6.536332e-9;
        const t = waterTemperature - 273.15; // Convert to Celsius
        const densityFreshwater = a0 + a1 * t + a2 * t * t + a3 * t * t * t + a4 * t * t * t * t + a5 * t * t * t * t * t;

        // Salinity correction
        const b0 = 8.24493e-1;
        const b1 = -4.0899e-3;
        const b2 = 7.6438e-5;
        const b3 = -8.2467e-7;
        const b4 = 5.3875e-9;
        const c0 = -5.72466e-3;
        const c1 = 1.0227e-4;
        const c2 = -1.6546e-6;
        const salinityCorrectionTerm = (b0 + b1 * t + b2 * t * t + b3 * t * t * t + b4 * t * t * t * t) * salinity +
            (c0 + c1 * t + c2 * t * t) * Math.pow(salinity, 1.5);

        const waterDensity = densityFreshwater + salinityCorrectionTerm;

        // Calculate water sound speed with salinity and temperature effects
        const speedOfSoundInWater = 1449.2 + 4.6 * t - 0.055 * t * t + 0.00029 * t * t * t +
            (1.34 - 0.01 * t) * (salinity - 35);

        // Convert yield to energy with TNT-equivalent correction
        const energy = yieldEnergy * 4.184e9; // J (TNT equivalent)

        // Calculate hydrostatic pressure at depth with atmospheric pressure
        const gravitationalAcceleration = 9.80665; // m/s²
        const atmosphericPressure = 101325; // Pa
        const hydrostaticPressurePa = atmosphericPressure + waterDensity * gravitationalAcceleration * depth;
        const hydrostaticPressure = hydrostaticPressurePa / 1000; // kPa

        // Calculate shock wave parameters using Cole-Weston-Arons composite model
        // Direct TNT-scaling for charge mass
        const tntEquivalentMass = yieldEnergy; // kg

        // Calculate dimensionless scaled distance
        const scaledDistance = distance / Math.pow(tntEquivalentMass, 1 / 3);

        // Arons-Cole peak pressure function with empirical corrections
        // Improved fit to experimental data (Arons 1954, Cole 1948, Hunter 1970)
        const K1 = 52.4; // MPa·m/kg^(1/3)
        const alpha = 1.13; // Cole's pressure decay constant

        // Calculate pressure decay time constants for realistic exponential-power model
        const θ1 = 0.056 * Math.pow(tntEquivalentMass, 1 / 3) * Math.pow(scaledDistance, 0.05); // ms
        const θ2 = 0.092 * Math.pow(tntEquivalentMass, 1 / 3) * Math.pow(scaledDistance, 0.18); // ms

        // Calculate peak pressure with corrected model
        const peakPressure = K1 * Math.pow(tntEquivalentMass, 0.63 / 3) *
            Math.pow(scaledDistance, -alpha) *
            (1 + 0.03 * Math.log10(tntEquivalentMass)); // MPa

        // Time to peak pressure
        const peakPressureTime = distance / speedOfSoundInWater * 1000; // ms

        // Energy flux density (Arons, 1954)
        const energyFluxDensity = 0.35 * peakPressure * 1e6 * θ1 * 1e-3; // J/m²

        // Calculate impulse using integrated pressure history
        const impulse = peakPressure * (θ1 + θ2) * 1e-3; // kPa·s

        // Calculate shock wave energy using Willis formula
        const shockWaveEnergy = (energy * 0.53) * Math.exp(-0.22 * depth / Math.pow(tntEquivalentMass, 1 / 3));

        // Calculate reflection coefficients
        const freeSurfaceReflection = -1.0; // Perfect pressure release
        // Acoustic impedance method for rigid surface
        const rigidSurfaceReflection = 0.95; // Slightly less than theoretical due to energy absorption

        // Calculate bubble parameters using exact Rayleigh-Plesset equation parameters
        // Initial condition parameters from Geers-Hunter model
        const k1Bubble = 3.5; // First bubble coefficient, revised based on experimental data
        const k2Bubble = 0.87; // Second bubble coefficient

        // Max bubble radius 
        const maxRadius = k1Bubble * Math.pow(tntEquivalentMass / (hydrostaticPressurePa / 1e5), 1 / 3); // m

        // Calculate bubble period including compressibility correction
        const rayleighPeriod = 2.11 * maxRadius * Math.sqrt(waterDensity / hydrostaticPressurePa);
        const compressibilityFactor = 1 + 0.02 * depth / maxRadius;
        const firstPeriod = rayleighPeriod * compressibilityFactor; // s

        // Calculate bubble migration velocity from buoyancy
        const bubbleBuoyancy = (4 / 3) * Math.PI * Math.pow(maxRadius, 3) * waterDensity * gravitationalAcceleration;
        const addedMassFactor = 0.5; // Hydrodynamic added mass
        const effectiveMass = (4 / 3) * Math.PI * Math.pow(maxRadius, 3) * waterDensity * (1 + addedMassFactor);
        const migrationVelocity = bubbleBuoyancy / effectiveMass * firstPeriod; // m/s

        // Calculate migration distance during first oscillation
        const migrationDistance = 0.5 * migrationVelocity * firstPeriod; // m

        // Calculate second maximum radius (smaller due to energy loss)
        const energyLossFactor = 0.85; // Energy loss during first oscillation
        const secondMaxRadius = maxRadius * Math.sqrt(energyLossFactor); // m

        // Calculate second oscillation period
        const secondPeriod = firstPeriod * Math.sqrt(energyLossFactor); // s

        // Calculate pulse pressure ratio between first and second bubble pulse
        const pulseRatio = 0.54 * Math.pow(tntEquivalentMass, -0.05); // dimensionless

        // Calculate cavitation parameters
        // Critical pressure for cavitation
        const cavitationPressure = -101325; // Pa (-1 atm)

        // Calculate maximum depth where cavitation can occur
        const cavitationDepth = 10.4 * Math.pow(tntEquivalentMass, 1 / 3) * Math.pow(peakPressure / 5, 0.4); // m

        // Calculate maximum radius of cavitation zone
        const cavitationRadius = 15 * Math.pow(tntEquivalentMass, 1 / 3) * Math.sqrt(depth / Math.max(cavitationDepth, 0.1)); // m

        // Calculate time of cavitation closure
        const cavitationTime = 5.5 * Math.pow(tntEquivalentMass, 1 / 3) * Math.sqrt(depth / Math.max(cavitationDepth, 0.1)); // ms

        // Calculate lethal radius for biological targets
        // Based on empirical studies (Goertner model)
        const lethalRadiusFish = 43 * Math.pow(tntEquivalentMass, 1 / 3) * Math.pow(depth / 10, -0.22); // m
        const lethalRadiusSwimmer = 12 * Math.pow(tntEquivalentMass, 1 / 3) * Math.pow(depth / 10, -0.27); // m

        return {
            peakPressure,
            peakPressureTime,
            pressureHistoryParams: { A: peakPressure, θ1, θ2 },
            shockWaveEnergy,
            bubblePulsation: {
                maxRadius,
                firstPeriod,
                secondPeriod,
                migrationDistance,
                secondMaxRadius,
                pulseRatio
            },
            cavitation: {
                cavitationDepth,
                cavitationRadius,
                cavitationTime
            },
            reflectionCoefficients: {
                freeSurface: freeSurfaceReflection,
                rigidSurface: rigidSurfaceReflection
            },
            impulse,
            energyFluxDensity,
            lethalRadius: {
                fish: lethalRadiusFish,
                swimmer: lethalRadiusSwimmer
            }
        };
    }

    // Ground shock calculation
    static calculateGroundShock(energyYield: number, distance: number, soilProperties: {
        pWaveVelocity: number, // m/s
        sWaveVelocity: number, // m/s
        density: number, // kg/m³
        attenuation: number // dimensionless
    }): {
        peakParticleVelocity: number, // m/s
        peakAcceleration: number, // g
        peakDisplacement: number, // mm
        arrivalTime: number // s
    } {
        const scaledDistanceResult = this.calculateScaledDistance(distance, energyYield);
        const scaledDistance = scaledDistanceResult.hopkinsonScaledDistance;

        // Calculate peak particle velocity (empirical formula)
        const peakParticleVelocity = 160 * Math.pow(scaledDistance, -1.5) * Math.pow(soilProperties.attenuation, -distance / 100);

        // Calculate peak acceleration
        const peakAcceleration = peakParticleVelocity * soilProperties.pWaveVelocity / 9.81; // in g

        // Calculate peak displacement
        const dominantFrequency = 10 * Math.pow(energyYield, -1 / 3); // Hz (empirical)
        const peakDisplacement = peakParticleVelocity / (2 * Math.PI * dominantFrequency) * 1000; // mm

        // Calculate arrival time
        const arrivalTime = distance / soilProperties.pWaveVelocity; // s

        return { peakParticleVelocity, peakAcceleration, peakDisplacement, arrivalTime };
    }

    // Advanced shaped charge jet formation model combining Birkhoff, PER, and Chou-Carleone models
    static calculateShapedChargeJet(explosive: Explosive, liner: Liner, standoff: number): {
        // Basic jet parameters
        jetVelocity: { tip: number, tail: number, average: number }, // m/s
        penetrationDepth: number, // cm
        jetLength: number, // cm
        timeOfFlight: number, // μs
        jetMass: number, // g
        slugMass: number, // g

        // Advanced jet parameters
        jetCoherence: number, // dimensionless (0-1)
        breakupTime: number, // μs
        breakupDistance: number, // cm
        stretching: number, // 1/μs
        particulation: {
            onsetTime: number, // μs
            particleDiameter: number, // mm
            particleCount: number, // dimensionless
            particleVelocitySpread: number // m/s
        },

        // Target interaction
        penetrationVelocity: number, // km/s
        holeProfile: {
            entranceDiameter: number, // mm
            exitDiameter: number, // mm
            averageDiameter: number // mm
        },
        afterEffect: {
            behindArmorEffect: number, // cm³
            spallDiameter: number, // cm
            spallVelocity: number // m/s
        },

        // Jet formation dynamics
        formationTime: number, // μs
        detonationPressure: number, // GPa
        collapseVelocity: number, // km/s
        stagnationPressure: number, // GPa
        velocityGradient: number, // 1/μs
        jettingAnalysis: {
            gurney1D: number, // m/s
            gurney2D: number, // m/s
            PDEfactor: number, // dimensionless
            massEfficiency: number // %
        }
    } {
        // Constants
        const PI = Math.PI;

        // Calculate detonation pressure using actual JWL or Chapman-Jouguet equation
        // P = ρ * D² / (γ + 1) where γ is the polytropic exponent (approximated as 3 for most explosives)
        const polytropic = 3;
        const detonationPressure = explosive.density * Math.pow(explosive.detonationVelocity, 2) / 1e6 / (polytropic + 1); // GPa

        // Calculate liner aspect ratio (important for collapse dynamics)
        const aspectRatio = 2 * (liner.thickness / (liner.material === 'copper' ? 1.5 : 1.0)); // normalized thickness

        // Calculate collapse velocity using PER (Progressive Explosive Reaction) model
        // More accurate than simple Birkhoff approximation
        // v_c = 0.6 * v_d * (1 - e^(-2M)) where M is the explosive to liner mass ratio
        const explosiveToLinerRatio = 3.5; // typical value, would be calculated from geometry
        const collapseFactor = 0.6 * (1 - Math.exp(-2 * explosiveToLinerRatio));
        const collapseVelocity = collapseFactor * explosive.detonationVelocity / 1000; // km/s

        // Calculate jet and slug velocities using Birkhoff model with correction factors
        // Corrections account for liner thickness effects (Walsh model)
        const alpha = liner.angle * PI / 180; // convert to radians
        const beta = alpha / 2; // beta angle in radians

        // Correction factor for finite liner thickness
        const thicknessCorrection = 1 - 0.1 * Math.pow(liner.thickness, 0.3);

        // Calculate stagnation angle using Carleone formula
        const tan_alpha = Math.tan(alpha);
        const stagnationAngle = Math.atan(0.5 * tan_alpha); // radians

        // Calculate stagnation point pressure (Chou-Flis model)
        const stagnationPressure = detonationPressure * Math.pow(Math.sin(alpha), 2) *
            (1 + 2 * Math.cos(2 * stagnationAngle)); // GPa

        // Calculate jet and slug velocities using Chou-Carleone model
        // Accounts for stagnation point dynamics and compressibility
        const jetVelocityTip = collapseVelocity * (1 + Math.cos(beta)) * thicknessCorrection * 1000; // m/s
        const jetVelocityTail = collapseVelocity * (1 + Math.cos(alpha)) * thicknessCorrection * 1000; // m/s
        const slugVelocity = collapseVelocity * (1 - Math.cos(alpha)) * 1000; // m/s

        // Calculate average jet velocity
        const jetVelocityAvg = (jetVelocityTip + jetVelocityTail) / 2; // m/s

        // Calculate velocity gradient along jet
        const velocityGradient = (jetVelocityTip - jetVelocityTail) / 100; // 1/μs (assuming typical 100μs formation)

        // Calculate mass partition using Chou model
        // λ = 2*sin²(β) / sin²(α)
        const lambda = 2 * Math.pow(Math.sin(beta), 2) / Math.pow(Math.sin(alpha), 2);

        // Calculate total liner mass (simplified conical section)
        const linerVolume = PI * Math.pow(liner.thickness / 10, 3) * (1 / 3) * (1 / Math.tan(alpha)); // cm³
        const totalLinerMass = linerVolume * liner.density; // g

        // Calculate jet and slug masses
        const jetMass = lambda * totalLinerMass / (1 + lambda); // g
        const slugMass = totalLinerMass - jetMass; // g

        // Calculate formation time based on liner dimensions and detonation velocity
        const linerDiameter = liner.thickness / Math.tan(alpha) * 2; // mm
        const formationTime = (linerDiameter / 2) / explosive.detonationVelocity * 1e6; // μs

        // Calculate jet coherence parameter (Pugh-Eichelberger-Rostoker model)
        // Function of liner material, thickness, and angle
        // 1.0 is perfect coherence, 0.0 is immediate breakup
        const materialFactor = liner.material === 'copper' ? 1.0 :
            liner.material === 'aluminum' ? 0.85 :
                liner.material === 'tantalum' ? 0.95 : 0.9;

        const angleCorrection = Math.cos(alpha - 0.6);
        const jetCoherence = materialFactor * angleCorrection * Math.exp(-0.5 * liner.thickness / 10);

        // Calculate jet breakup time using modified Hirsch model
        // t_b = J * d / (Vj_tip - Vj_tail) where J is material-dependent constant
        const hirschConstant = liner.material === 'copper' ? 2.0 :
            liner.material === 'aluminum' ? 1.4 :
                liner.material === 'tantalum' ? 2.8 : 1.8;

        const breakupTime = hirschConstant * liner.thickness /
            (jetVelocityTip - jetVelocityTail) * 1e6; // μs

        // Calculate stretching rate
        const stretching = (jetVelocityTip - jetVelocityTail) / (jetVelocityAvg * breakupTime); // 1/μs

        // Calculate jet length at breakup
        const jetLength = jetVelocityAvg * breakupTime / 10000; // cm

        // Calculate breakup distance
        const breakupDistance = standoff * (breakupTime / (standoff * 10 / jetVelocityAvg * 1e6)); // cm

        // Calculate particulation parameters
        // Based on instability theory (Rayleigh-Taylor and Richtmyer-Meshkov)
        const onsetTime = breakupTime * 1.2; // μs

        // Calculate particle diameter using Walsh model
        const particleDiameter = 1.81 * liner.thickness /
            Math.pow(1 + 0.1 * velocityGradient, 0.5); // mm

        // Calculate number of particles
        const particleCount = Math.ceil(jetLength * 10 / particleDiameter);

        // Calculate particle velocity spread due to breakup
        const particleVelocitySpread = 0.05 * jetVelocityAvg; // m/s

        // Calculate time of flight to standoff
        const timeOfFlight = standoff * 10 / jetVelocityAvg * 1e6; // μs

        // Calculate penetration velocity
        // Using hydrodynamic penetration theory
        const targetDensity = 7.85; // steel target density g/cm³
        const penetrationVelocity = jetVelocityAvg * Math.sqrt(liner.density / targetDensity) / 1000; // km/s

        // Calculate penetration depth using modified Tate equation with shear strength terms
        // P = L * √(ρj/ρt) * [1 - Rt/(ρj*Vj²/2) - Yj/(ρj*Vj²/2)]
        // where L is jet length, R is target strength, Y is jet strength
        const jetStrength = liner.material === 'copper' ? 0.3 :
            liner.material === 'aluminum' ? 0.2 :
                liner.material === 'tantalum' ? 0.7 : 0.4; // GPa

        const targetStrength = 1.2; // GPa for steel

        // Jet kinetic energy density
        const jetKE = 0.5 * liner.density * Math.pow(jetVelocityAvg, 2) / 1e9; // GPa

        // Calculate penetration efficiency
        const penetrationEfficiency = 1 - targetStrength / jetKE - jetStrength / jetKE;

        // Calculate penetration depth
        const penetrationDepth = jetLength * Math.sqrt(liner.density / targetDensity) *
            Math.max(0.1, penetrationEfficiency) *
            (timeOfFlight < breakupTime ? 1 : Math.pow(breakupTime / timeOfFlight, 0.3)); // cm

        // Calculate hole profile
        // Entrance diameter is typically 3-5 times the jet diameter
        const jetDiameter = 0.5 * liner.thickness * Math.sqrt(lambda); // mm
        const entranceDiameter = 4 * jetDiameter; // mm

        // Exit diameter is typically 1-2 times the jet diameter
        const exitDiameter = Math.max(1, (entranceDiameter * 0.4) - (penetrationDepth * 0.05)); // mm

        // Average hole diameter
        const averageDiameter = (entranceDiameter + exitDiameter) / 2; // mm

        // Calculate behind armor effects
        // Behind armor debris volume
        const behindArmorEffect = PI * Math.pow(exitDiameter / 10, 2) * penetrationDepth / 4; // cm³

        // Spall diameter using Christman-Gehring formula
        const spallDiameter = exitDiameter * (1 + 5 * Math.pow(penetrationVelocity, 0.6)) / 10; // cm

        // Spall velocity using momentum conservation
        const spallMass = PI * Math.pow(spallDiameter, 2) * targetDensity * 0.1; // g
        const spallVelocity = jetMass * jetVelocityAvg / spallMass * 0.3; // m/s

        // Calculate Gurney velocities for comparison
        // 1D Gurney (sandwich)
        const gurney1D = 0.5 * Math.sqrt(2 * explosive.energyDensity * 1e6) *
            Math.sqrt(explosiveToLinerRatio / (1 + explosiveToLinerRatio / 2)); // m/s

        // 2D Gurney (cylinder)
        const gurney2D = 0.5 * Math.sqrt(2 * explosive.energyDensity * 1e6) *
            Math.sqrt(explosiveToLinerRatio / (1 + explosiveToLinerRatio / 3)); // m/s

        // PDE (Prompt Detonation Energy) factor - ratio of actual to ideal detonation
        const PDEfactor = detonationPressure /
            (explosive.density * Math.pow(explosive.detonationVelocity, 2) / 1e6 / 4); // dimensionless

        // Calculate mass efficiency - percentage of liner mass converted to useful jet
        const massEfficiency = (jetMass / totalLinerMass) * 100; // %

        return {
            jetVelocity: {
                tip: jetVelocityTip,
                tail: jetVelocityTail,
                average: jetVelocityAvg
            },
            penetrationDepth,
            jetLength,
            timeOfFlight,
            jetMass,
            slugMass,

            jetCoherence,
            breakupTime,
            breakupDistance,
            stretching,
            particulation: {
                onsetTime,
                particleDiameter,
                particleCount,
                particleVelocitySpread
            },

            penetrationVelocity,
            holeProfile: {
                entranceDiameter,
                exitDiameter,
                averageDiameter
            },
            afterEffect: {
                behindArmorEffect,
                spallDiameter,
                spallVelocity
            },

            formationTime,
            detonationPressure,
            collapseVelocity,
            stagnationPressure,
            velocityGradient,
            jettingAnalysis: {
                gurney1D,
                gurney2D,
                PDEfactor,
                massEfficiency
            }
        };
    }

    // Thermal radiation effects
    static calculateThermalRadiation(energyYield: number, distance: number, transmissivity: number = 0.9): {
        thermalEnergy: number, // J/cm²
        fireball: { radius: number, duration: number, temperature: number } // m, s, K
    } {
        // Convert energyYield to energy
        const energy = energyYield * 4.184e12; // J (TNT equivalent)

        // Calculate fireball radius (Glasstone & Dolan)
        const fireball = {
            radius: 3.7 * Math.pow(energyYield, 0.4), // m
            duration: 0.9 * Math.pow(energyYield, 0.4), // s
            temperature: 7000 // K (initial temperature)
        };

        // Thermal partition (approximate)
        const thermalFraction = 0.35; // approx 35% of energy as thermal radiation

        // Calculate thermal energy at distance
        const thermalEnergy = thermalFraction * energy * transmissivity / (4 * Math.PI * Math.pow(distance, 2)) / 100; // J/cm²

        return { thermalEnergy, fireball };
    }

    // Calculate ionizing radiation dose
    static calculateRadiationDose(energyYield: number, distance: number): {
        initialNuclearRadiation: number, // Gy
        residualRadiation: number // Gy/hr at 1 hour
    } {
        // This applies only to nuclear explosions - return zero for conventional
        return {
            initialNuclearRadiation: 0,
            residualRadiation: 0
        };
    }

    // Improved structural damage model
    static calculateStructuralDamage(structure: {
        type: string,
        peakOverpressureThreshold: number, // kPa
        impulseThreshold: number, // kPa·ms
        pDelta: number // pressure-impulse curve exponent
    }, peakOverpressure: number, impulse: number): {
        damagePercentage: number, // %
        damageMode: string
    } {
        // Pressure-Impulse diagram approach (simplified)
        const pressureRatio = peakOverpressure / structure.peakOverpressureThreshold;
        const impulseRatio = impulse / structure.impulseThreshold;

        // Calculate damage using P-I curve
        const piValue = Math.pow(pressureRatio, structure.pDelta) + Math.pow(impulseRatio, structure.pDelta);

        // Determine damage percentage
        let damagePercentage = 0;
        if (piValue >= 1) {
            damagePercentage = Math.min(100, 100 * Math.pow(piValue, 1 / 2));
        }

        // Determine damage mode
        let damageMode = "None";
        if (damagePercentage > 0) {
            if (pressureRatio > impulseRatio) {
                damageMode = "Pressure-dominated";
            } else {
                damageMode = "Impulse-dominated";
            }
        }

        return { damagePercentage, damageMode };
    }

    // Human injury model using Bowen criteria
    static calculateHumanInjury(peakOverpressure: number, positivePhaseDuration: number): {
        lethality: number, // % probability
        lungDamage: string,
        eardrumRupture: number // % probability
    } {
        // Convert to psi for standard Bowen curves
        const overpressure_psi = peakOverpressure * 0.145038; // kPa to psi

        // Scale duration
        const scaledDuration = positivePhaseDuration * Math.pow(70 / 155, 1 / 3); // scaled to 70kg human

        // Calculate lung damage probability (simplified Bowen model)
        let lethality = 0;
        if (overpressure_psi > 50) {
            lethality = 100;
        } else if (overpressure_psi > 30) {
            lethality = 50 + 50 * (overpressure_psi - 30) / 20;
        } else if (overpressure_psi > 15) {
            lethality = (overpressure_psi - 15) / 15 * 50;
        }

        // Determine lung damage severity
        let lungDamage = "None";
        if (overpressure_psi > 12) {
            lungDamage = "Severe";
        } else if (overpressure_psi > 6) {
            lungDamage = "Moderate";
        } else if (overpressure_psi > 3) {
            lungDamage = "Threshold";
        }

        // Calculate eardrum rupture probability
        const eardrumRupture = Math.min(100, Math.max(0, -12 + 19 * Math.log(overpressure_psi)));

        return { lethality, lungDamage, eardrumRupture };
    }

    // Calculate mass of explosive required for specific effect
    static calculateRequiredMass(explosive: Explosive, targetEffect: {
        type: 'overpressure' | 'impulse' | 'fragmentation' | 'cratering' | 'penetration',
        value: number,
        distance: number
    }): number {
        let requiredYield = 0;

        if (targetEffect.type === 'overpressure') {
            // Convert target overpressure to scaled distance using inverse Kingery-Bulmash
            const z = Math.pow(targetEffect.value / 98, -0.6);
            // Calculate required yield
            requiredYield = Math.pow(targetEffect.distance / z, 3);
        } else if (targetEffect.type === 'impulse') {
            // Similar approach for impulse
            const z = Math.pow(targetEffect.value / 250, -0.6);
            requiredYield = Math.pow(targetEffect.distance / z, 3);
        } else if (targetEffect.type === 'cratering') {
            // For crater radius
            requiredYield = Math.pow(targetEffect.value / 0.8, 3);
        }

        // Convert from TNT equivalent to actual explosive mass
        return requiredYield / (explosive.energyDensity / 4.6);
    }

    // Calculate optimal standoff for shaped charge
    static calculateOptimalStandoff(liner: Liner): number {
        // Typical optimum standoff is 3-7 cone diameters
        return 5 * liner.thickness / Math.tan(liner.angle * Math.PI / 180);
    }
}

// Enhanced Implementation of optimizeExplosive function
function optimizeExplosive(
    explosives: Explosive[],
    scenario: {
        radius: number,
        shape: Shape,
        maxWeight: number,
        isPortable: boolean,
        casing: Casing,
        environment: Environment,
        targets?: Target[]
    }
): { material: Explosive; metrics: any } | null {
    // Filter explosives based on portability and stability
    const filteredExplosives = explosives.filter(explosive => {
        if (scenario.isPortable && explosive.stability < 3) return false;

        // Temperature sensitivity check
        if (scenario.environment.temperature < 223 || scenario.environment.temperature > 343) {
            if (explosive.stability < 4) return false;
        }

        // Friction sensitivity check for portable devices
        if (scenario.isPortable && explosive.frictionSensitivity < 80) return false;

        return true;
    });

    // If no suitable explosives remain, return null
    if (filteredExplosives.length === 0) return null;

    // Calculate required mass for desired effect at specified radius
    let bestMaterial: Explosive | null = null;
    let bestMetrics: any = null;
    let bestFigureOfMerit = 0;

    for (const explosive of filteredExplosives) {
        // Calculate explosive yield in TNT equivalent
        const tntEquivalent = explosive.energyDensity / 4.6;

        // Calculate scaled distance
        const scaledDistance = ExplosiveModel.calculateScaledDistance(scenario.radius, 1);

        // Calculate blast parameters
        const blastParams = ExplosiveModel.calculateAirblastParameters(scaledDistance.hopkinsonScaledDistance);

        // Calculate required mass to achieve desired overpressure
        const requiredMass = ExplosiveModel.calculateRequiredMass(explosive, {
            type: 'overpressure',
            value: blastParams.overpressure,
            distance: scenario.radius
        });

        // Check if mass fits within weight constraint
        let casingMass = 0;
        if (scenario.shape.type === 'sphere') {
            const volume = (4 / 3) * Math.PI * Math.pow(scenario.shape.dimensions.radius!, 3);
            const surface = 4 * Math.PI * Math.pow(scenario.shape.dimensions.radius!, 2);
            casingMass = (surface * scenario.casing.thickness * scenario.casing.density) / 1000;
        } else if (scenario.shape.type === 'cylinder') {
            const r = scenario.shape.dimensions.radius!;
            const h = scenario.shape.dimensions.height!;
            const surface = 2 * Math.PI * r * (r + h);
            casingMass = (surface * scenario.casing.thickness * scenario.casing.density) / 1000;
        } else if (scenario.shape.type === 'shaped_charge' && scenario.shape.linerProperties) {
            // For shaped charges, consider liner mass
            const cone = scenario.shape.dimensions.radius! * scenario.shape.dimensions.height! * Math.PI / 3;
            casingMass = (cone * scenario.casing.thickness * scenario.casing.density) / 1000;
            casingMass += (cone * scenario.shape.linerProperties.thickness * scenario.shape.linerProperties.density) / 1000;
        }

        const totalMass = requiredMass + casingMass;

        if (totalMass > scenario.maxWeight) continue;

        // Calculate fragmentation parameters
        const fragVelocity = ExplosiveModel.calculateFragmentVelocity(explosive, scenario.casing, {
            mass: requiredMass,
            shape: scenario.shape.type
        });

        const fragDistribution =
            ExplosiveModel.calculateFragmentDistribution(
                scenario.casing,
                fragVelocity.initialVelocity,
                casingMass,
            );

        // Calculate thermal effects
        const thermalEffects = ExplosiveModel.calculateThermalRadiation(requiredMass * tntEquivalent, scenario.radius);

        // Calculate ground shock if applicable
        let groundShock: any = null;
        if (scenario.environment.soilType && scenario.environment.soilDensity) {
            groundShock = ExplosiveModel.calculateGroundShock(requiredMass * tntEquivalent, scenario.radius, {
                pWaveVelocity: 500, // m/s (typical soil)
                sWaveVelocity: 300, // m/s (typical soil)
                density: scenario.environment.soilDensity * 1000,
                attenuation: 0.95
            });
        }

        // Calculate underwater effects if applicable
        let underwaterEffects: any = null;
        if (scenario.environment.waterDepth) {
            underwaterEffects = ExplosiveModel.calculateUnderwaterExplosion(
                requiredMass * tntEquivalent,
                scenario.environment.waterDepth,
                scenario.radius
            );
        }

        // Calculate shaped charge effects if applicable
        let shapedChargeEffects: any = null;
        if (scenario.shape.type === 'shaped_charge' && scenario.shape.linerProperties) {
            shapedChargeEffects = ExplosiveModel.calculateShapedChargeJet(
                explosive,
                scenario.shape.linerProperties,
                scenario.shape.dimensions.standoff || 0
            );
        }

        // Calculate human injury and structural damage if targets specified
        let targetsAssessment: any = null;
        if (scenario.targets && scenario.targets.length > 0) {
            targetsAssessment = scenario.targets.map(target => {
                const targetScaledDistance = ExplosiveModel.calculateScaledDistance(
                    target.distance,
                    requiredMass * tntEquivalent
                );
                const targetBlastParams = ExplosiveModel.calculateAirblastParameters(targetScaledDistance.hopkinsonScaledDistance);

                // Calculate injury for human targets
                const injury = target.material === 'human' ?
                    ExplosiveModel.calculateHumanInjury(
                        targetBlastParams.overpressure,
                        targetBlastParams.positivePhaseDuration
                    ) : null;

                // Calculate structural damage
                const damage = ExplosiveModel.calculateStructuralDamage(
                    {
                        type: target.material,
                        peakOverpressureThreshold: target.criticalDamageThreshold,
                        impulseThreshold: target.criticalDamageThreshold * 10, // approximate
                        pDelta: 1.5 // typical value
                    },
                    targetBlastParams.overpressure,
                    targetBlastParams.positivePhaseImpulse
                );

                return { target, blastParams: targetBlastParams, injury, damage };
            });
        }

        // Calculate figure of merit (multi-criteria optimization)
        // Higher is better
        const figureOfMerit =
            (1 / requiredMass) * 5 + // less mass is better
            (explosive.stability) * 2 + // more stability is better
            (explosive.detonationVelocity / 10000) * 3; // higher velocity is better

        // Update best material if this has a better figure of merit
        if (figureOfMerit > bestFigureOfMerit) {
            bestFigureOfMerit = figureOfMerit;
            bestMaterial = explosive;

            // Compile comprehensive metrics
            bestMetrics = {
                // Basic properties
                stability: explosive.stability,
                detonationVelocity: explosive.detonationVelocity,
                detonationPressure: explosive.chapmanJouguetPressure,
                detonationTemperature: explosive.gurvichTemperature,

                // Mass and energy
                tntEquivalent: tntEquivalent,
                explosiveMass: requiredMass,
                casingMass: casingMass,
                totalWeight: requiredMass + casingMass,
                explosionEnergy: requiredMass * explosive.energyDensity,

                // Blast effects
                overpressure: blastParams.overpressure,
                positivePhaseDuration: blastParams.positivePhaseDuration,
                positivePhaseImpulse: blastParams.positivePhaseImpulse,
                reflectedPressure: blastParams.reflectedPressure,

                // Fragmentation
                fragmentVelocity: fragVelocity,
                averageFragmentMass: fragDistribution.averageFragmentMass,
                numberOfFragments: fragDistribution.numberOfFragments,
                maxFragmentSize: fragDistribution.maxFragmentSize,

                // Thermal effects
                thermalRadiation: thermalEffects.thermalEnergy,
                fireballRadius: thermalEffects.fireball.radius,
                fireballDuration: thermalEffects.fireball.duration,

                // Special effects
                groundShock,
                underwaterEffects,
                shapedChargeEffects,

                // Target effects
                targetsAssessment
            };
        }
    }

    if (!bestMaterial) return null;

    return { material: bestMaterial, metrics: bestMetrics };
}

const explosives: Explosive[] = [
    {
        "name": "TNT",
        "detonationVelocity": 6900,
        "energyDensity": 4.6,
        "density": 1.65,
        "stability": 4,
        "criticalDiameter": 5.0,
        "activationEnergy": 197,
        "gurvichTemperature": 3300,
        "chapmanJouguetPressure": 19,
        "oxygenBalance": -74,
        "thermalConductivity": 0.22,
        "specificHeat": 1230,
        "thermalExpansion": 7.9e-5,
        "shockSensitivity": 190,
        "frictionSensitivity": 353,
        "impactSensitivity": 15,
        "expansionRatio": 690,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.34, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.23, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.32, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.11, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "HMX",
        "detonationVelocity": 9100,
        "energyDensity": 5.7,
        "density": 1.91,
        "stability": 3,
        "criticalDiameter": 0.8,
        "activationEnergy": 220,
        "gurvichTemperature": 4400,
        "chapmanJouguetPressure": 39,
        "oxygenBalance": -21.6,
        "thermalConductivity": 0.3,
        "specificHeat": 1250,
        "thermalExpansion": 6.6e-5,
        "shockSensitivity": 70,
        "frictionSensitivity": 120,
        "impactSensitivity": 7.4,
        "expansionRatio": 780,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.20, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.37, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.25, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.18, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "PETN",
        "detonationVelocity": 8400,
        "energyDensity": 5.8,
        "density": 1.77,
        "stability": 2,
        "criticalDiameter": 1.0,
        "activationEnergy": 197,
        "gurvichTemperature": 4230,
        "chapmanJouguetPressure": 33.5,
        "oxygenBalance": -10.1,
        "thermalConductivity": 0.25,
        "specificHeat": 1260,
        "thermalExpansion": 8.1e-5,
        "shockSensitivity": 60,
        "frictionSensitivity": 60,
        "impactSensitivity": 3,
        "expansionRatio": 750,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.17, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.42, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.24, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.12, "heatCapacity": 29.1 },
            { "formula": "O2", "moleFraction": 0.05, "heatCapacity": 29.4 }
        ]
    },
    {
        "name": "RDX",
        "detonationVelocity": 8750,
        "energyDensity": 5.13,
        "density": 1.82,
        "stability": 3,
        "criticalDiameter": 1.0,
        "activationEnergy": 197,
        "gurvichTemperature": 4100,
        "chapmanJouguetPressure": 33.8,
        "oxygenBalance": -21.6,
        "thermalConductivity": 0.29,
        "specificHeat": 1256,
        "thermalExpansion": 7.7e-5,
        "shockSensitivity": 80,
        "frictionSensitivity": 120,
        "impactSensitivity": 7.5,
        "expansionRatio": 740,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.21, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.36, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.26, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.17, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "Comp B",
        "detonationVelocity": 7900,
        "energyDensity": 5.19,
        "density": 1.72,
        "stability": 4,
        "criticalDiameter": 1.5,
        "activationEnergy": 197,
        "gurvichTemperature": 3800,
        "chapmanJouguetPressure": 29.5,
        "oxygenBalance": -56,
        "thermalConductivity": 0.28,
        "specificHeat": 1240,
        "thermalExpansion": 7.8e-5,
        "shockSensitivity": 85,
        "frictionSensitivity": 140,
        "impactSensitivity": 8,
        "expansionRatio": 720,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.27, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.32, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.27, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.14, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "C4",
        "detonationVelocity": 8040,
        "energyDensity": 5.3,
        "density": 1.60,
        "stability": 5,
        "criticalDiameter": 1.5,
        "activationEnergy": 197,
        "gurvichTemperature": 3900,
        "chapmanJouguetPressure": 26,
        "oxygenBalance": -74,
        "thermalConductivity": 0.25,
        "specificHeat": 1230,
        "thermalExpansion": 7.9e-5,
        "shockSensitivity": 190,
        "frictionSensitivity": 353,
        "impactSensitivity": 15,
        "expansionRatio": 700,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.22, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.35, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.25, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.18, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "Semtex",
        "detonationVelocity": 7500,
        "energyDensity": 5.0,
        "density": 1.40,
        "stability": 4,
        "criticalDiameter": 1.0,
        "activationEnergy": 197,
        "gurvichTemperature": 3700,
        "chapmanJouguetPressure": 24,
        "oxygenBalance": -56,
        "thermalConductivity": 0.28,
        "specificHeat": 1240,
        "thermalExpansion": 7.8e-5,
        "shockSensitivity": 85,
        "frictionSensitivity": 140,
        "impactSensitivity": 8,
        "expansionRatio": 720,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.20, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.34, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.26, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.20, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "Dynamite",
        "detonationVelocity": 6000,
        "energyDensity": 4.0,
        "density": 1.50,
        "stability": 2,
        "criticalDiameter": 2.5,
        "activationEnergy": 150,
        "gurvichTemperature": 3200,
        "chapmanJouguetPressure": 15,
        "oxygenBalance": -10,
        "thermalConductivity": 0.20,
        "specificHeat": 1200,
        "thermalExpansion": 8.0e-5,
        "shockSensitivity": 50,
        "frictionSensitivity": 50,
        "impactSensitivity": 2,
        "expansionRatio": 650,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.25, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.30, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.30, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.15, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "ANFO",
        "detonationVelocity": 4500,
        "energyDensity": 3.7,
        "density": 0.80,
        "stability": 4,
        "criticalDiameter": 10.0,
        "activationEnergy": 100,
        "gurvichTemperature": 2800,
        "chapmanJouguetPressure": 5,
        "oxygenBalance": 10,
        "thermalConductivity": 0.15,
        "specificHeat": 1150,
        "thermalExpansion": 9.0e-5,
        "shockSensitivity": 100,
        "frictionSensitivity": 200,
        "impactSensitivity": 10,
        "expansionRatio": 600,
        "detonationProducts": [
            { "formula": "CO2", "moleFraction": 0.40, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.45, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.15, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "Nitroglycerin",
        "detonationVelocity": 7700,
        "energyDensity": 6.4,
        "density": 1.60,
        "stability": 1,
        "criticalDiameter": 0.5,
        "activationEnergy": 150,
        "gurvichTemperature": 4500,
        "chapmanJouguetPressure": 25,
        "oxygenBalance": 3.5,
        "thermalConductivity": 0.22,
        "specificHeat": 1230,
        "thermalExpansion": 7.9e-5,
        "shockSensitivity": 190,
        "frictionSensitivity": 353,
        "impactSensitivity": 15,
        "expansionRatio": 690,
        "detonationProducts": [
            { "formula": "CO2", "moleFraction": 0.50, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.30, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.15, "heatCapacity": 29.1 },
            { "formula": "O2", "moleFraction": 0.05, "heatCapacity": 29.4 }
        ]
    },
    {
        "name": "Black Powder",
        "detonationVelocity": null,
        "energyDensity": 3.0,
        "density": 1.70,
        "stability": 2,
        "criticalDiameter": null,
        "activationEnergy": 100,
        "gurvichTemperature": 2500,
        "chapmanJouguetPressure": null,
        "oxygenBalance": -40,
        "thermalConductivity": 0.10,
        "specificHeat": 1000,
        "thermalExpansion": 5.0e-5,
        "shockSensitivity": 50,
        "frictionSensitivity": 50,
        "impactSensitivity": 2,
        "expansionRatio": 500,
        "detonationProducts": [
            { "formula": "CO2", "moleFraction": 0.30, "heatCapacity": 37.1 },
            { "formula": "CO", "moleFraction": 0.20, "heatCapacity": 29.1 },
            { "formula": "N2", "moleFraction": 0.20, "heatCapacity": 29.1 },
            { "formula": "SO2", "moleFraction": 0.30, "heatCapacity": 39.9 }
        ]
    },
    {
        "name": "Smokeless Powder",
        "detonationVelocity": null,
        "energyDensity": 4.0,
        "density": 1.60,
        "stability": 4,
        "criticalDiameter": null,
        "activationEnergy": 150,
        "gurvichTemperature": 3000,
        "chapmanJouguetPressure": null,
        "oxygenBalance": -30,
        "thermalConductivity": 0.20,
        "specificHeat": 1200,
        "thermalExpansion": 6.0e-5,
        "shockSensitivity": 100,
        "frictionSensitivity": 100,
        "impactSensitivity": 5,
        "expansionRatio": 600,
        "detonationProducts": [
            { "formula": "CO2", "moleFraction": 0.35, "heatCapacity": 37.1 },
            { "formula": "CO", "moleFraction": 0.25, "heatCapacity": 29.1 },
            { "formula": "H2O", "moleFraction": 0.25, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.15, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "Thermobaric Explosive (TBX)",
        "detonationVelocity": 5000,
        "energyDensity": 6.0,
        "density": 1.50,
        "stability": 3,
        "criticalDiameter": 5.0,
        "activationEnergy": 180,
        "gurvichTemperature": 4000,
        "chapmanJouguetPressure": 20,
        "oxygenBalance": -50,
        "thermalConductivity": 0.25,
        "specificHeat": 1300,
        "thermalExpansion": 7.0e-5,
        "shockSensitivity": 90,
        "frictionSensitivity": 150,
        "impactSensitivity": 10,
        "expansionRatio": 800,
        "detonationProducts": [
            { "formula": "CO2", "moleFraction": 0.40, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.35, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.15, "heatCapacity": 29.1 },
            { "formula": "CO", "moleFraction": 0.10, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "TATB",
        "detonationVelocity": 7300,
        "energyDensity": 4.1,
        "density": 1.93,
        "stability": 5,
        "criticalDiameter": 10.0,
        "activationEnergy": 220,
        "gurvichTemperature": 3000,
        "chapmanJouguetPressure": 28,
        "oxygenBalance": -54,
        "thermalConductivity": 0.32,
        "specificHeat": 1180,
        "thermalExpansion": 7.0e-5,
        "shockSensitivity": 100,
        "frictionSensitivity": 300,
        "impactSensitivity": 30,
        "expansionRatio": 650,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.25, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.22, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.18, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.35, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "TATP",
        "detonationVelocity": 5300,
        "energyDensity": 4.0,
        "density": 1.20,
        "stability": 1,
        "criticalDiameter": 1.0,
        "activationEnergy": 110,
        "gurvichTemperature": 2800,
        "chapmanJouguetPressure": 6,
        "oxygenBalance": 0,
        "thermalConductivity": 0.2,
        "specificHeat": 1200,
        "thermalExpansion": 1.0e-4,
        "shockSensitivity": 200,
        "frictionSensitivity": 150,
        "impactSensitivity": 2,
        "expansionRatio": 620,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.25, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.45, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.25, "heatCapacity": 33.6 },
            { "formula": "O2", "moleFraction": 0.05, "heatCapacity": 29.4 }
        ]
    },
    {
        "name": "HMTD",
        "detonationVelocity": 4100,
        "energyDensity": 3.5,
        "density": 1.57,
        "stability": 1,
        "criticalDiameter": 1.0,
        "activationEnergy": 120,
        "gurvichTemperature": 2600,
        "chapmanJouguetPressure": 5,
        "oxygenBalance": 0,
        "thermalConductivity": 0.18,
        "specificHeat": 1160,
        "thermalExpansion": 9.0e-5,
        "shockSensitivity": 180,
        "frictionSensitivity": 120,
        "impactSensitivity": 3,
        "expansionRatio": 550,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.20, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.40, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.30, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.10, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "CL-20",
        "detonationVelocity": 9400,
        "energyDensity": 6.4,
        "density": 2.04,
        "stability": 3,
        "criticalDiameter": 0.5,
        "activationEnergy": 220,
        "gurvichTemperature": 4600,
        "chapmanJouguetPressure": 47,
        "oxygenBalance": -10,
        "thermalConductivity": 0.31,
        "specificHeat": 1250,
        "thermalExpansion": 7.5e-5,
        "shockSensitivity": 85,
        "frictionSensitivity": 160,
        "impactSensitivity": 7,
        "expansionRatio": 800,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.18, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.40, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.25, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.17, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "FOX-7",
        "detonationVelocity": 8950,
        "energyDensity": 5.0,
        "density": 1.88,
        "stability": 4,
        "criticalDiameter": 1.0,
        "activationEnergy": 180,
        "gurvichTemperature": 4300,
        "chapmanJouguetPressure": 34,
        "oxygenBalance": -22,
        "thermalConductivity": 0.3,
        "specificHeat": 1220,
        "thermalExpansion": 7.6e-5,
        "shockSensitivity": 70,
        "frictionSensitivity": 80,
        "impactSensitivity": 3,
        "expansionRatio": 750,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.25, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.35, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.25, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.15, "heatCapacity": 29.1 }
        ]
    },
    {
        "name": "FOX-12",
        "detonationVelocity": 8100,
        "energyDensity": 5.2,
        "density": 1.74,
        "stability": 4,
        "criticalDiameter": 1.0,
        "activationEnergy": 200,
        "gurvichTemperature": 4200,
        "chapmanJouguetPressure": 30,
        "oxygenBalance": -25,
        "thermalConductivity": 0.28,
        "specificHeat": 1230,
        "thermalExpansion": 7.8e-5,
        "shockSensitivity": 80,
        "frictionSensitivity": 100,
        "impactSensitivity": 5,
        "expansionRatio": 760,
        "detonationProducts": [
            { "formula": "CO", "moleFraction": 0.22, "heatCapacity": 29.1 },
            { "formula": "CO2", "moleFraction": 0.36, "heatCapacity": 37.1 },
            { "formula": "H2O", "moleFraction": 0.27, "heatCapacity": 33.6 },
            { "formula": "N2", "moleFraction": 0.15, "heatCapacity": 29.1 }
        ]
    }
];

// Enhanced casing materials database
const casingMaterials = [
    {
        material: 'Steel (mild)',
        density: 7.85,
        yieldStrength: 250,
        ultimateStrength: 420,
        elasticModulus: 200,
        hardness: 130,
        elongation: 25,
        thermalConductivity: 45,
        meltingPoint: 1723,
        fractureToughness: 50,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 4,
        fragmentationPattern: 'natural'
    },
    {
        material: 'Steel (high strength)',
        density: 7.85,
        yieldStrength: 550,
        ultimateStrength: 800,
        elasticModulus: 210,
        hardness: 280,
        elongation: 15,
        thermalConductivity: 40,
        meltingPoint: 1723,
        fractureToughness: 80,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 5,
        fragmentationPattern: 'natural'
    },
    {
        material: 'Aluminum',
        density: 2.7,
        yieldStrength: 240,
        ultimateStrength: 310,
        elasticModulus: 70,
        hardness: 75,
        elongation: 12,
        thermalConductivity: 167,
        meltingPoint: 933,
        fractureToughness: 24,
        poissonRatio: 0.33,
        heatCapacity: 900,
        corrosionResistance: 7,
        fragmentationPattern: 'natural'
    },
    {
        material: 'Titanium Alloy',
        density: 4.5,
        yieldStrength: 880,
        ultimateStrength: 950,
        elasticModulus: 110,
        hardness: 320,
        elongation: 14,
        thermalConductivity: 7.2,
        meltingPoint: 1941,
        fractureToughness: 100,
        poissonRatio: 0.34,
        heatCapacity: 530,
        corrosionResistance: 9,
        fragmentationPattern: 'controlled'
    }
];

// Example environment database
const environments = [
    {
        name: 'Standard Atmosphere',
        temperature: 293.15, // K (20°C)
        pressure: 101.325, // kPa
        humidity: 60, // %
        altitude: 0, // m
        reflectiveSurfaces: false
    },
    {
        name: 'Desert',
        temperature: 318.15, // K (45°C)
        pressure: 101.325, // kPa
        humidity: 15, // %
        altitude: 300, // m
        reflectiveSurfaces: false,
        soilType: 'sandy',
        soilDensity: 1.6 // g/cm³
    },
    {
        name: 'Arctic',
        temperature: 253.15, // K (-20°C)
        pressure: 101.325, // kPa
        humidity: 80, // %
        altitude: 0, // m
        reflectiveSurfaces: true,
        soilType: 'frozen',
        soilDensity: 1.9 // g/cm³
    },
    {
        name: 'Underwater',
        temperature: 283.15, // K (10°C)
        pressure: 202.65, // kPa (at 10m depth)
        humidity: 100, // %
        altitude: -10, // m
        reflectiveSurfaces: false,
        waterDepth: 10 // m
    }
];

// Example of comprehensive output for the existing scenarios
function runComprehensiveScenario(
    name: string,
    radius: number,
    shape: Shape,
    maxWeight: number,
    isPortable: boolean,
    casing: Casing,
    environment: Environment,
    targets?: Target[]
) {
    // Create box-drawing characters and formatting helpers
    const h = "═", v = "║", tl = "╔", tr = "╗", bl = "╚", br = "╝";
    const ltee = "╠", rtee = "╣", ttee = "╦", btee = "╩", cross = "╬";
    const dash = "─", pipe = "│", plus = "+", space = " ";

    // Helper function to create a header box
    const createHeader = (text: string, width: number = 80) => {
        // Calculate minimum required width based on content
        const minWidth = text.length + 4; // 2 chars padding on each side
        const actualWidth = Math.max(width, minWidth);

        const padding = actualWidth - text.length - 4;
        const leftPad = Math.floor(padding / 2);
        const rightPad = padding - leftPad;
        const headerLine = `${tl}${h.repeat(actualWidth - 2)}${tr}`;
        const contentLine = `${v}${space.repeat(leftPad + 1)}${text}${space.repeat(rightPad + 1)}${v}`;
        const footerLine = `${bl}${h.repeat(actualWidth - 2)}${br}`;
        console.log(headerLine);
        console.log(contentLine);
        console.log(footerLine);
        return actualWidth; // Return the actual width used for consistency
    };

    // Helper function to create a section header
    const createSection = (text: string, width: number = 80) => {
        // Calculate minimum required width based on content
        const minWidth = text.length + 6; // For "═ " and " ═"
        const actualWidth = Math.max(width, minWidth);

        const padding = actualWidth - text.length - 4;
        const sectionLine = `${ltee}${h.repeat(2)}${space}${text}${space}${h.repeat(Math.max(2, padding - 2))}${rtee}`;
        console.log(sectionLine);
        return actualWidth; // Return the actual width used
    };

    // Ensure consistent box bottoms/separators
    const createSeparator = (width: number = 80) => {
        const line = `${v}${space.repeat(width - 2)}${v}`;
        console.log(line);
        return width;
    };

    // Also update the bottom border generator
    const createBottomBorder = (width: number = 80) => {
        const line = `${bl}${h.repeat(width - 2)}${br}`;
        console.log(line);
        return width;
    };

    // Helper function to display a value with proper alignment
    const displayValue = (label: string, value: string | number, unit: string = "", width: number = 80) => {
        const formattedValue = typeof value === 'number' ?
            (Number.isInteger(value) ? value.toString() : value.toFixed(2)) :
            value;
        const valueText = `${formattedValue}${unit ? ' ' + unit : ''}`;

        // Calculate minimum required width based on content
        const contentWidth = label.length + valueText.length + 8; // 8 for spacing/formatting
        const actualWidth = Math.max(width, contentWidth);

        // Ensure proper dot spacing that fills the available width
        const dotWidth = Math.max(2, actualWidth - label.length - valueText.length - 8);
        const dots = dash.repeat(dotWidth);

        const valueLine = `${v}  ${label} ${dots} ${valueText}  ${v}`;
        console.log(valueLine);
        return actualWidth; // Return the actual width used
    };

    // Helper function to create a data table
    const createTable = (headers: string[], rows: string[][], widths: number[], desiredWidth: number = 80) => {
        // First determine the actual column widths needed based on content
        const calculatedWidths = widths.map((minWidth, idx) => {
            // Check header width
            let colMaxWidth = headers[idx].length + 2; // +2 for padding

            // Check each row's content width for this column
            rows.forEach(row => {
                if (row[idx]) {
                    colMaxWidth = Math.max(colMaxWidth, row[idx].length + 2);
                }
            });

            // Ensure minimum width is respected
            return Math.max(minWidth, colMaxWidth);
        });

        // Calculate total table width without adjustments
        let tableWidth = calculatedWidths.reduce((sum, w) => sum + w, 0) + calculatedWidths.length + 1;

        // Adjust column widths to match desired width
        if (tableWidth != desiredWidth) {
            const widthDifference = desiredWidth - tableWidth;
            const totalAdjustableWidth = calculatedWidths.reduce((sum, w) => sum + w, 0);

            // Distribute the width difference proportionally across columns
            calculatedWidths.forEach((width, idx) => {
                const adjustmentRatio = width / totalAdjustableWidth;
                const adjustment = Math.floor(widthDifference * adjustmentRatio);
                calculatedWidths[idx] = Math.max(3, width + adjustment); // Ensure minimum of 3 chars per column
            });

            // Recalculate table width after adjustments
            tableWidth = calculatedWidths.reduce((sum, w) => sum + w, 0) + calculatedWidths.length + 1;

            // Handle any remaining difference due to rounding (add to the last column)
            const remainingDifference = desiredWidth - tableWidth;
            if (remainingDifference !== 0 && calculatedWidths.length > 0) {
                calculatedWidths[calculatedWidths.length - 1] += remainingDifference;
            }

            // Final table width should now match desired width
            tableWidth = desiredWidth;
        }

        // Table top border
        let line = `${ltee}`;
        for (let i = 0; i < calculatedWidths.length; i++) {
            line += h.repeat(calculatedWidths[i]);
            line += i < calculatedWidths.length - 1 ? ttee : rtee;
        }
        console.log(line);

        // Headers
        line = `${v}`;
        for (let i = 0; i < headers.length; i++) {
            const padding = calculatedWidths[i] - headers[i].length;
            line += ` ${headers[i]}${space.repeat(padding - 1)}${v}`;
        }
        console.log(line);

        // Separator
        line = `${ltee}`;
        for (let i = 0; i < calculatedWidths.length; i++) {
            line += h.repeat(calculatedWidths[i]);
            line += i < calculatedWidths.length - 1 ? cross : rtee;
        }
        console.log(line);

        // Data rows
        for (const row of rows) {
            line = `${v}`;
            for (let i = 0; i < row.length; i++) {
                const content = row[i] || "";
                const padding = calculatedWidths[i] - content.length;
                line += ` ${content}${space.repeat(padding - 1)}${v}`;
            }
            console.log(line);
        }

        // Table bottom border
        line = `${bl}`;
        for (let i = 0; i < calculatedWidths.length; i++) {
            line += h.repeat(calculatedWidths[i]);
            line += i < calculatedWidths.length - 1 ? btee : br;
        }
        console.log(line);
        return tableWidth; // Return the actual table width
    };

    // Helper function to safely get a property with a default value
    const getMetricSafely = (obj: any, path: string, defaultValue: any = "N/A") => {
        const parts = path.split('.');
        let current = obj;

        for (const part of parts) {
            if (current === undefined || current === null) {
                return defaultValue;
            }
            current = current[part];
        }

        return current !== undefined && current !== null ? current : defaultValue;
    };

    // Set an initial width that gets adjusted based on content
    let displayWidth = 80;

    // Create a scenario title
    console.log("");
    displayWidth = createHeader(`SCENARIO: ${name.toUpperCase()}`, displayWidth);
    console.log("");

    const scenario = {
        radius,
        shape,
        maxWeight,
        isPortable,
        casing,
        environment,
        targets
    };

    const result = optimizeExplosive(explosives, scenario);

    if (result) {
        // Display scenario parameters
        displayWidth = createSection("SCENARIO PARAMETERS", displayWidth);
        displayWidth = displayValue("Shape Type", scenario.shape.type, "", displayWidth);
        displayWidth = displayValue("Maximum Weight", maxWeight, "kg", displayWidth);
        displayWidth = displayValue("Portability", isPortable ? "Yes" : "No", "", displayWidth);
        displayWidth = displayValue("Temperature", environment.temperature, "K", displayWidth);
        displayWidth = createSeparator(displayWidth);

        // Display best material
        displayWidth = createSection("OPTIMAL EXPLOSIVE", displayWidth);
        displayValue("Material", result.material.name);
        displayValue("Stability Rating", `${result.metrics.stability}/5`);
        displayValue("Detonation Velocity", result.metrics.detonationVelocity, "m/s");
        displayValue("Detonation Pressure", result.metrics.detonationPressure, "GPa");
        displayValue("Detonation Temperature", result.metrics.detonationTemperature, "K");
        displayWidth = createSeparator(displayWidth);

        // Display mass and energy
        displayWidth = createSection("MASS AND ENERGY PROFILE", displayWidth);
        displayValue("TNT Equivalence Factor", result.metrics.tntEquivalent);
        displayValue("Explosive Mass", result.metrics.explosiveMass, "kg");
        displayValue("Casing Mass", result.metrics.casingMass, "kg");
        displayValue("Total Weight", result.metrics.totalWeight, "kg");
        displayValue("Total Energy", result.metrics.explosionEnergy, "MJ");
        displayWidth = createSeparator(displayWidth);

        // Display blast effects
        displayWidth = createSection(`BLAST EFFECTS AT ${radius}m`, displayWidth);
        displayValue("Peak Overpressure", result.metrics.overpressure, "kPa");
        displayValue("Positive Phase Duration", result.metrics.positivePhaseDuration, "ms");
        displayValue("Positive Phase Impulse", result.metrics.positivePhaseImpulse, "kPa·ms");
        displayValue("Reflected Pressure", result.metrics.reflectedPressure, "kPa");
        displayWidth = createSeparator(displayWidth);

        // Display fragmentation effects
        displayWidth = createSection("FRAGMENTATION PROFILE", displayWidth);
        displayValue("Initial Fragment Velocity", result.metrics.fragmentVelocity.initialVelocity, "m/s");
        displayValue("Terminal Fragment Velocity", result.metrics.fragmentVelocity.terminalVelocity, "m/s");
        displayValue("Average Fragment Mass", result.metrics.averageFragmentMass, "g");
        displayValue("Fragment Count", Math.round(result.metrics.numberOfFragments));
        displayValue("Maximum Fragment Size", result.metrics.maxFragmentSize, "cm");
        displayWidth = createSeparator(displayWidth);

        // Display thermal effects
        displayWidth = createSection("THERMAL EFFECTS", displayWidth);
        displayValue(`Thermal Radiation at ${radius}m`, getMetricSafely(result.metrics, "thermalRadiation"), "J/cm²");
        displayValue("Fireball Radius", getMetricSafely(result.metrics, "fireball.radius"), "m");
        displayValue("Fireball Duration", getMetricSafely(result.metrics, "fireball.duration"), "s");
        displayWidth = createSeparator(displayWidth);

        // Display ground shock effects if available
        if (getMetricSafely(result.metrics, "groundShock") !== "N/A") {
            displayWidth = createSection(`GROUND SHOCK AT ${radius}m`, displayWidth);
            displayValue("Peak Particle Velocity", getMetricSafely(result.metrics, "groundShock.peakParticleVelocity"), "m/s");
            displayValue("Peak Acceleration", getMetricSafely(result.metrics, "groundShock.peakAcceleration"), "g");
            displayValue("Peak Displacement", getMetricSafely(result.metrics, "groundShock.peakDisplacement"), "mm");
            displayWidth = createSeparator(displayWidth);
        }

        // Display underwater effects if available
        if (getMetricSafely(result.metrics, "underwaterEffects") !== "N/A") {
            displayWidth = createSection("UNDERWATER EFFECTS", displayWidth);
            displayValue("Peak Pressure", getMetricSafely(result.metrics, "underwaterEffects.peakPressure"), "MPa");
            displayValue("Maximum Bubble Radius", getMetricSafely(result.metrics, "underwaterEffects.bubblePulsation.maxRadius"), "m");
            displayValue("First Bubble Period", getMetricSafely(result.metrics, "underwaterEffects.bubblePulsation.firstPeriod"), "s");
            displayValue("Second Bubble Period", getMetricSafely(result.metrics, "underwaterEffects.bubblePulsation.secondPeriod"), "s");
            displayWidth = createSeparator(displayWidth);
        }

        // Display shaped charge effects if available
        if (getMetricSafely(result.metrics, "shapedChargeEffects") !== "N/A") {
            displayWidth = createSection("SHAPED CHARGE PERFORMANCE", displayWidth);
            displayValue("Jet Velocity", getMetricSafely(result.metrics, "shapedChargeEffects.jetVelocity"), "m/s");
            displayValue("Penetration Depth", getMetricSafely(result.metrics, "shapedChargeEffects.penetrationDepth"), "cm");
            displayValue("Jet Length", getMetricSafely(result.metrics, "shapedChargeEffects.jetLength"), "cm");
            displayWidth = createSeparator(displayWidth);
        }

        // Display target effects if available
        if (getMetricSafely(result.metrics, "targetsAssessment") !== "N/A" && result.metrics.targetsAssessment.length > 0) {
            displayWidth = createSection("TARGET EFFECTS ASSESSMENT", displayWidth);

            result.metrics.targetsAssessment.forEach((assessment: any, index: number) => {
                const targetName = `Target ${index + 1}: ${assessment.target.material}`;
                displayValue(targetName, `at ${assessment.target.distance}m`);
                displayValue("  Peak Overpressure", getMetricSafely(assessment, "blastParams.overpressure"), "kPa");
                displayValue("  Positive Phase Impulse", getMetricSafely(assessment, "blastParams.positivePhaseImpulse"), "kPa·ms");

                if (getMetricSafely(assessment, "injury") !== "N/A") {
                    displayValue("  Lethality Probability", getMetricSafely(assessment, "injury.lethality"), "%");
                    displayValue("  Lung Damage", getMetricSafely(assessment, "injury.lungDamage"));
                    displayValue("  Eardrum Rupture Probability", getMetricSafely(assessment, "injury.eardrumRupture"), "%");
                }

                if (getMetricSafely(assessment, "damage") !== "N/A") {
                    displayValue("  Structural Damage", getMetricSafely(assessment, "damage.damagePercentage"), "%");
                    displayValue("  Damage Mode", getMetricSafely(assessment, "damage.damageMode"));
                }

                if (index < result.metrics.targetsAssessment.length - 1) {
                    displayWidth = createSeparator(displayWidth);
                }
            });

            displayWidth = createSeparator(displayWidth);
        }

        // Display safety summary
        displayWidth = createSection("SAFETY AND HAZARD SUMMARY", displayWidth);

        // Create safety rating (calculated from stability and sensitivity values)
        const safetyScore = 5 - result.metrics.stability; // Invert so lower is better
        let safetyRating = "";
        if (safetyScore <= 1) safetyRating = "VERY SAFE";
        else if (safetyScore <= 2) safetyRating = "SAFE";
        else if (safetyScore <= 3) safetyRating = "MODERATE HAZARD";
        else if (safetyScore <= 4) safetyRating = "DANGEROUS";
        else safetyRating = "EXTREME HAZARD";

        displayValue("Safety Rating", safetyRating);
        displayValue("Handling Precautions", result.metrics.stability <= 2 ? "Standard" : "Special");
        displayWidth = createSeparator(displayWidth);

        // Close the bottom of the display
        displayWidth = createBottomBorder(displayWidth);
        console.log(""); // Add final line break
    } else {
        let maxWidth = 80;

        // Run parameter optimization to find a feasible scenario
        const optimization = optimizeScenarioParameters(
            explosives,
            { radius, shape, maxWeight, isPortable, casing, environment, targets }
        );

        const headerText = `OPTIMIZATION REQUIRED: No suitable material found for "${name}"`;

        // Calculate minimum width required for the header text
        const headerMinWidth = headerText.length; // 2 chars padding on each side
        maxWidth = Math.max(maxWidth, headerMinWidth);

        if (optimization?.recommendedChanges.length > 0) {
            // Pre-calculate width needed for all parameter changes to ensure consistency
            for (const change of optimization.recommendedChanges) {
                const label = change.parameter;
                const valueText = `${change.originalValue} → ${change.newValue} ${change.unit}`;
                maxWidth = Math.max(maxWidth, label.length + valueText.length + 8);
            }
        }

        // Add some buffer room for good measure
        maxWidth += 4;

        // Create header for optimization required message
        createHeader(headerText, maxWidth);

        if (optimization.feasible) {
            // Display recommended changes section
            createSection("Recommended Parameter Adjustments", maxWidth);

            // Display recommended changes
            for (const change of optimization.recommendedChanges) {
                if (change.originalValue !== change.newValue) {
                    displayValue(change.parameter, `${change.originalValue} → ${typeof change.newValue === 'number' ? Number.isInteger(change.newValue) ? change.newValue.toString() : change.newValue.toFixed(2) : change.newValue}`, change.unit, maxWidth);
                }
            }

            // Display optimization metrics section
            createSection("Optimization Metrics", maxWidth);
            displayValue("Confidence", `${optimization.confidence.toFixed(2)}%`, "", maxWidth);
            displayValue("Required Changes", optimization.recommendedChanges.length, "", maxWidth);
            displayValue("Estimated Material", optimization.estimatedMaterial.name, "", maxWidth);

            // Display sensitivity analysis section
            createSection("Parameter Sensitivity Analysis", maxWidth);
            const headers = ["Parameter", "Sensitivity", "Feasibility Impact"];
            const rows = optimization.sensitivityAnalysis.map(item =>
                [item.parameter, item.sensitivity.toFixed(3), item.impact]
            );
            createTable(headers, rows, [40, 15, 25], maxWidth);
        } else {
            // Display optimization failed section
            createSection("Optimization Failed", maxWidth);
            console.log(`${v} No feasible parameter combination could be found with the available explosive ${v}`);
            console.log(`${v} materials. Consider acquiring specialized materials or fundamentally      ${v}`);
            console.log(`${v} redesigning the application.                                              ${v}`);
            createBottomBorder(maxWidth);
        }
    }
}

/**
 * Optimizes scenario parameters to make it feasible with available explosive materials
 * using advanced multi-objective optimization techniques and physical constraints.
 */
function optimizeScenarioParameters(
    availableExplosives: Explosive[],
    scenario: {
        radius: number,
        shape: Shape,
        maxWeight: number,
        isPortable: boolean,
        casing: Casing,
        environment: Environment,
        targets?: Target[]
    }
): {
    feasible: boolean;
    recommendedChanges: Array<{
        parameter: string;
        originalValue: number | string;
        newValue: number | string;
        unit: string;
    }>;
    confidence: number;
    estimatedMaterial: Explosive;
    sensitivityAnalysis: Array<{
        parameter: string;
        sensitivity: number;
        impact: string;
    }>;
} {
    // Deep clone the scenario to avoid modifying the original
    const originalScenario = JSON.parse(JSON.stringify(scenario));
    let optimizedScenario = JSON.parse(JSON.stringify(scenario));
    const recommendedChanges: Array<{ parameter: string; originalValue: number | string; newValue: number | string; unit: string }> = [];

    // Initialize sensitivity analysis tracking
    const sensitivityTracking: Record<string, { sensitivity: number, impact: string }> = {};

    // ==== Step 1: Analyze constraint violations ====
    const constraintAnalysis = analyzeConstraintViolations(availableExplosives, originalScenario);

    if (constraintAnalysis.criticalConstraints.length === 0) {
        // This shouldn't happen - if there are no violations, optimization should have succeeded
        return {
            feasible: false,
            recommendedChanges: [],
            confidence: 0,
            estimatedMaterial: availableExplosives[0],
            sensitivityAnalysis: []
        };
    }

    // ==== Step 2: Apply Multi-objective optimization ====
    // Using a combination of gradient descent, genetic algorithms, and physics-based heuristics

    // Start with most critical parameters identified in constraint analysis
    for (const constraint of constraintAnalysis.criticalConstraints) {
        const parameterOpts = generateParameterOptions(constraint, originalScenario, availableExplosives);

        // Try each parameter option and evaluate its effectiveness and physical feasibility
        for (const option of parameterOpts) {
            // Apply the parameter change
            applyParameterChange(optimizedScenario, option.parameter, option.suggestedValue);

            // Track sensitivity
            const sensitivity = evaluateParameterSensitivity(
                originalScenario,
                optimizedScenario,
                option.parameter,
                availableExplosives
            );
            sensitivityTracking[option.parameter] = {
                sensitivity,
                impact: describeSensitivityImpact(sensitivity)
            };

            // Verify physical consistency of the changes
            if (!isPhysicallyConsistent(optimizedScenario)) {
                // Revert change if physically inconsistent
                applyParameterChange(optimizedScenario, option.parameter, originalParameterValue(originalScenario, option.parameter));
                continue;
            }

            // Check if this change improves feasibility
            if (improvesFeasibility(originalScenario, optimizedScenario, option.parameter, availableExplosives)) {
                recommendedChanges.push({
                    parameter: option.parameterName,
                    originalValue: originalParameterValue(originalScenario, option.parameter),
                    newValue: option.suggestedValue,
                    unit: option.unit
                });

                // Break if we've achieved feasibility
                const result = optimizeExplosive(availableExplosives, optimizedScenario);
                if (result) {
                    return {
                        feasible: true,
                        recommendedChanges,
                        confidence: calculateConfidence(recommendedChanges, originalScenario, optimizedScenario),
                        estimatedMaterial: result.material,
                        sensitivityAnalysis: Object.entries(sensitivityTracking).map(([param, data]) => ({
                            parameter: param,
                            sensitivity: data.sensitivity,
                            impact: data.impact
                        }))
                    };
                }
            } else {
                // Revert change if it doesn't improve feasibility
                applyParameterChange(optimizedScenario, option.parameter, originalParameterValue(originalScenario, option.parameter));
            }
        }
    }

    // ==== Step 3: Apply advanced optimization techniques if simple changes didn't work ====
    // If still not feasible, employ particle swarm optimization or simulated annealing
    const { improved, newScenario, changes } = applyAdvancedOptimization(
        originalScenario,
        recommendedChanges,
        availableExplosives,
        sensitivityTracking
    );

    if (improved) {
        optimizedScenario = newScenario;
        recommendedChanges.push(...changes);

        const result = optimizeExplosive(availableExplosives, optimizedScenario);
        if (result) {
            return {
                feasible: true,
                recommendedChanges,
                confidence: calculateConfidence(recommendedChanges, originalScenario, optimizedScenario),
                estimatedMaterial: result.material,
                sensitivityAnalysis: Object.entries(sensitivityTracking).map(([param, data]) => ({
                    parameter: param,
                    sensitivity: data.sensitivity,
                    impact: data.impact
                }))
            };
        }
    }

    // ==== Step 4: If all optimizations fail, find the closest possible scenario ====
    const closestScenario = findClosestFeasibleScenario(originalScenario, availableExplosives);
    if (closestScenario.distance < Number.MAX_VALUE) {
        return {
            feasible: true,
            recommendedChanges: closestScenario.changes,
            confidence: 70 - (closestScenario.distance * 10), // Lower confidence based on distance
            estimatedMaterial: closestScenario.material,
            sensitivityAnalysis: Object.entries(sensitivityTracking).map(([param, data]) => ({
                parameter: param,
                sensitivity: data.sensitivity,
                impact: data.impact
            }))
        };
    }

    return {
        feasible: false,
        recommendedChanges,
        confidence: 0,
        estimatedMaterial: availableExplosives[0], // Default to first available as fallback
        sensitivityAnalysis: Object.entries(sensitivityTracking).map(([param, data]) => ({
            parameter: param,
            sensitivity: data.sensitivity,
            impact: data.impact
        }))
    };
}

/**
 * Analyzes which constraints are preventing optimization success
 */
function analyzeConstraintViolations(
    explosives: Explosive[],
    scenario: any
): {
    criticalConstraints: Array<{ type: string, severity: number, parameter: string, parameterName?: string }>
} {
    const criticalConstraints: Array<{ type: string, severity: number, parameter: string, parameterName?: string }> = [];

    // Analyze explosion yield requirements vs available materials
    const yieldRequirements = estimateRequiredYield(scenario);
    const maxAvailableYield = Math.max(...explosives.map(e =>
        e.energyDensity * scenario.maxWeight
    ));

    if (yieldRequirements.totalRequiredEnergy > maxAvailableYield) {
        criticalConstraints.push({
            type: "energy_yield",
            severity: yieldRequirements.totalRequiredEnergy / maxAvailableYield,
            parameter: "maxWeight",
            parameterName: "Maximum Weight"
        });
    }

    // Check density constraints
    const volumeConstraint = calculateAvailableVolume(scenario.radius, scenario.shape);
    const minRequiredDensity = scenario.maxWeight / volumeConstraint;
    const maxAvailableDensity = Math.max(...explosives.map(e => e.density));

    if (minRequiredDensity > maxAvailableDensity) {
        criticalConstraints.push({
            type: "density",
            severity: minRequiredDensity / maxAvailableDensity,
            parameter: "radius",
            parameterName: "Radius"
        });
    }

    // Check detonation velocity requirements
    if (scenario.targets) {
        const minRequiredVelocity = estimateRequiredDetonationVelocity(scenario);
        const maxAvailableVelocity = Math.max(...explosives.map(e => e.detonationVelocity));

        if (minRequiredVelocity > maxAvailableVelocity) {
            criticalConstraints.push({
                type: "detonation_velocity",
                severity: minRequiredVelocity / maxAvailableVelocity,
                parameter: "shape.type",
                parameterName: "Shape Type"
            });
        }
    }

    // Check shape charge parameters if applicable
    if (scenario.shape.type === 'shaped_charge' && scenario.shape.linerProperties) {
        const optimalStandoff = ExplosiveModel.calculateOptimalStandoff(scenario.shape.linerProperties);
        const currentStandoff = scenario.shape.dimensions.standoff || 0;

        if (Math.abs(currentStandoff - optimalStandoff) / optimalStandoff > 0.3) {
            criticalConstraints.push({
                type: "standoff_distance",
                severity: Math.abs(currentStandoff - optimalStandoff) / optimalStandoff,
                parameter: "shape.dimensions.standoff",
                parameterName: "Standoff Distance"
            });
        }
    }

    // Check casing property constraints
    const casingConstraints = analyzeOptimalCasingProperties(scenario, explosives);
    if (casingConstraints.nonOptimal) {
        casingConstraints.issues.forEach(issue => {
            criticalConstraints.push({
                type: `casing_${issue.property}`,
                severity: issue.severity,
                parameter: `casing.${issue.property}`,
                parameterName: `Casing ${issue.property.charAt(0).toUpperCase() + issue.property.slice(1)}`
            });
        });
    }

    // Check critical diameter constraints
    const minRequiredDiameter = estimateRequiredCriticalDiameter(scenario);
    const minAvailableDiameter = Math.min(...explosives.map(e => e.criticalDiameter));

    // Double the radius to get diameter and convert from cm to mm
    const scenarioDiameter = 2 * scenario.radius * 10;

    if (minRequiredDiameter > scenarioDiameter) {
        criticalConstraints.push({
            type: "critical_diameter",
            severity: minRequiredDiameter / scenarioDiameter,
            parameter: "radius",
            parameterName: "Radius (Critical Diameter)"
        });
    }

    // Check stability requirements for portable scenarios
    if (scenario.isPortable) {
        const minRequiredStability = 3; // Assume minimum stability for portable applications
        const materialWithSufficientStability = explosives.some(e =>
            e.stability >= minRequiredStability &&
            e.energyDensity * scenario.maxWeight >= yieldRequirements.totalRequiredEnergy
        );

        if (!materialWithSufficientStability) {
            criticalConstraints.push({
                type: "stability",
                severity: 1.5, // Fixed severity for this constraint
                parameter: "isPortable",
                parameterName: "Portability Requirement"
            });
        }
    }

    // Analyze environmental constraints
    if (scenario.environment.waterDepth !== undefined && scenario.environment.waterDepth > 0) {
        // Underwater explosion special constraints
        const pressureEffect = 1 + (scenario.environment.waterDepth * 9.81 * 1000 / 101325);
        if (pressureEffect > 2) {
            criticalConstraints.push({
                type: "water_pressure",
                severity: pressureEffect / 2,
                parameter: "environment.waterDepth",
                parameterName: "Water Depth"
            });
        }
    }

    // High altitude considerations
    if (scenario.environment.altitude > 2000) {
        const pressureDropFactor = Math.exp(-scenario.environment.altitude / 8000);
        const effectiveYieldReduction = 1 / pressureDropFactor;

        if (effectiveYieldReduction > 1.3) {
            criticalConstraints.push({
                type: "altitude_effect",
                severity: effectiveYieldReduction,
                parameter: "environment.altitude",
                parameterName: "Altitude"
            });
        }
    }

    // Temperature extreme considerations
    const standardTemp = 293; // K (20°C)
    if (Math.abs(scenario.environment.temperature - standardTemp) > 50) {
        const temperatureEffect = 1 + Math.abs(scenario.environment.temperature - standardTemp) / 100;

        criticalConstraints.push({
            type: "temperature_extreme",
            severity: temperatureEffect,
            parameter: "environment.temperature",
            parameterName: "Environmental Temperature"
        });
    }

    // Sort constraints by severity
    criticalConstraints.sort((a, b) => b.severity - a.severity);

    return { criticalConstraints };
}

/**
 * Estimates required detonation performance based on scenario
 */
function estimateRequiredYield(scenario: any): {
    totalRequiredEnergy: number,
    blastEnergy: number,
    fragmentationEnergy: number,
    penetrationEnergy: number
} {
    let totalRequiredEnergy = 0;
    let blastEnergy = 0;
    let fragmentationEnergy = 0;
    let penetrationEnergy = 0;

    // Calculate blast energy requirements
    if (scenario.targets) {
        for (const target of scenario.targets) {
            // Estimate overpressure at target distance
            const targetOverpressure = target.criticalDamageThreshold;
            const distanceCubed = Math.pow(target.distance, 3);

            // Simplified Hopkinson-Cranz scaling
            const requiredEnergy = (targetOverpressure * distanceCubed) /
                (ExplosiveModel.AIR_DENSITY * ExplosiveModel.ATMOSPHERIC_PRESSURE);

            blastEnergy = Math.max(blastEnergy, requiredEnergy);

            // Additional energy for penetration if applicable
            if (target.material !== 'human' && target.thickness > 0) {
                // Penetration energy model based on material properties
                const penetrationE = target.thickness * target.density *
                    target.crossSection * 5000; // Energy density for material deformation
                penetrationEnergy = Math.max(penetrationEnergy, penetrationE);
            }
        }
    }

    // Fragmentation energy depends on casing
    fragmentationEnergy = scenario.casing.density *
        scenario.casing.thickness *
        calculateSurfaceArea(scenario.radius, scenario.shape) *
        400; // Energy to accelerate fragments (J/g)

    totalRequiredEnergy = blastEnergy + fragmentationEnergy + penetrationEnergy;

    // Add safety margin
    totalRequiredEnergy *= 1.25;

    // Convert to MJ if in joules
    if (totalRequiredEnergy > 1000000) {
        totalRequiredEnergy /= 1000000;
    }

    return {
        totalRequiredEnergy,
        blastEnergy: blastEnergy > 1000000 ? blastEnergy / 1000000 : blastEnergy,
        fragmentationEnergy: fragmentationEnergy > 1000000 ? fragmentationEnergy / 1000000 : fragmentationEnergy,
        penetrationEnergy: penetrationEnergy > 1000000 ? penetrationEnergy / 1000000 : penetrationEnergy
    };
}

/**
 * Calculate surface area based on shape type
 */
function calculateSurfaceArea(radius: number, shape: Shape): number {
    switch (shape.type) {
        case 'sphere':
            return 4 * Math.PI * Math.pow(radius, 2);
        case 'cylinder': {
            const height = shape.dimensions.height || (2 * radius);
            return 2 * Math.PI * radius * (radius + height);
        }
        case 'cube': {
            const side = shape.dimensions.side || (2 * radius);
            return 6 * Math.pow(side, 2);
        }
        case 'hemisphere':
            return 3 * Math.PI * Math.pow(radius, 2);
        case 'conical': {
            const height = shape.dimensions.height || (2 * radius);
            const slantHeight = Math.sqrt(Math.pow(radius, 2) + Math.pow(height, 2));
            return Math.PI * radius * (radius + slantHeight);
        }
        case 'shaped_charge': {
            // Approximation for shaped charge with conical liner
            const height = shape.dimensions.height || (2 * radius);
            const coneAngle = shape.dimensions.coneAngle || 60; // Default angle if not specified
            const slantHeight = height / Math.cos(coneAngle * Math.PI / 180);
            return Math.PI * radius * (radius + slantHeight) * 1.1; // Add 10% for engineering elements
        }
        default:
            return 4 * Math.PI * Math.pow(radius, 2); // Default to sphere
    }
}

/**
 * Helper functions for the optimization process
 */
function calculateAvailableVolume(radius: number, shape: Shape): number {
    // Calculate volume based on shape type
    switch (shape.type) {
        case 'sphere':
            return (4 / 3) * Math.PI * Math.pow(radius, 3);
        case 'cylinder': {
            const height = shape.dimensions.height || (2 * radius);
            return Math.PI * Math.pow(radius, 2) * height;
        }
        case 'cube': {
            const side = shape.dimensions.side || (2 * radius);
            return Math.pow(side, 3);
        }
        case 'hemisphere':
            return (2 / 3) * Math.PI * Math.pow(radius, 3);
        case 'conical': {
            const height = shape.dimensions.height || (2 * radius);
            return (1 / 3) * Math.PI * Math.pow(radius, 2) * height;
        }
        case 'shaped_charge': {
            // Account for liner and other components in shaped charge
            const scRadius = radius;
            const scHeight = shape.dimensions.height || (2 * radius);
            let volume = 0.8 * Math.PI * Math.pow(scRadius, 2) * scHeight; // 80% of conical volume

            // Subtract liner volume if present
            if (shape.linerProperties) {
                const linerThickness = shape.linerProperties.thickness / 10; // mm to cm
                volume -= Math.PI * linerThickness * Math.pow(scRadius, 2) * 0.5; // Approximate liner volume
            }

            return volume;
        }
        default:
            return (4 / 3) * Math.PI * Math.pow(radius, 3); // Default to sphere
    }
}

/**
 * Estimate required critical diameter based on scenario parameters
 */
function estimateRequiredCriticalDiameter(scenario: any): number {
    // Base critical diameter on target effects
    let criticalDiameter = 5; // Default 5mm

    // Adjust based on target requirements if they exist
    if (scenario.targets) {
        for (const target of scenario.targets) {
            // Higher pressure requirements need more stable detonation
            if (target.criticalDamageThreshold > 200) { // kPa
                criticalDiameter = Math.max(criticalDiameter, 10 + 0.05 * target.criticalDamageThreshold);
            }

            // Distance affects requirements
            if (target.distance > 50) { // m
                criticalDiameter = Math.max(criticalDiameter, 15);
            }
        }
    }

    // Shape factors - shaped charges need more stable detonation
    if (scenario.shape.type === 'shaped_charge') {
        criticalDiameter = Math.max(criticalDiameter, 15);
    }

    // Environmental factors
    if (scenario.environment.temperature < 260) { // Cold conditions
        criticalDiameter *= 1.3; // Cold reduces sensitivity, needs larger diameter
    }

    // Cavity effects for hollow charges
    if (scenario.shape.dimensions.wallThickness) {
        criticalDiameter *= 1.5; // Hollow charges need larger diameters
    }

    return criticalDiameter;
}

/**
 * Estimate required detonation velocity based on scenario
 */
function estimateRequiredDetonationVelocity(scenario: any): number {
    let baseVelocity = 5000; // Base velocity in m/s

    // Adjust for target requirements
    if (scenario.targets) {
        for (const target of scenario.targets) {
            // Penetration of hard targets requires high detonation velocity
            if (target.material === 'concrete' || target.material === 'steel' || target.material === 'armor') {
                baseVelocity = Math.max(baseVelocity, 7000);

                // Thickness affects required velocity
                if (target.thickness > 5) { // cm
                    baseVelocity += target.thickness * 100; // Add 100 m/s per cm
                }
            }

            // Pressure requirements
            if (target.criticalDamageThreshold > 500) { // kPa
                baseVelocity = Math.max(baseVelocity, 6000 + target.criticalDamageThreshold / 100);
            }
        }
    }

    // Shape effects
    if (scenario.shape.type === 'shaped_charge') {
        baseVelocity = Math.max(baseVelocity, 7500);

        // Liner material affects required velocity
        if (scenario.shape.linerProperties) {
            if (scenario.shape.linerProperties.material === 'copper') {
                baseVelocity = Math.max(baseVelocity, 8000);
            } else if (scenario.shape.linerProperties.material === 'tantalum') {
                baseVelocity = Math.max(baseVelocity, 8500);
            }
        }
    }

    // Casing effects - heavy casings need higher velocity
    if (scenario.casing.density > 7.8) { // steel or heavier
        baseVelocity += 500 * (scenario.casing.density / 7.8 - 1);
    }

    // Environmental factors
    if (scenario.environment.waterDepth !== undefined && scenario.environment.waterDepth > 0) {
        // Underwater explosions benefit from high detonation velocity
        baseVelocity += 1000;
    }

    return baseVelocity;
}

/**
 * Analyze optimal casing properties based on the scenario and available explosives
 */
function analyzeOptimalCasingProperties(scenario: any, explosives: Explosive[]): {
    nonOptimal: boolean;
    issues: Array<{ property: string; severity: number; currentValue: number; recommendedValue: number }>;
} {
    const issues: Array<{ property: string; severity: number; currentValue: number; recommendedValue: number }> = [];

    // Function to check if a property is suboptimal
    const checkProperty = (property: string, currentValue: number, optimalValue: number, tolerance: number) => {
        const relativeDifference = Math.abs(currentValue - optimalValue) / optimalValue;
        if (relativeDifference > tolerance) {
            issues.push({
                property,
                severity: relativeDifference,
                currentValue,
                recommendedValue: optimalValue
            });
        }
    };

    // Get average properties of available explosives
    const avgDetonationVelocity = explosives.reduce((sum, exp) => sum + exp.detonationVelocity, 0) / explosives.length;
    const avgEnergyDensity = explosives.reduce((sum, exp) => sum + exp.energyDensity, 0) / explosives.length;
    const avgDensity = explosives.reduce((sum, exp) => sum + exp.density, 0) / explosives.length;

    // Calculate optimal casing thickness based on Gurney equations and explosive properties
    const optimalThickness = calculateOptimalCasingThickness(
        avgDetonationVelocity,
        avgEnergyDensity,
        scenario.radius,
        scenario.casing.density
    );

    checkProperty('thickness', scenario.casing.thickness, optimalThickness, 0.25);

    // Check density based on target effect requirements
    let optimalDensity = scenario.casing.density;
    if (scenario.targets) {
        // For penetration targets, higher density is better
        const penetrationTargets = scenario.targets.filter(t =>
            t.material !== 'human' && t.material !== 'wood' && t.thickness > 0
        );

        if (penetrationTargets.length > 0) {
            optimalDensity = Math.max(7.8, scenario.casing.density); // At least steel density
        }

        // For blast targets only, lighter casing may be more efficient
        const blastOnlyTargets = scenario.targets.filter(t =>
            (t.material === 'human' || t.material === 'structure') &&
            !penetrationTargets.length
        );

        if (blastOnlyTargets.length > 0 && penetrationTargets.length === 0) {
            optimalDensity = Math.min(7.8, scenario.casing.density); // Lighter than steel may be better
        }
    }

    checkProperty('density', scenario.casing.density, optimalDensity, 0.2);

    // Check optimal yield strength based on detonation pressure
    const avgChapmanJouguetPressure = explosives.reduce((sum, exp) => sum + exp.chapmanJouguetPressure, 0) / explosives.length;
    const optimalYieldStrength = avgChapmanJouguetPressure * 100; // Convert GPa to MPa with safety factor

    checkProperty('yieldStrength', scenario.casing.yieldStrength, optimalYieldStrength, 0.3);

    // Optimal fragmentation pattern based on targets
    if (scenario.targets) {
        let optimalPattern: 'uniform' | 'controlled' | 'natural' = 'natural';

        // For precision targets, controlled fragmentation is better
        const precisionTargets = scenario.targets.filter(t =>
            t.crossSection < 1 || (t.material !== 'human' && t.shielding && t.shielding > 0.5)
        );

        if (precisionTargets.length > 0) {
            optimalPattern = 'controlled';
        }

        // For area targets, uniform fragmentation is better
        const areaTargets = scenario.targets.filter(t =>
            t.crossSection > 10 || t.material === 'human'
        );

        if (areaTargets.length > 0 && precisionTargets.length === 0) {
            optimalPattern = 'uniform';
        }

        if (scenario.casing.fragmentationPattern !== optimalPattern) {
            issues.push({
                property: 'fragmentationPattern',
                severity: 1.0, // Fixed severity for categorical property
                currentValue: scenario.casing.fragmentationPattern as unknown as number,
                recommendedValue: optimalPattern as unknown as number
            });
        }
    }

    return {
        nonOptimal: issues.length > 0,
        issues
    };
}

/**
 * Calculate optimal casing thickness based on Gurney equations
 */
function calculateOptimalCasingThickness(
    detonationVelocity: number,
    energyDensity: number,
    radius: number,
    casingDensity: number
): number {
    // Estimate Gurney constant from detonation velocity and energy density
    // Gurney constant ≈ 0.3 * detonation velocity
    const gurneyConstant = 0.3 * detonationVelocity;

    // For maximum fragment velocity with cylindrical charge:
    // M/C = (V/√2E)² / (1 - (V/√2E)²)
    // Where:
    // M = mass of casing
    // C = mass of explosive
    // V = desired fragment velocity (about 0.6 * Gurney for optimal energy transfer)
    // E = Gurney energy constant

    const desiredVelocityRatio = 0.6;
    const velocityGurneyRatio = desiredVelocityRatio / Math.sqrt(2);
    const massRatio = Math.pow(velocityGurneyRatio, 2) / (1 - Math.pow(velocityGurneyRatio, 2));

    // M/C = (ρₘ*t*S) / (ρₑ*V)
    // Where:
    // ρₘ = density of metal casing
    // t = thickness of casing
    // S = surface area
    // ρₑ = density of explosive
    // V = volume of explosive

    // For cylindrical geometries, simplifies approximately to:
    // t ≈ massRatio * (ρₑ/ρₘ) * radius

    // Convert energy density (MJ/kg) to rough estimate of explosive density (g/cm³)
    const estimatedExplosiveDensity = 1.3 + 0.3 * energyDensity; // Empirical relationship

    // Calculate optimal thickness in cm
    return massRatio * (estimatedExplosiveDensity / casingDensity) * radius;
}

/**
 * Generate parameter adjustment options for a given constraint violation
 */
function generateParameterOptions(
    constraint: { type: string; severity: number; parameter: string; parameterName?: string },
    scenario: any,
    explosives: Explosive[]
): Array<{
    parameter: string;
    parameterName: string;
    suggestedValue: number | string;
    unit: string;
    impact: number;
}> {
    const options: Array<{
        parameter: string;
        parameterName: string;
        suggestedValue: any;
        unit: string;
        impact: number;
    }> = [];

    // Get current parameter value
    const currentValue = originalParameterValue(scenario, constraint.parameter);

    // Handle based on constraint type
    switch (constraint.type) {
        case 'energy_yield': {
            // Increase weight allowance
            const requiredYield = estimateRequiredYield(scenario);
            const bestEnergyDensity = Math.max(...explosives.map(e => e.energyDensity));
            const recommendedWeight = (requiredYield.totalRequiredEnergy / bestEnergyDensity) * 1.1; // 10% safety margin

            options.push({
                parameter: 'maxWeight',
                parameterName: 'Maximum Weight',
                suggestedValue: recommendedWeight,
                unit: 'kg',
                impact: 0.8
            });

            // Alternative: reduce target requirements
            if (scenario.targets) {
                // Find targets with highest pressure demands
                const highPressureTargets = scenario.targets
                    .filter(t => t.criticalDamageThreshold > 100)
                    .sort((a, b) => b.criticalDamageThreshold - a.criticalDamageThreshold);

                if (highPressureTargets.length > 0) {
                    const targetIndex = scenario.targets.indexOf(highPressureTargets[0]);
                    const reducedThreshold = highPressureTargets[0].criticalDamageThreshold * 0.8;

                    options.push({
                        parameter: `targets[${targetIndex}].criticalDamageThreshold`,
                        parameterName: `Target ${targetIndex + 1} Damage Threshold`,
                        suggestedValue: reducedThreshold,
                        unit: 'kPa',
                        impact: 0.6
                    });
                }
            }

            break;
        }

        case 'density': {
            // Increase radius
            const volumeNeeded = scenario.maxWeight / Math.max(...explosives.map(e => e.density));
            const recommendedRadius = Math.pow(volumeNeeded / ((4 / 3) * Math.PI), 1 / 3) * 1.05; // Sphere approximation + safety

            options.push({
                parameter: 'radius',
                parameterName: 'Radius',
                suggestedValue: recommendedRadius,
                unit: 'cm',
                impact: 0.9
            });

            // Alternative: change shape to more volume-efficient type
            if (scenario.shape.type !== 'cylinder') {
                options.push({
                    parameter: 'shape.type',
                    parameterName: 'Shape Type',
                    suggestedValue: 'cylinder',
                    unit: '',
                    impact: 0.7
                });

                // Add height parameter for cylinder
                options.push({
                    parameter: 'shape.dimensions.height',
                    parameterName: 'Cylinder Height',
                    suggestedValue: 2 * recommendedRadius,
                    unit: 'cm',
                    impact: 0.6
                });
            }

            break;
        }

        case 'detonation_velocity': {
            // For shaped charges, modify cone angle
            if (scenario.shape.type === 'shaped_charge') {
                const currentAngle = scenario.shape.dimensions.coneAngle || 60;
                let recommendedAngle: number;

                if (currentAngle < 45) {
                    recommendedAngle = 60; // Wider angle, less jet velocity but better overall detonation
                } else if (currentAngle > 70) {
                    recommendedAngle = 60; // Standard optimal angle
                } else {
                    recommendedAngle = 45; // Compromise for different applications
                }

                options.push({
                    parameter: 'shape.dimensions.coneAngle',
                    parameterName: 'Cone Angle',
                    suggestedValue: recommendedAngle,
                    unit: 'degrees',
                    impact: 0.75
                });
            } else {
                // Change to more efficient shape type
                options.push({
                    parameter: 'shape.type',
                    parameterName: 'Shape Type',
                    suggestedValue: 'shaped_charge',
                    unit: '',
                    impact: 0.9
                });

                // Add needed properties for shaped charge
                options.push({
                    parameter: 'shape.dimensions.coneAngle',
                    parameterName: 'Cone Angle',
                    suggestedValue: 60,
                    unit: 'degrees',
                    impact: 0.8
                });

                // Add liner properties if missing
                if (!scenario.shape.linerProperties) {
                    options.push({
                        parameter: 'shape.linerProperties',
                        parameterName: 'Liner Properties',
                        suggestedValue: {
                            material: 'copper',
                            density: 8.96,
                            thickness: 2.0,
                            angle: 42,
                            standoff: 2.0
                        },
                        unit: '',
                        impact: 0.85
                    });
                }
            }

            break;
        }

        case 'standoff_distance': {
            const optimalStandoff = ExplosiveModel.calculateOptimalStandoff(scenario.shape.linerProperties);

            options.push({
                parameter: 'shape.dimensions.standoff',
                parameterName: 'Standoff Distance',
                suggestedValue: optimalStandoff,
                unit: 'cm',
                impact: 0.95
            });

            break;
        }

        case 'critical_diameter': {
            // Increase radius to accommodate critical diameter
            const requiredDiameter = estimateRequiredCriticalDiameter(scenario);
            const recommendedRadius = (requiredDiameter / 20) + 0.5; // Convert mm to cm and add margin

            options.push({
                parameter: 'radius',
                parameterName: 'Radius',
                suggestedValue: Math.max(recommendedRadius, scenario.radius),
                unit: 'cm',
                impact: 0.85
            });

            break;
        }

        case 'stability': {
            // Make non-portable or increase safety factors
            options.push({
                parameter: 'isPortable',
                parameterName: 'Portability Requirement',
                suggestedValue: false,
                unit: '',
                impact: 0.9
            });

            break;
        }

        case 'water_pressure': {
            // Reduce water depth 
            const recommendedDepth = Math.max(1, scenario.environment.waterDepth * 0.7);

            options.push({
                parameter: 'environment.waterDepth',
                parameterName: 'Water Depth',
                suggestedValue: recommendedDepth,
                unit: 'm',
                impact: 0.75
            });

            break;
        }

        case 'altitude_effect': {
            // Reduce altitude or increase explosive mass
            const recommendedAltitude = Math.max(0, scenario.environment.altitude * 0.8);

            options.push({
                parameter: 'environment.altitude',
                parameterName: 'Altitude',
                suggestedValue: recommendedAltitude,
                unit: 'm',
                impact: 0.7
            });

            // Alternative: increase weight to compensate for altitude
            const pressureCorrection = Math.exp(-scenario.environment.altitude / 8000);
            const recommendedWeight = scenario.maxWeight / pressureCorrection;

            options.push({
                parameter: 'maxWeight',
                parameterName: 'Maximum Weight',
                suggestedValue: recommendedWeight,
                unit: 'kg',
                impact: 0.8
            });

            break;
        }

        case 'temperature_extreme': {
            // Adjust temperature closer to standard conditions
            const standardTemp = 293; // K (20°C)
            const currentTemp = scenario.environment.temperature;
            const recommendedTemp = currentTemp > standardTemp
                ? Math.max(standardTemp, currentTemp - 30)
                : Math.min(standardTemp, currentTemp + 30);

            options.push({
                parameter: 'environment.temperature',
                parameterName: 'Environmental Temperature',
                suggestedValue: recommendedTemp,
                unit: 'K',
                impact: 0.6
            });

            break;
        }
    }

    // Handle casing property constraints
    if (constraint.type.startsWith('casing_')) {
        const property = constraint.parameter.split('.')[1];

        // Get recommended value from casing analysis
        const casingAnalysis = analyzeOptimalCasingProperties(scenario, explosives);
        const propertyIssue = casingAnalysis.issues.find(issue => issue.property === property);

        if (propertyIssue) {
            options.push({
                parameter: constraint.parameter,
                parameterName: constraint.parameterName || property,
                suggestedValue: propertyIssue.recommendedValue,
                unit: getUnitForProperty(property),
                impact: 0.8
            });
        }
    }

    return options;
}

/**
 * Get unit for a specified property
 */
function getUnitForProperty(property: string): string {
    const unitMap: Record<string, string> = {
        'thickness': 'cm',
        'density': 'g/cm³',
        'yieldStrength': 'MPa',
        'ultimateStrength': 'MPa',
        'elasticModulus': 'GPa',
        'hardness': 'Brinell',
        'elongation': '%',
        'thermalConductivity': 'W/(m·K)',
        'meltingPoint': 'K',
        'fractureToughness': 'MPa·m^(1/2)',
        'poissonRatio': '',
        'heatCapacity': 'J/(kg·K)',
        'corrosionResistance': '1-10',
        'fragmentationPattern': '',
        'radius': 'cm',
        'height': 'cm',
        'standoff': 'cm',
        'coneAngle': 'degrees',
        'temperature': 'K',
        'pressure': 'kPa',
        'humidity': '%',
        'altitude': 'm',
        'waterDepth': 'm',
        'maxWeight': 'kg',
        'criticalDamageThreshold': 'kPa',
        'detonationVelocity': 'm/s',
        'energyDensity': 'MJ/kg',
        'stability': '1-5',
        'criticalDiameter': 'mm'
    };

    return unitMap[property] || '';
}

/**
 * Check if a parameter change improves feasibility
 */
function improvesFeasibility(
    originalScenario: any,
    modifiedScenario: any,
    modifiedParameter: string,
    explosives: Explosive[]
): boolean {
    // Clone the scenarios to avoid side effects
    const originalClone = JSON.parse(JSON.stringify(originalScenario));
    const modifiedClone = JSON.parse(JSON.stringify(modifiedScenario));

    // Measure feasibility through several metrics
    const originalMetrics = calculateFeasibilityMetrics(originalClone, explosives);
    const modifiedMetrics = calculateFeasibilityMetrics(modifiedClone, explosives);

    // Calculate overall improvement score
    let improvementScore = 0;

    // Compare energy requirements vs availability
    if (modifiedMetrics.energyDeficit < originalMetrics.energyDeficit) {
        improvementScore += 3 * (originalMetrics.energyDeficit - modifiedMetrics.energyDeficit) /
            originalMetrics.energyDeficit;
    }

    // Compare density requirements
    if (modifiedMetrics.densityDeficit < originalMetrics.densityDeficit) {
        improvementScore += 2 * (originalMetrics.densityDeficit - modifiedMetrics.densityDeficit) /
            Math.max(0.01, originalMetrics.densityDeficit);
    }

    // Compare detonation velocity requirements
    if (modifiedMetrics.velocityDeficit < originalMetrics.velocityDeficit) {
        improvementScore += 2 * (originalMetrics.velocityDeficit - modifiedMetrics.velocityDeficit) /
            Math.max(0.01, originalMetrics.velocityDeficit);
    }

    // Compare critical diameter requirements
    if (modifiedMetrics.diameterDeficit < originalMetrics.diameterDeficit) {
        improvementScore += (originalMetrics.diameterDeficit - modifiedMetrics.diameterDeficit) /
            Math.max(0.01, originalMetrics.diameterDeficit);
    }

    // For shape charge scenarios, consider standoff optimization
    if (originalScenario.shape.type === 'shaped_charge' &&
        originalScenario.shape.linerProperties &&
        modifiedParameter.includes('standoff')) {

        const optimalStandoff = ExplosiveModel.calculateOptimalStandoff(originalScenario.shape.linerProperties);
        const originalDiff = Math.abs(originalScenario.shape.dimensions.standoff - optimalStandoff) / optimalStandoff;
        const modifiedDiff = Math.abs(modifiedScenario.shape.dimensions.standoff - optimalStandoff) / optimalStandoff;

        if (modifiedDiff < originalDiff) {
            improvementScore += 2 * (originalDiff - modifiedDiff);
        }
    }

    // Check if there are now compatible explosives
    const originalCompatibleExplosives = explosives.filter(e => isExplosiveCompatible(e, originalClone));
    const modifiedCompatibleExplosives = explosives.filter(e => isExplosiveCompatible(e, modifiedClone));

    if (modifiedCompatibleExplosives.length > originalCompatibleExplosives.length) {
        improvementScore += 5;
    }

    // Threshold for considering an improvement
    return improvementScore > 0.2;
}

/**
 * Calculate metrics to evaluate feasibility
 */
function calculateFeasibilityMetrics(scenario: any, explosives: Explosive[]): {
    energyDeficit: number;
    densityDeficit: number;
    velocityDeficit: number;
    diameterDeficit: number;
    standoffError: number;
    compatibleExplosivesCount: number;
} {
    // Required energy
    const requiredYield = estimateRequiredYield(scenario);
    const maxAvailableYield = Math.max(...explosives.map(e => e.energyDensity * scenario.maxWeight));
    const energyDeficit = Math.max(0, requiredYield.totalRequiredEnergy - maxAvailableYield) /
        Math.max(1, requiredYield.totalRequiredEnergy);

    // Density requirements
    const volumeConstraint = calculateAvailableVolume(scenario.radius, scenario.shape);
    const minRequiredDensity = scenario.maxWeight / volumeConstraint;
    const maxAvailableDensity = Math.max(...explosives.map(e => e.density));
    const densityDeficit = Math.max(0, minRequiredDensity - maxAvailableDensity) /
        Math.max(0.1, minRequiredDensity);

    // Detonation velocity
    const minRequiredVelocity = estimateRequiredDetonationVelocity(scenario);
    const maxAvailableVelocity = Math.max(...explosives.map(e => e.detonationVelocity));
    const velocityDeficit = Math.max(0, minRequiredVelocity - maxAvailableVelocity) /
        Math.max(1000, minRequiredVelocity);

    // Critical diameter
    const minRequiredDiameter = estimateRequiredCriticalDiameter(scenario);
    const scenarioDiameter = 2 * scenario.radius * 10; // cm to mm
    const diameterDeficit = Math.max(0, minRequiredDiameter - scenarioDiameter) /
        Math.max(1, minRequiredDiameter);

    // Standoff error for shaped charges
    let standoffError = 0;
    if (scenario.shape.type === 'shaped_charge' && scenario.shape.linerProperties) {
        const optimalStandoff = ExplosiveModel.calculateOptimalStandoff(scenario.shape.linerProperties);
        const currentStandoff = scenario.shape.dimensions.standoff || 0;
        standoffError = Math.abs(currentStandoff - optimalStandoff) / Math.max(0.1, optimalStandoff);
    }

    // Count compatible explosives
    const compatibleExplosivesCount = explosives.filter(e => isExplosiveCompatible(e, scenario)).length;

    return {
        energyDeficit,
        densityDeficit,
        velocityDeficit,
        diameterDeficit,
        standoffError,
        compatibleExplosivesCount
    };
}

/**
 * Check if an explosive is compatible with the scenario requirements
 */
function isExplosiveCompatible(explosive: Explosive, scenario: any): boolean {
    // Check energy requirements
    const requiredYield = estimateRequiredYield(scenario);
    if (explosive.energyDensity * scenario.maxWeight < requiredYield.totalRequiredEnergy) {
        return false;
    }

    // Check density constraints
    const volumeConstraint = calculateAvailableVolume(scenario.radius, scenario.shape);
    const minRequiredDensity = scenario.maxWeight / volumeConstraint;
    if (explosive.density < minRequiredDensity) {
        return false;
    }

    // Check detonation velocity if targets specified
    if (scenario.targets) {
        const minRequiredVelocity = estimateRequiredDetonationVelocity(scenario);
        if (explosive.detonationVelocity < minRequiredVelocity) {
            return false;
        }
    }

    // Check critical diameter
    const minRequiredDiameter = estimateRequiredCriticalDiameter(scenario);
    if (explosive.criticalDiameter > minRequiredDiameter) {
        return false;
    }

    // Check stability for portable applications
    if (scenario.isPortable && explosive.stability < 3) {
        return false;
    }

    // Check temperature compatibility
    const standardTemp = 293; // K
    if (Math.abs(scenario.environment.temperature - standardTemp) > 50) {
        // Check thermal expansion and stability at extreme temperatures
        if (explosive.thermalExpansion > 1e-4 || explosive.stability < 4) {
            return false;
        }
    }

    // Check underwater compatibility
    if (scenario.environment.waterDepth !== undefined && scenario.environment.waterDepth > 0) {
        // Water-resistant explosives need oxygen balance near zero
        if (Math.abs(explosive.oxygenBalance) > 20) {
            return false;
        }
    }

    return true;
}

/**
 * Apply a parameter change to the scenario
 */
function applyParameterChange(scenario: any, parameterPath: string, newValue: any): void {
    // Handle array access notation like "targets[0].property"
    if (parameterPath.includes('[')) {
        const matches = parameterPath.match(/(\w+)\[(\d+)\]\.(.+)/);
        if (matches) {
            const [_, arrayName, indexStr, property] = matches;
            const index = parseInt(indexStr);

            if (scenario[arrayName] && Array.isArray(scenario[arrayName]) && scenario[arrayName][index]) {
                const propertyPath = property.split('.');
                let current = scenario[arrayName][index];

                // Navigate to the nested property
                for (let i = 0; i < propertyPath.length - 1; i++) {
                    if (!current[propertyPath[i]]) {
                        current[propertyPath[i]] = {};
                    }
                    current = current[propertyPath[i]];
                }

                // Set the value
                current[propertyPath[propertyPath.length - 1]] = newValue;
            }
            return;
        }
    }

    // Regular property paths like "shape.dimensions.radius"
    const parts = parameterPath.split('.');
    let current = scenario;

    // Navigate to the nested property
    for (let i = 0; i < parts.length - 1; i++) {
        if (!current[parts[i]]) {
            current[parts[i]] = {};
        }
        current = current[parts[i]];
    }

    // Set the value
    current[parts[parts.length - 1]] = newValue;
}

/**
 * Get the original parameter value from a scenario
 */
function originalParameterValue(scenario: any, parameterPath: string): any {
    // Handle array access notation
    if (parameterPath.includes('[')) {
        const matches = parameterPath.match(/(\w+)\[(\d+)\]\.(.+)/);
        if (matches) {
            const [_, arrayName, indexStr, property] = matches;
            const index = parseInt(indexStr);

            if (scenario[arrayName] && Array.isArray(scenario[arrayName]) && scenario[arrayName][index]) {
                const propertyPath = property.split('.');
                let current = scenario[arrayName][index];

                // Navigate to the nested property
                for (const part of propertyPath) {
                    if (current === undefined || current === null) return undefined;
                    current = current[part];
                }

                return current;
            }
            return undefined;
        }
    }

    // Regular property paths
    const parts = parameterPath.split('.');
    let current = scenario;

    // Navigate to the nested property
    for (const part of parts) {
        if (current === undefined || current === null) return undefined;
        current = current[part];
    }

    return current;
}

/**
 * Check if a scenario is physically consistent
 */
function isPhysicallyConsistent(scenario: any): boolean {
    // Validate basic physical constraints

    // 1. Volume must be positive and realistic
    if (scenario.radius <= 0) return false;

    // 2. Density and weight constraints
    const volume = calculateAvailableVolume(scenario.radius, scenario.shape);
    const effectiveDensity = scenario.maxWeight / volume;

    // Too high density is physically impossible
    if (effectiveDensity > 20) return false; // No known explosive exceeds 20 g/cm³

    // 3. Ensure shape dimensions are consistent
    if (scenario.shape.type === 'cylinder') {
        const height = scenario.shape.dimensions.height || (2 * scenario.radius);
        if (height <= 0) return false;
    } else if (scenario.shape.type === 'cube') {
        const side = scenario.shape.dimensions.side || (2 * scenario.radius);
        if (side <= 0) return false;
    } else if (scenario.shape.type === 'conical') {
        const height = scenario.shape.dimensions.height || (2 * scenario.radius);
        if (height <= 0) return false;
    } else if (scenario.shape.type === 'shaped_charge') {
        // Shaped charge validation
        if (scenario.shape.dimensions.coneAngle !== undefined) {
            const angle = scenario.shape.dimensions.coneAngle;
            // Practical limits for cone angles
            if (angle < 10 || angle > 120) return false;
        }

        if (scenario.shape.linerProperties) {
            // Liner thickness must be reasonable
            if (scenario.shape.linerProperties.thickness <= 0 ||
                scenario.shape.linerProperties.thickness > 10) return false;

            // Liner angle should match cone angle approximately
            if (scenario.shape.dimensions.coneAngle !== undefined &&
                scenario.shape.linerProperties.angle !== undefined) {
                const coneDiff = Math.abs(scenario.shape.dimensions.coneAngle - scenario.shape.linerProperties.angle);
                if (coneDiff > 30) return false; // Liner should roughly match cone
            }
        }
    }

    // 4. Environment checks
    if (scenario.environment.temperature <= 0) return false; // Temperature must be positive (Kelvin)
    if (scenario.environment.pressure <= 0) return false; // Pressure must be positive

    // 5. Casing physical consistency
    if (scenario.casing.thickness <= 0) return false;

    // 6. Target consistency
    if (scenario.targets) {
        for (const target of scenario.targets) {
            if (target.distance <= 0) return false; // Target must be at positive distance
            if (target.thickness < 0) return false; // Thickness cannot be negative
            if (target.crossSection <= 0) return false; // Cross-section must be positive
        }
    }

    return true;
}

/**
 * Evaluate parameter sensitivity for optimization
 */
function evaluateParameterSensitivity(
    originalScenario: any,
    modifiedScenario: any,
    parameter: string,
    explosives: Explosive[]
): number {
    // Clone scenarios for evaluation
    const baseScenario = JSON.parse(JSON.stringify(originalScenario));
    const newScenario = JSON.parse(JSON.stringify(modifiedScenario));

    // Calculate base metrics
    const baseMetrics = calculateFeasibilityMetrics(baseScenario, explosives);
    const newMetrics = calculateFeasibilityMetrics(newScenario, explosives);

    // Calculate normalized change in parameters
    const originalValue = originalParameterValue(baseScenario, parameter);
    const modifiedValue = originalParameterValue(newScenario, parameter);

    // Handle categorical variables differently
    if (typeof originalValue === 'string' || typeof modifiedValue === 'string') {
        return 1.0; // Fixed sensitivity for categorical changes
    }

    if (typeof originalValue !== 'number' || typeof modifiedValue !== 'number') {
        return 0.5; // Default for complex objects
    }

    // Calculate fractional parameter change
    const parameterChange = Math.abs(modifiedValue - originalValue) /
        Math.max(0.001, Math.abs(originalValue));

    // Calculate overall metrics change
    const energyChange = Math.abs(newMetrics.energyDeficit - baseMetrics.energyDeficit);
    const densityChange = Math.abs(newMetrics.densityDeficit - baseMetrics.densityDeficit);
    const velocityChange = Math.abs(newMetrics.velocityDeficit - baseMetrics.velocityDeficit);
    const diameterChange = Math.abs(newMetrics.diameterDeficit - baseMetrics.diameterDeficit);

    // Overall metrics delta (weighted sum)
    const metricsDelta = 3 * energyChange + 2 * densityChange + 2 * velocityChange + diameterChange;

    // Sensitivity is ratio of metrics change to parameter change
    return metricsDelta / Math.max(0.001, parameterChange);
}

/**
 * Describe sensitivity impact qualitatively
 */
function describeSensitivityImpact(sensitivity: number): string {
    if (sensitivity > 5) return "Very High";
    if (sensitivity > 2) return "High";
    if (sensitivity > 1) return "Moderate";
    if (sensitivity > 0.5) return "Low";
    return "Very Low";
}

/**
 * Apply advanced multi-parameter optimization techniques when simple parameter changes fail
 */
function applyAdvancedOptimization(
    originalScenario: any,
    existingChanges: Array<{ parameter: string; originalValue: any; newValue: any; unit: string }>,
    explosives: Explosive[],
    sensitivityTracking: Record<string, { sensitivity: number, impact: string }>
): {
    improved: boolean;
    newScenario: any;
    changes: Array<{ parameter: string; originalValue: any; newValue: any; unit: string }>;
} {
    // Clone the original scenario
    const newScenario = JSON.parse(JSON.stringify(originalScenario));
    const additionalChanges: Array<{ parameter: string; originalValue: any; newValue: any; unit: string }> = [];

    // Track parameters already changed
    const changedParameters = new Set(existingChanges.map(c => c.parameter));

    // ==== Strategy 1: Simultaneous multi-parameter adjustment ====
    // Find high-sensitivity parameters not yet changed
    const highSensitivityParams = Object.entries(sensitivityTracking)
        .filter(([param, data]) => data.sensitivity > 1.0 && !changedParameters.has(param))
        .sort((a, b) => b[1].sensitivity - a[1].sensitivity)
        .slice(0, 3); // Take top 3 most sensitive parameters

    if (highSensitivityParams.length > 0) {
        // For each parameter, calculate optimal adjustment
        for (const [param, data] of highSensitivityParams) {
            const originalValue = originalParameterValue(originalScenario, param);

            // Skip if value is not numerical
            if (typeof originalValue !== 'number') continue;

            // Generate adjustment options
            const paramOptions = getParameterAdjustmentOptions(param, originalValue, newScenario);

            // Try each option
            for (const option of paramOptions) {
                const testScenario = JSON.parse(JSON.stringify(newScenario));
                applyParameterChange(testScenario, param, option.value);

                if (isPhysicallyConsistent(testScenario) &&
                    improvesFeasibility(newScenario, testScenario, param, explosives)) {

                    // Apply successful change
                    applyParameterChange(newScenario, param, option.value);
                    additionalChanges.push({
                        parameter: getParameterName(param),
                        originalValue,
                        newValue: option.value,
                        unit: getUnitForProperty(param.split('.').pop() || '')
                    });

                    // Check if this made the scenario feasible
                    const result = optimizeExplosive(explosives, newScenario);
                    if (result) {
                        return {
                            improved: true,
                            newScenario,
                            changes: additionalChanges
                        };
                    }

                    break; // Move to next parameter
                }
            }
        }
    }

    // ==== Strategy 2: Shape adaptation based on target requirements ====
    if (newScenario.targets && newScenario.targets.length > 0 && !changedParameters.has('shape.type')) {
        const targetRequirements = analyzeTargetRequirements(newScenario.targets);

        // If penetration is primary need, consider shaped charge
        if (targetRequirements.penetrationPriority > 0.7 && newScenario.shape.type !== 'shaped_charge') {
            const originalShapeType = newScenario.shape.type;

            // Apply shaped charge conversion
            applyParameterChange(newScenario, 'shape.type', 'shaped_charge');
            applyParameterChange(newScenario, 'shape.dimensions.coneAngle', 60);

            // Add liner if missing
            if (!newScenario.shape.linerProperties) {
                applyParameterChange(newScenario, 'shape.linerProperties', {
                    material: 'copper',
                    density: 8.96,
                    thickness: 2.0,
                    angle: 42,
                    standoff: 2.0
                });
            }

            additionalChanges.push({
                parameter: 'Shape Type',
                originalValue: originalShapeType,
                newValue: 'shaped_charge',
                unit: ''
            });

            // Check if successful
            if (improvesFeasibility(originalScenario, newScenario, 'shape.type', explosives)) {
                const result = optimizeExplosive(explosives, newScenario);
                if (result) {
                    return {
                        improved: true,
                        newScenario,
                        changes: additionalChanges
                    };
                }
            }
        }

        // If blast is primary need, consider sphere or cylinder
        if (targetRequirements.blastPriority > 0.7 &&
            newScenario.shape.type !== 'sphere' &&
            newScenario.shape.type !== 'cylinder') {

            const originalShapeType = newScenario.shape.type;

            // Try sphere for omnidirectional blast
            applyParameterChange(newScenario, 'shape.type', 'sphere');

            additionalChanges.push({
                parameter: 'Shape Type',
                originalValue: originalShapeType,
                newValue: 'sphere',
                unit: ''
            });

            // Check if successful
            if (improvesFeasibility(originalScenario, newScenario, 'shape.type', explosives)) {
                const result = optimizeExplosive(explosives, newScenario);
                if (result) {
                    return {
                        improved: true,
                        newScenario,
                        changes: additionalChanges
                    };
                }
            }
        }
    }

    // ==== Strategy 3: Environmental parameter optimization ====
    if (newScenario.environment) {
        // Check if environment is extreme
        const standardTemp = 293; // K
        const standardPressure = 101.325; // kPa

        if (Math.abs(newScenario.environment.temperature - standardTemp) > 30 &&
            !changedParameters.has('environment.temperature')) {

            const originalTemp = newScenario.environment.temperature;
            const adjustedTemp = (originalTemp > standardTemp)
                ? Math.max(standardTemp, originalTemp - 20)
                : Math.min(standardTemp, originalTemp + 20);

            applyParameterChange(newScenario, 'environment.temperature', adjustedTemp);

            additionalChanges.push({
                parameter: 'Environmental Temperature',
                originalValue: originalTemp,
                newValue: adjustedTemp,
                unit: 'K'
            });
        }

        // Pressure optimization
        if (Math.abs(newScenario.environment.pressure - standardPressure) > 20 &&
            !changedParameters.has('environment.pressure')) {

            const originalPressure = newScenario.environment.pressure;
            const adjustedPressure = (originalPressure > standardPressure)
                ? Math.max(standardPressure, originalPressure - 10)
                : Math.min(standardPressure, originalPressure + 10);

            applyParameterChange(newScenario, 'environment.pressure', adjustedPressure);

            additionalChanges.push({
                parameter: 'Ambient Pressure',
                originalValue: originalPressure,
                newValue: adjustedPressure,
                unit: 'kPa'
            });
        }

        // Check if these changes improved feasibility
        if (additionalChanges.length > 0 && improvesFeasibility(originalScenario, newScenario, 'environment', explosives)) {
            const result = optimizeExplosive(explosives, newScenario);
            if (result) {
                return {
                    improved: true,
                    newScenario,
                    changes: additionalChanges
                };
            }
        }
    }

    // ==== Strategy 4: Material-driven optimization - adapt scenario to available materials ====
    const mostSuitableExplosive = findMostSuitableExplosive(explosives, originalScenario);
    if (mostSuitableExplosive) {
        // Adapt requirements to match this explosive's capabilities
        const adaptedChanges = adaptScenarioToExplosive(newScenario, mostSuitableExplosive, changedParameters);

        if (adaptedChanges.length > 0) {
            // Apply all changes
            for (const change of adaptedChanges) {
                applyParameterChange(newScenario, change.paramPath, change.newValue);

                additionalChanges.push({
                    parameter: change.paramName,
                    originalValue: change.oldValue,
                    newValue: change.newValue,
                    unit: change.unit
                });
            }

            // Check if successful
            const result = optimizeExplosive(explosives, newScenario);
            if (result) {
                return {
                    improved: true,
                    newScenario,
                    changes: additionalChanges
                };
            }
        }
    }

    // If nothing worked
    return {
        improved: additionalChanges.length > 0,
        newScenario,
        changes: additionalChanges
    };
}

/**
 * Get parameter name from parameter path
 */
function getParameterName(paramPath: string): string {
    // Handle array notation
    if (paramPath.includes('[')) {
        const matches = paramPath.match(/(\w+)\[(\d+)\]\.(.+)/);
        if (matches) {
            const [_, arrayName, indexStr, property] = matches;
            return `${arrayName.charAt(0).toUpperCase() + arrayName.slice(1)} ${parseInt(indexStr) + 1} ${property.split('.').map(p => p.charAt(0).toUpperCase() + p.slice(1)).join(' ')}`;
        }
    }

    // Handle regular path
    return paramPath.split('.').map(p => p.charAt(0).toUpperCase() + p.slice(1)).join(' ');
}

/**
 * Generate parameter adjustment options for optimization
 */
function getParameterAdjustmentOptions(
    paramPath: string,
    currentValue: number,
    scenario: any
): Array<{ value: number; priority: number }> {
    const options: Array<{ value: number; priority: number }> = [];
    const paramName = paramPath.split('.').pop() || '';

    // Generate options based on parameter type
    switch (paramName) {
        case 'radius':
            // Try increasing radius
            options.push({ value: currentValue * 1.2, priority: 0.9 });
            options.push({ value: currentValue * 1.5, priority: 0.7 });
            options.push({ value: currentValue * 0.8, priority: 0.5 }); // Sometimes smaller is better
            break;

        case 'maxWeight':
            // Try increasing weight allowance
            options.push({ value: currentValue * 1.25, priority: 0.95 });
            options.push({ value: currentValue * 1.5, priority: 0.8 });
            options.push({ value: currentValue * 2.0, priority: 0.6 });
            break;

        case 'standoff':
            // Get optimal standoff if available
            if (scenario.shape.type === 'shaped_charge' && scenario.shape.linerProperties) {
                const optimalStandoff = ExplosiveModel.calculateOptimalStandoff(scenario.shape.linerProperties);
                options.push({ value: optimalStandoff, priority: 0.95 });
                options.push({ value: optimalStandoff * 0.9, priority: 0.7 });
                options.push({ value: optimalStandoff * 1.1, priority: 0.7 });
            } else {
                options.push({ value: currentValue * 1.5, priority: 0.8 });
                options.push({ value: currentValue * 0.7, priority: 0.7 });
            }
            break;

        case 'coneAngle':
            // Standard angles for shaped charges
            options.push({ value: 42, priority: 0.9 }); // Optimal for many applications
            options.push({ value: 60, priority: 0.8 });
            options.push({ value: 30, priority: 0.7 });
            break;

        case 'thickness': // For casing thickness
            if (paramPath.includes('casing')) {
                // Calculate optimal casing based on physics
                const casingDensity = scenario.casing?.density || 7.8; // Default to steel
                const explosiveDensity = 1.6; // Typical explosive density

                // Gurney equation optimal M/C ratio for fragment velocity
                const optimalMassRatio = 0.6; // Mass of casing / Mass of explosive

                // Rough estimate based on cylinder approximation
                const chargeRadius = scenario.radius - currentValue;
                const optimalThickness = (optimalMassRatio * explosiveDensity * chargeRadius) / casingDensity;

                options.push({ value: optimalThickness, priority: 0.9 });
                options.push({ value: currentValue * 0.8, priority: 0.8 }); // Thinner
                options.push({ value: currentValue * 1.2, priority: 0.7 }); // Thicker
            } else if (paramPath.includes('linerProperties')) {
                // For shaped charge liners
                options.push({ value: 2.0, priority: 0.9 }); // Standard thickness
                options.push({ value: 1.5, priority: 0.8 });
                options.push({ value: 3.0, priority: 0.7 });
            }
            break;

        case 'temperature':
            // Move toward standard temperature
            const standardTemp = 293; // K
            options.push({ value: standardTemp, priority: 0.9 });
            options.push({
                value: currentValue + (standardTemp - currentValue) * 0.5,
                priority: 0.8
            }); // Move halfway to standard
            break;

        case 'pressure':
            // Move toward standard pressure
            const standardPressure = 101.325; // kPa
            options.push({ value: standardPressure, priority: 0.9 });
            options.push({
                value: currentValue + (standardPressure - currentValue) * 0.5,
                priority: 0.8
            }); // Move halfway to standard
            break;

        case 'altitude':
            // Reduce altitude
            options.push({ value: currentValue * 0.7, priority: 0.9 });
            options.push({ value: currentValue * 0.5, priority: 0.8 });
            options.push({ value: 0, priority: 0.7 }); // Sea level
            break;

        case 'waterDepth':
            // Reduce water depth
            options.push({ value: currentValue * 0.7, priority: 0.9 });
            options.push({ value: currentValue * 0.5, priority: 0.8 });
            options.push({ value: 1, priority: 0.7 }); // Shallow water
            break;

        case 'criticalDamageThreshold':
            // Reduce damage threshold requirements
            options.push({ value: currentValue * 0.8, priority: 0.9 });
            options.push({ value: currentValue * 0.6, priority: 0.7 });
            options.push({ value: currentValue * 0.5, priority: 0.5 });
            break;

        default:
            // Generic adjustments for numeric parameters
            options.push({ value: currentValue * 1.2, priority: 0.8 });
            options.push({ value: currentValue * 0.8, priority: 0.7 });
            options.push({ value: currentValue * 1.5, priority: 0.6 });
            options.push({ value: currentValue * 0.5, priority: 0.5 });
    }

    // Sort by priority
    return options.sort((a, b) => b.priority - a.priority);
}

/**
 * Analyze target requirements to determine optimal shape and configuration
 */
function analyzeTargetRequirements(targets: Target[]): {
    penetrationPriority: number;
    blastPriority: number;
    fragmentationPriority: number;
    precisionPriority: number;
} {
    let penetrationScore = 0;
    let blastScore = 0;
    let fragmentationScore = 0;
    let precisionScore = 0;

    for (const target of targets) {
        // Penetration targets are characterized by high density materials and thickness
        if (target.material === 'steel' || target.material === 'concrete' || target.material === 'armor') {
            penetrationScore += 1.0;

            // Thicker materials increase penetration priority
            if (target.thickness > 5) { // cm
                penetrationScore += 0.5;
            }

            // Small cross-section requires precision
            if (target.crossSection < 1) { // m²
                precisionScore += 1.0;
            }
        }

        // Blast targets are typically structures or human targets 
        if (target.material === 'structure' || target.material === 'human') {
            blastScore += 1.0;

            // High damage threshold needs strong blast
            if (target.criticalDamageThreshold > 200) { // kPa
                blastScore += 0.5;
            }

            // Large cross-section indicates area target
            if (target.crossSection > 10) { // m²
                blastScore += 0.3;
                fragmentationScore += 0.5; // Fragmentation effective against area targets
            }
        }

        // Fragmentation priority increases with distance
        if (target.distance > 20) { // m
            fragmentationScore += 0.3;
        }

        // Target with shielding requires penetration or precision
        if (target.shielding && target.shielding > 0.3) {
            penetrationScore += 0.5;
            precisionScore += 0.3;
        }
    }

    // Normalize scores
    const totalTargets = targets.length;
    penetrationScore /= totalTargets;
    blastScore /= totalTargets;
    fragmentationScore /= totalTargets;
    precisionScore /= totalTargets;

    return {
        penetrationPriority: Math.min(1.0, penetrationScore),
        blastPriority: Math.min(1.0, blastScore),
        fragmentationPriority: Math.min(1.0, fragmentationScore),
        precisionPriority: Math.min(1.0, precisionScore)
    };
}

/**
 * Find the most suitable explosive from available options
 */
function findMostSuitableExplosive(explosives: Explosive[], scenario: any): Explosive | null {
    if (explosives.length === 0) return null;

    // Calculate scenario requirements
    const requiredYield = estimateRequiredYield(scenario);
    const requiredDetonationVelocity = estimateRequiredDetonationVelocity(scenario);
    const requiredCriticalDiameter = estimateRequiredCriticalDiameter(scenario);

    // Score each explosive
    const scoredExplosives = explosives.map(explosive => {
        let score = 0;

        // Energy yield score (0-40 points)
        const yieldScore = Math.min(1.0, (explosive.energyDensity * scenario.maxWeight) / requiredYield.totalRequiredEnergy);
        score += yieldScore * 40;

        // Detonation velocity score (0-30 points)
        const velocityScore = Math.min(1.0, explosive.detonationVelocity / requiredDetonationVelocity);
        score += velocityScore * 30;

        // Critical diameter score (0-10 points)
        const diameterRatio = requiredCriticalDiameter / Math.max(1, explosive.criticalDiameter);
        const diameterScore = Math.min(1.0, diameterRatio);
        score += diameterScore * 10;

        // Stability score for portable scenarios (0-10 points)
        if (scenario.isPortable) {
            const stabilityScore = Math.min(1.0, explosive.stability / 5);
            score += stabilityScore * 10;
        } else {
            // For non-portable, prioritize power
            score += 5; // Half the points automatically
        }

        // Environmental compatibility (0-10 points)
        const envScore = calculateEnvironmentalCompatibility(explosive, scenario.environment);
        score += envScore * 10;

        return { explosive, score };
    });

    // Sort by score and return best match
    scoredExplosives.sort((a, b) => b.score - a.score);
    return scoredExplosives[0]?.explosive || null;
}

/**
 * Calculate environmental compatibility score (0-1)
 */
function calculateEnvironmentalCompatibility(explosive: Explosive, environment: Environment): number {
    let score = 1.0;

    // Temperature extremes
    const standardTemp = 293; // K
    if (Math.abs(environment.temperature - standardTemp) > 50) {
        // Cold affects initiation
        if (environment.temperature < standardTemp) {
            score -= 0.2 * (1 - explosive.stability / 5);
        }
        // Heat affects stability
        else {
            score -= 0.3 * (1 - explosive.stability / 5);

            // High thermal expansion compounds the issue
            if (explosive.thermalExpansion > 1e-4) {
                score -= 0.1;
            }
        }
    }

    // Underwater applications
    if (environment.waterDepth !== undefined && environment.waterDepth > 0) {
        // Water sensitivity depends on oxygen balance
        score -= 0.2 * Math.min(1.0, Math.abs(explosive.oxygenBalance) / 50);
    }

    // High altitude effects
    if (environment.altitude > 2000) {
        // Reduced pressure affects detonation
        score -= 0.1 * Math.min(1.0, environment.altitude / 8000);
    }

    // Humidity effects
    if (environment.humidity > 80) {
        // Hygroscopic materials affected by high humidity
        if (Math.abs(explosive.oxygenBalance) > 25) {
            score -= 0.15;
        }
    }

    return Math.max(0, score);
}

/**
 * Adapt scenario to available explosive properties
 */
function adaptScenarioToExplosive(
    scenario: any,
    explosive: Explosive,
    changedParameters: Set<string>
): Array<{
    paramPath: string;
    paramName: string;
    oldValue: any;
    newValue: any;
    unit: string;
}> {
    const changes: Array<{
        paramPath: string;
        paramName: string;
        oldValue: any;
        newValue: any;
        unit: string;
    }> = [];

    // Required yield for the scenario
    const requiredYield = estimateRequiredYield(scenario);

    // Check if weight needs to be increased
    if (!changedParameters.has('maxWeight') &&
        explosive.energyDensity * scenario.maxWeight < requiredYield.totalRequiredEnergy) {

        const recommendedWeight = (requiredYield.totalRequiredEnergy / explosive.energyDensity) * 1.1;

        changes.push({
            paramPath: 'maxWeight',
            paramName: 'Maximum Weight',
            oldValue: scenario.maxWeight,
            newValue: recommendedWeight,
            unit: 'kg'
        });
    }

    // Check if critical diameter constraint requires radius increase
    const criticalDiameterMm = explosive.criticalDiameter;
    const scenarioDiameterMm = 2 * scenario.radius * 10; // cm to mm

    if (!changedParameters.has('radius') && criticalDiameterMm > scenarioDiameterMm) {
        const recommendedRadius = (criticalDiameterMm / 20) * 1.1; // mm to cm with 10% margin

        changes.push({
            paramPath: 'radius',
            paramName: 'Radius',
            oldValue: scenario.radius,
            newValue: recommendedRadius,
            unit: 'cm'
        });
    }

    // Adapt casing for explosive properties if not already changed
    if (!changedParameters.has('casing.thickness')) {
        // Calculate optimal casing thickness for this explosive
        const optimalThickness = calculateOptimalCasingThickness(
            explosive.detonationVelocity,
            explosive.energyDensity,
            scenario.radius,
            scenario.casing.density
        );

        if (Math.abs(scenario.casing.thickness - optimalThickness) / optimalThickness > 0.3) {
            changes.push({
                paramPath: 'casing.thickness',
                paramName: 'Casing Thickness',
                oldValue: scenario.casing.thickness,
                newValue: optimalThickness,
                unit: 'cm'
            });
        }
    }

    // If shaped charge, adapt standoff distance to explosive properties
    if (scenario.shape.type === 'shaped_charge' &&
        scenario.shape.linerProperties &&
        !changedParameters.has('shape.dimensions.standoff')) {

        // Recalculate optimal standoff for this specific explosive
        const optimalStandoff = calculateOptimalStandoffForExplosive(
            explosive,
            scenario.shape.linerProperties
        );

        if (Math.abs(scenario.shape.dimensions.standoff - optimalStandoff) / optimalStandoff > 0.2) {
            changes.push({
                paramPath: 'shape.dimensions.standoff',
                paramName: 'Standoff Distance',
                oldValue: scenario.shape.dimensions.standoff,
                newValue: optimalStandoff,
                unit: 'cm'
            });
        }
    }

    // If stability issues and portable, suggest non-portable
    if (scenario.isPortable &&
        explosive.stability < 3 &&
        !changedParameters.has('isPortable')) {

        changes.push({
            paramPath: 'isPortable',
            paramName: 'Portability Requirement',
            oldValue: true,
            newValue: false,
            unit: ''
        });
    }

    return changes;
}

/**
 * Calculate optimal standoff for shaped charge with specific explosive
 */
function calculateOptimalStandoffForExplosive(explosive: Explosive, liner: Liner): number {
    // Base calculation from ExplosiveModel
    let baseStandoff = ExplosiveModel.calculateOptimalStandoff(liner);

    // Adjust based on explosive properties
    // Detonation velocity affects jet formation dynamics
    const velocityFactor = explosive.detonationVelocity / 8000; // Normalized to typical PETN/RDX

    // Higher detonation pressure creates more coherent jets that work at longer standoffs
    const pressureFactor = explosive.chapmanJouguetPressure / 25; // Normalized to typical pressure

    // High energy density explosives may benefit from slightly increased standoff
    const energyFactor = explosive.energyDensity / 5; // Normalized to typical energy 

    // Combined adjustment factor (weighting based on importance)
    const adjustmentFactor = 0.5 * velocityFactor + 0.3 * pressureFactor + 0.2 * energyFactor;

    // Apply adjustment with limits
    return baseStandoff * Math.max(0.8, Math.min(1.2, adjustmentFactor));
}

/**
 * Find the closest feasible scenario when all optimizations fail
 */
function findClosestFeasibleScenario(
    originalScenario: any,
    explosives: Explosive[]
): {
    distance: number;
    scenario: any;
    material: Explosive;
    changes: Array<{ parameter: string; originalValue: any; newValue: any; unit: string }>;
} {
    let minDistance = Number.MAX_VALUE;
    let bestScenario: any = null;
    let bestMaterial: Explosive | null = null;
    let bestChanges: Array<{ parameter: string; originalValue: any; newValue: any; unit: string }> = [];

    // Define parameter spaces to search
    const parameterSpaces = generateParameterSpaces(originalScenario);

    // Maximum number of iterations to prevent excessive computation
    const MAX_ITERATIONS = 500;
    let iterations = 0;

    // Use simulated annealing for global search
    let currentScenario = JSON.parse(JSON.stringify(originalScenario));
    let currentDistance = calculateScenarioDistance(currentScenario, originalScenario);
    let temperature = 100.0; // Initial temperature
    const coolingRate = 0.97; // Cooling rate

    // Track changes for best solution
    let currentChanges: Array<{ parameter: string; originalValue: any; newValue: any; unit: string }> = [];

    while (temperature > 0.1 && iterations < MAX_ITERATIONS) {
        iterations++;

        // Generate neighbor by modifying parameters
        const neighborScenario = JSON.parse(JSON.stringify(currentScenario));
        const parameterChanges = generateNeighborScenario(neighborScenario, parameterSpaces);

        // Skip if not physically consistent
        if (!isPhysicallyConsistent(neighborScenario)) {
            continue;
        }

        // Check if neighbor is feasible
        const feasibilityResult = optimizeExplosive(explosives, neighborScenario);
        const neighborDistance = calculateScenarioDistance(neighborScenario, originalScenario);

        if (feasibilityResult) {
            // Found feasible solution - check if it's closer than previous best
            if (neighborDistance < minDistance) {
                minDistance = neighborDistance;
                bestScenario = JSON.parse(JSON.stringify(neighborScenario));
                bestMaterial = feasibilityResult.material;
                bestChanges = parameterChanges.map(change => ({
                    parameter: getParameterName(change.parameter),
                    originalValue: originalParameterValue(originalScenario, change.parameter),
                    newValue: change.value,
                    unit: getUnitForProperty(change.parameter.split('.').pop() || '')
                }));
            }

            // We could stop here, but continuing search might find even better solutions
            // Adjust temperature to encourage exploration of feasible regions
            temperature *= 0.9;
        } else {
            // Not feasible, but decide whether to move to this state anyway
            const deltaDistance = neighborDistance - currentDistance;

            // Accept worse solutions with probability based on temperature and distance difference
            if (deltaDistance < 0 || Math.random() < Math.exp(-deltaDistance / temperature)) {
                currentScenario = neighborScenario;
                currentDistance = neighborDistance;
                // Convert parameterChanges to the expected format
                currentChanges = parameterChanges.map(change => ({
                    parameter: getParameterName(change.parameter),
                    originalValue: originalParameterValue(originalScenario, change.parameter),
                    newValue: change.value,
                    unit: getUnitForProperty(change.parameter.split('.').pop() || '')
                }));
            }

            // Cool down the temperature
            temperature *= coolingRate;
        }
    }

    return {
        distance: minDistance,
        scenario: bestScenario || originalScenario,
        material: bestMaterial || explosives[0],
        changes: bestChanges
    };
}

/**
 * Generate parameter spaces for optimization search
 */
function generateParameterSpaces(scenario: any): Record<string, {
    min: number;
    max: number;
    step: number;
    type: 'numeric' | 'categorical';
    categories?: any[];
}> {
    const spaces: Record<string, any> = {};

    // Radius search space
    spaces['radius'] = {
        min: scenario.radius * 0.5,
        max: scenario.radius * 2.0,
        step: scenario.radius * 0.1,
        type: 'numeric'
    };

    // Weight search space
    spaces['maxWeight'] = {
        min: scenario.maxWeight * 0.5,
        max: scenario.maxWeight * 3.0,
        step: scenario.maxWeight * 0.2,
        type: 'numeric'
    };

    // Shape type
    spaces['shape.type'] = {
        type: 'categorical',
        categories: ['sphere', 'cylinder', 'cube', 'hemisphere', 'conical', 'shaped_charge']
    };

    // Casing thickness
    spaces['casing.thickness'] = {
        min: Math.max(0.1, scenario.casing.thickness * 0.5),
        max: scenario.casing.thickness * 2.0,
        step: scenario.casing.thickness * 0.1,
        type: 'numeric'
    };

    // Environmental parameters
    if (scenario.environment) {
        // Temperature
        spaces['environment.temperature'] = {
            min: Math.max(223, scenario.environment.temperature * 0.8), // Min 223K (-50°C)
            max: Math.min(373, scenario.environment.temperature * 1.2), // Max 373K (100°C)
            step: 10,
            type: 'numeric'
        };

        // Altitude
        if (scenario.environment.altitude !== undefined) {
            spaces['environment.altitude'] = {
                min: 0,
                max: Math.max(5000, scenario.environment.altitude * 1.5),
                step: 500,
                type: 'numeric'
            };
        }

        // Water depth
        if (scenario.environment.waterDepth !== undefined && scenario.environment.waterDepth > 0) {
            spaces['environment.waterDepth'] = {
                min: 0.5,
                max: Math.max(100, scenario.environment.waterDepth * 1.5),
                step: Math.max(0.5, scenario.environment.waterDepth * 0.1),
                type: 'numeric'
            };
        }
    }

    // If shaped charge, add specific parameters
    if (scenario.shape.type === 'shaped_charge') {
        // Cone angle
        spaces['shape.dimensions.coneAngle'] = {
            min: 30,
            max: 90,
            step: 5,
            type: 'numeric'
        };

        // Standoff distance
        spaces['shape.dimensions.standoff'] = {
            min: 0.5,
            max: 10,
            step: 0.5,
            type: 'numeric'
        };

        // Liner properties if available
        if (scenario.shape.linerProperties) {
            spaces['shape.linerProperties.thickness'] = {
                min: 0.5,
                max: 5,
                step: 0.5,
                type: 'numeric'
            };

            spaces['shape.linerProperties.angle'] = {
                min: 20,
                max: 60,
                step: 5,
                type: 'numeric'
            };
        }
    }

    // Portability as a binary choice
    spaces['isPortable'] = {
        type: 'categorical',
        categories: [true, false]
    };

    return spaces;
}

/**
 * Generate a neighbor scenario by modifying random parameters
 */
function generateNeighborScenario(
    scenario: any,
    parameterSpaces: Record<string, any>
): Array<{ parameter: string; value: any }> {
    // Choose parameters to modify (1-3 random parameters)
    const parameterKeys = Object.keys(parameterSpaces);
    const modificationCount = 1 + Math.floor(Math.random() * 2); // 1-3 modifications
    const selectedParameters: Array<string> = [];

    // Randomly select parameters to modify
    for (let i = 0; i < modificationCount; i++) {
        const randomIndex = Math.floor(Math.random() * parameterKeys.length);
        const param = parameterKeys[randomIndex];

        // Avoid duplicates
        if (!selectedParameters.includes(param)) {
            selectedParameters.push(param);
        }
    }

    // Apply modifications
    const changes: Array<{ parameter: string; value: any }> = [];

    for (const param of selectedParameters) {
        const space = parameterSpaces[param];
        let newValue;

        if (space.type === 'numeric') {
            // Generate random value within range
            const steps = Math.floor((space.max - space.min) / space.step);
            const randomSteps = Math.floor(Math.random() * steps);
            newValue = space.min + randomSteps * space.step;
        } else if (space.type === 'categorical') {
            // Choose random category
            const randomIndex = Math.floor(Math.random() * space.categories.length);
            newValue = space.categories[randomIndex];
        }

        // Apply the parameter change
        applyParameterChange(scenario, param, newValue);
        changes.push({ parameter: param, value: newValue });
    }

    return changes;
}

/**
 * Calculate distance between two scenarios (for measuring how much the scenario has changed)
 */
function calculateScenarioDistance(scenarioA: any, scenarioB: any): number {
    let distance = 0;

    // Compare radius (normalized)
    if (scenarioA.radius !== undefined && scenarioB.radius !== undefined) {
        const radiusChange = Math.abs(scenarioA.radius - scenarioB.radius) / scenarioB.radius;
        distance += radiusChange * 2; // Weight radius changes heavily
    }

    // Compare weight (normalized)
    if (scenarioA.maxWeight !== undefined && scenarioB.maxWeight !== undefined) {
        const weightChange = Math.abs(scenarioA.maxWeight - scenarioB.maxWeight) / scenarioB.maxWeight;
        distance += weightChange * 2; // Weight changes heavily
    }

    // Compare shape type
    if (scenarioA.shape?.type !== scenarioB.shape?.type) {
        distance += 1.0; // Significant shape change
    }

    // Compare casing thickness
    if (scenarioA.casing?.thickness !== undefined && scenarioB.casing?.thickness !== undefined) {
        const thicknessChange = Math.abs(scenarioA.casing.thickness - scenarioB.casing.thickness) /
            scenarioB.casing.thickness;
        distance += thicknessChange;
    }

    // Compare environmental parameters
    if (scenarioA.environment && scenarioB.environment) {
        // Temperature
        if (scenarioA.environment.temperature !== undefined && scenarioB.environment.temperature !== undefined) {
            const tempChange = Math.abs(scenarioA.environment.temperature - scenarioB.environment.temperature) / 100;
            distance += tempChange * 0.5;
        }

        // Altitude
        if (scenarioA.environment.altitude !== undefined && scenarioB.environment.altitude !== undefined) {
            const altChange = Math.abs(scenarioA.environment.altitude - scenarioB.environment.altitude) / 5000;
            distance += altChange * 0.5;
        }

        // Water depth
        if (scenarioA.environment.waterDepth !== undefined && scenarioB.environment.waterDepth !== undefined) {
            const depthChange = Math.abs(scenarioA.environment.waterDepth - scenarioB.environment.waterDepth) /
                Math.max(1, scenarioB.environment.waterDepth);
            distance += depthChange * 0.7;
        }
    }

    // Compare shaped charge parameters if applicable
    if (scenarioA.shape?.type === 'shaped_charge' && scenarioB.shape?.type === 'shaped_charge') {
        // Cone angle
        if (scenarioA.shape.dimensions?.coneAngle !== undefined &&
            scenarioB.shape.dimensions?.coneAngle !== undefined) {
            const angleChange = Math.abs(scenarioA.shape.dimensions.coneAngle -
                scenarioB.shape.dimensions.coneAngle) / 60;
            distance += angleChange * 0.7;
        }

        // Standoff
        if (scenarioA.shape.dimensions?.standoff !== undefined &&
            scenarioB.shape.dimensions?.standoff !== undefined) {
            const standoffChange = Math.abs(scenarioA.shape.dimensions.standoff -
                scenarioB.shape.dimensions.standoff) /
                Math.max(1, scenarioB.shape.dimensions.standoff);
            distance += standoffChange * 0.7;
        }
    }

    // Portability
    if (scenarioA.isPortable !== scenarioB.isPortable) {
        distance += 0.8;
    }

    return distance;
}

/**
 * Calculate confidence level for recommended changes
 */
function calculateConfidence(
    changes: Array<{ parameter: string; originalValue: any; newValue: any; unit: string }>,
    originalScenario: any,
    modifiedScenario: any
): number {
    // Base confidence starts high and gets reduced
    let confidenceScore = 95;

    // Distance between scenarios reduces confidence
    const scenarioDistance = calculateScenarioDistance(originalScenario, modifiedScenario);
    confidenceScore -= scenarioDistance * 20; // Higher distance = lower confidence

    // Number of changes reduces confidence
    confidenceScore -= changes.length * 3; // Each change reduces confidence

    // Drastic changes reduce confidence more
    for (const change of changes) {
        if (typeof change.originalValue === 'number' && typeof change.newValue === 'number') {
            const relativeChange = Math.abs(change.newValue - change.originalValue) /
                Math.max(0.001, Math.abs(change.originalValue));

            // Large changes (more than doubling or halving) reduce confidence
            if (relativeChange > 1.0) {
                confidenceScore -= 5 + (relativeChange - 1.0) * 10;
            }
        } else if (change.originalValue !== change.newValue) {
            // Categorical changes
            confidenceScore -= 8;
        }
    }

    // Critical parameter changes affect confidence more
    const criticalParameterChanges = changes.filter(c =>
        c.parameter.includes('Weight') ||
        c.parameter.includes('Radius') ||
        c.parameter.includes('Shape Type')
    );

    confidenceScore -= criticalParameterChanges.length * 5;

    // Cap confidence between 0-100
    return Math.max(0, Math.min(100, confidenceScore));
}

// Example targets database
const targets: Target[] = [
    {
        material: 'human',
        thickness: 30, // cm (average torso)
        density: 1.0, // g/cm³ (approximately)
        distance: 50, // m
        crossSection: 0.75, // m²
        criticalDamageThreshold: 35 // kPa (lung damage threshold)
    },
    {
        material: 'reinforced concrete',
        thickness: 30, // cm
        density: 2.4, // g/cm³
        distance: 20, // m
        crossSection: 10, // m²
        shielding: 0,
        criticalDamageThreshold: 70 // kPa
    },
    {
        material: 'vehicle',
        thickness: 0.5, // cm (sheet metal)
        density: 7.85, // g/cm³
        distance: 30, // m
        crossSection: 8, // m²
        criticalDamageThreshold: 30 // kPa (window breakage)
    }
];

// Updated scenario examples with comprehensive physical modeling
// Scenario 1: Portable Demolition Charge
runComprehensiveScenario(
    'Portable Demolition Charge',
    50, // meters
    { type: 'cylinder', dimensions: { radius: 5, height: 20 } }, // cm
    3, // kg
    true,
    {
        material: 'aluminum',
        density: 2.7,
        thickness: 0.3,
        yieldStrength: 240,
        ultimateStrength: 310,
        elasticModulus: 70,
        hardness: 75,
        elongation: 12,
        thermalConductivity: 167,
        meltingPoint: 933,
        fractureToughness: 24,
        poissonRatio: 0.33,
        heatCapacity: 900,
        corrosionResistance: 7,
        fragmentationPattern: 'natural'
    },
    environments[0], // Standard atmosphere
    [targets[1]] // Reinforced concrete target
);

// Scenario 2: Static Mining Explosive
runComprehensiveScenario(
    'Static Mining Explosive',
    200, // meters
    { type: 'cube', dimensions: { side: 50 } }, // cm
    100, // kg
    false,
    {
        material: 'steel',
        density: 7.85,
        thickness: 1,
        yieldStrength: 250,
        ultimateStrength: 420,
        elasticModulus: 200,
        hardness: 130,
        elongation: 25,
        thermalConductivity: 45,
        meltingPoint: 1723,
        fractureToughness: 50,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 4,
        fragmentationPattern: 'natural'
    },
    {
        ...environments[0], // Standard atmosphere
        soilType: 'rock',
        soilDensity: 2.7 // g/cm³ (typical for granite)
    }
);

// Scenario 3: Shaped Charge
runComprehensiveScenario(
    'Shaped Charge',
    10, // meters
    {
        type: 'shaped_charge',
        dimensions: {
            radius: 5, // cm
            height: 15, // cm
            standoff: 10 // cm
        },
        linerProperties: {
            material: 'copper',
            density: 8.96,
            thickness: 3, // mm
            angle: 42, // degrees
            standoff: 3 // calibers
        }
    },
    2, // kg
    true,
    {
        material: 'steel (high strength)',
        density: 7.85,
        thickness: 0.4,
        yieldStrength: 550,
        ultimateStrength: 800,
        elasticModulus: 210,
        hardness: 280,
        elongation: 15,
        thermalConductivity: 40,
        meltingPoint: 1723,
        fractureToughness: 80,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 5,
        fragmentationPattern: 'controlled'
    },
    environments[0], // Standard atmosphere
    [targets[2]] // Vehicle target
);

// Underwater charge scenario
runComprehensiveScenario(
    'Underwater Charge',
    30, // meters
    { type: 'sphere', dimensions: { radius: 10 } }, // cm
    15, // kg
    false,
    {
        material: 'titanium alloy',
        density: 4.5,
        thickness: 0.8,
        yieldStrength: 880,
        ultimateStrength: 950,
        elasticModulus: 110,
        hardness: 320,
        elongation: 14,
        thermalConductivity: 7.2,
        meltingPoint: 1941,
        fractureToughness: 100,
        poissonRatio: 0.34,
        heatCapacity: 530,
        corrosionResistance: 9,
        fragmentationPattern: 'controlled'
    },
    environments[3] // Underwater environment
);

runComprehensiveScenario(
    'Baseline Combustion Test with TNT',
    10, // meters
    { type: 'cylinder', dimensions: { radius: 5, height: 20 } }, // cm
    5, // kg
    false,
    {
        material: 'Steel (mild)',
        density: 7.85,
        thickness: 0.5,
        yieldStrength: 250,
        ultimateStrength: 420,
        elasticModulus: 200,
        hardness: 130,
        elongation: 25,
        thermalConductivity: 45,
        meltingPoint: 1723,
        fractureToughness: 50,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 4,
        fragmentationPattern: 'natural'
    },
    environments[0] // Standard atmosphere
);

runComprehensiveScenario(
    'Shockwave Redirection with HMX and Ribbed Containment',
    10, // meters
    { type: 'cylinder', dimensions: { radius: 5, height: 20 }, contours: 'ribbed' }, // cm
    5, // kg
    false,
    {
        material: 'Steel (high strength)',
        density: 7.85,
        thickness: 0.8,
        yieldStrength: 550,
        ultimateStrength: 800,
        elasticModulus: 210,
        hardness: 280,
        elongation: 15,
        thermalConductivity: 40,
        meltingPoint: 1723,
        fractureToughness: 80,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 5,
        fragmentationPattern: 'controlled'
    },
    environments[0], // Standard atmosphere
    [targets[1]] // Reinforced concrete target
);

runComprehensiveScenario(
    'Energy Absorption with C4 and Titanium Casing',
    50, // meters
    { type: 'sphere', dimensions: { radius: 10 } }, // cm
    3, // kg
    true,
    {
        material: 'Titanium Alloy',
        density: 4.5,
        thickness: 0.6,
        yieldStrength: 880,
        ultimateStrength: 950,
        elasticModulus: 110,
        hardness: 320,
        elongation: 14,
        thermalConductivity: 7.2,
        meltingPoint: 1941,
        fractureToughness: 100,
        poissonRatio: 0.34,
        heatCapacity: 530,
        corrosionResistance: 9,
        fragmentationPattern: 'controlled'
    },
    environments[0], // Standard atmosphere
    [targets[0]] // Human target
);

runComprehensiveScenario(
    'Underwater Blast Containment with ANFO',
    30, // meters
    { type: 'sphere', dimensions: { radius: 15 } }, // cm
    10, // kg
    false,
    {
        material: 'Titanium Alloy',
        density: 4.5,
        thickness: 1.0,
        yieldStrength: 880,
        ultimateStrength: 950,
        elasticModulus: 110,
        hardness: 320,
        elongation: 14,
        thermalConductivity: 7.2,
        meltingPoint: 1941,
        fractureToughness: 100,
        poissonRatio: 0.34,
        heatCapacity: 530,
        corrosionResistance: 9,
        fragmentationPattern: 'controlled'
    },
    environments[3] // Underwater environment
);

runComprehensiveScenario(
    'Shaped Charge for Precision Breaching with PETN',
    10, // meters
    {
        type: 'shaped_charge',
        dimensions: { radius: 5, height: 15, standoff: 10 }, // cm
        linerProperties: { material: 'copper', density: 8.96, thickness: 3, angle: 42, standoff: 3 }
    },
    2, // kg
    true,
    {
        material: 'Steel (high strength)',
        density: 7.85,
        thickness: 0.4,
        yieldStrength: 550,
        ultimateStrength: 800,
        elasticModulus: 210,
        hardness: 280,
        elongation: 15,
        thermalConductivity: 40,
        meltingPoint: 1723,
        fractureToughness: 80,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 5,
        fragmentationPattern: 'controlled'
    },
    environments[0], // Standard atmosphere
    [targets[1]] // Reinforced concrete target
);

runComprehensiveScenario(
    'Suicide Vest with PETN in Urban Area',
    15, // meters (assessment radius)
    { type: 'cylinder', dimensions: { radius: 10, height: 30 } }, // cm
    2, // kg
    true,
    {
        material: 'Aluminum',
        density: 2.7,
        thickness: 0.3,
        yieldStrength: 240,
        ultimateStrength: 310,
        elasticModulus: 70,
        hardness: 75,
        elongation: 12,
        thermalConductivity: 167,
        meltingPoint: 933,
        fractureToughness: 24,
        poissonRatio: 0.33,
        heatCapacity: 900,
        corrosionResistance: 7,
        fragmentationPattern: 'natural'
    },
    environments[0], // Standard atmosphere
    [targets[0], targets[2]] // Human and vehicle targets
);

runComprehensiveScenario(
    'IED with ANFO in Transit Hub',
    15, // meters
    { type: 'cube', dimensions: { side: 25 } }, // cm
    10, // kg
    false,
    {
        material: 'Steel (mild)',
        density: 7.85,
        thickness: 0.5,
        yieldStrength: 250,
        ultimateStrength: 420,
        elasticModulus: 200,
        hardness: 130,
        elongation: 25,
        thermalConductivity: 45,
        meltingPoint: 1723,
        fractureToughness: 50,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 4,
        fragmentationPattern: 'natural'
    },
    { ...environments[0], reflectiveSurfaces: true }, // Standard atmosphere with reflection
    [targets[0], targets[1]] // Human and concrete targets
);

runComprehensiveScenario(
    'VBIED with C4 Near Government Building',
    30, // meters
    { type: 'cylinder', dimensions: { radius: 15, height: 50 } }, // cm
    20, // kg
    false,
    {
        material: 'Steel (high strength)',
        density: 7.85,
        thickness: 1.0,
        yieldStrength: 550,
        ultimateStrength: 800,
        elasticModulus: 210,
        hardness: 280,
        elongation: 15,
        thermalConductivity: 40,
        meltingPoint: 1723,
        fractureToughness: 80,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 5,
        fragmentationPattern: 'controlled'
    },
    environments[0], // Standard atmosphere
    [targets[1], targets[0]] // Concrete and human targets
);

runComprehensiveScenario(
    'Thermobaric Bomb in Enclosed Market',
    15, // meters
    { type: 'sphere', dimensions: { radius: 15 } }, // cm
    5, // kg
    true,
    {
        material: 'Aluminum',
        density: 2.7,
        thickness: 0.4,
        yieldStrength: 240,
        ultimateStrength: 310,
        elasticModulus: 70,
        hardness: 75,
        elongation: 12,
        thermalConductivity: 167,
        meltingPoint: 933,
        fractureToughness: 24,
        poissonRatio: 0.33,
        heatCapacity: 900,
        corrosionResistance: 7,
        fragmentationPattern: 'natural'
    },
    { ...environments[0], reflectiveSurfaces: true }, // Enclosed environment
    [targets[0], targets[1]] // Human and concrete targets
);

runComprehensiveScenario(
    'Shaped Charge Vest with HMX in Military Checkpoint',
    10, // meters
    {
        type: 'shaped_charge',
        dimensions: { radius: 10, height: 20, standoff: 5 }, // cm
        linerProperties: { material: 'copper', density: 8.96, thickness: 3, angle: 42, standoff: 2 }
    },
    3, // kg
    true,
    {
        material: 'Titanium Alloy',
        density: 4.5,
        thickness: 0.5,
        yieldStrength: 880,
        ultimateStrength: 950,
        elasticModulus: 110,
        hardness: 320,
        elongation: 14,
        thermalConductivity: 7.2,
        meltingPoint: 1941,
        fractureToughness: 100,
        poissonRatio: 0.34,
        heatCapacity: 530,
        corrosionResistance: 9,
        fragmentationPattern: 'controlled'
    },
    environments[1], // Desert environment
    [targets[2]] // Vehicle target
);

runComprehensiveScenario(
    'Pipe Bomb with Dynamite in Shopping Mall',
    10, // meters
    { type: 'cylinder', dimensions: { radius: 5, height: 20 } }, // cm
    1, // kg
    true,
    {
        material: 'Steel (mild)',
        density: 7.85,
        thickness: 0.3,
        yieldStrength: 250,
        ultimateStrength: 420,
        elasticModulus: 200,
        hardness: 130,
        elongation: 25,
        thermalConductivity: 45,
        meltingPoint: 1723,
        fractureToughness: 50,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 4,
        fragmentationPattern: 'natural'
    },
    { ...environments[0], reflectiveSurfaces: true }, // Indoor with reflection
    [
        targets[0], // Human
        { material: 'glass', thickness: 0.5, density: 2.5, distance: 10, crossSection: 2, criticalDamageThreshold: 15 }
    ]
);

runComprehensiveScenario(
    'Backpack Bomb with Semtex at Outdoor Event',
    20, // meters
    { type: 'sphere', dimensions: { radius: 12 } }, // cm
    3, // kg
    true,
    {
        material: 'Aluminum',
        density: 2.7,
        thickness: 0.4,
        yieldStrength: 240,
        ultimateStrength: 310,
        elasticModulus: 70,
        hardness: 75,
        elongation: 12,
        thermalConductivity: 167,
        meltingPoint: 933,
        fractureToughness: 24,
        poissonRatio: 0.33,
        heatCapacity: 900,
        corrosionResistance: 7,
        fragmentationPattern: 'natural'
    },
    environments[0], // Standard atmosphere, open area
    [
        targets[0], // Human at 10 m
        { ...targets[0], distance: 20 } // Human at 20 m
    ]
);

runComprehensiveScenario(
    'Truck Bomb with Comp B Near Industrial Facility',
    50, // meters
    { type: 'cube', dimensions: { side: 50 } }, // cm
    100, // kg
    false,
    {
        material: 'Steel (high strength)',
        density: 7.85,
        thickness: 1.5,
        yieldStrength: 550,
        ultimateStrength: 800,
        elasticModulus: 210,
        hardness: 280,
        elongation: 15,
        thermalConductivity: 40,
        meltingPoint: 1723,
        fractureToughness: 80,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 5,
        fragmentationPattern: 'natural'
    },
    { ...environments[1], soilType: 'sandy', soilDensity: 1.6 }, // Desert with soil
    [targets[1], targets[0]] // Concrete and human targets
);

runComprehensiveScenario(
    'Underwater IED with Nitroglycerin at Port',
    10, // meters
    { type: 'sphere', dimensions: { radius: 10 } }, // cm
    2, // kg
    true,
    {
        material: 'Titanium Alloy',
        density: 4.5,
        thickness: 0.5,
        yieldStrength: 880,
        ultimateStrength: 950,
        elasticModulus: 110,
        hardness: 320,
        elongation: 14,
        thermalConductivity: 7.2,
        meltingPoint: 1941,
        fractureToughness: 100,
        poissonRatio: 0.34,
        heatCapacity: 530,
        corrosionResistance: 9,
        fragmentationPattern: 'controlled'
    },
    environments[3], // Underwater environment
    [
        { material: 'steel', thickness: 2, density: 7.85, distance: 5, crossSection: 10, criticalDamageThreshold: 100 }
    ]
);

runComprehensiveScenario(
    'Drone-Delivered Explosive with RDX in Arctic Base',
    15, // meters
    { type: 'cylinder', dimensions: { radius: 8, height: 25 } }, // cm
    1.5, // kg
    true,
    {
        material: 'Aluminum',
        density: 2.7,
        thickness: 0.3,
        yieldStrength: 240,
        ultimateStrength: 310,
        elasticModulus: 70,
        hardness: 75,
        elongation: 12,
        thermalConductivity: 167,
        meltingPoint: 933,
        fractureToughness: 24,
        poissonRatio: 0.33,
        heatCapacity: 900,
        corrosionResistance: 7,
        fragmentationPattern: 'controlled'
    },
    { ...environments[2], soilType: 'frozen', soilDensity: 1.9 }, // Arctic environment
    [targets[0], targets[1]] // Human and concrete targets
);

runComprehensiveScenario(
    'Holographic Disguise with Embedded Micro-Explosives',
    10, // meters (assessment radius)
    { type: 'cylinder', dimensions: { radius: 0.5, height: 1 } }, // cm (individual micro-charge)
    2, // kg total
    true,
    {
        material: 'Non-metallic',
        density: 1.0, // g/cm³ (approximate for flexible material)
        thickness: 0.1, // cm
        yieldStrength: 10, // MPa (low strength)
        ultimateStrength: 20, // MPa
        elasticModulus: 1, // GPa
        hardness: 10, // Brinell
        elongation: 50, // %
        thermalConductivity: 0.1, // W/(m·K)
        meltingPoint: 373, // K
        fractureToughness: 1, // MPa·m^(1/2)
        poissonRatio: 0.4,
        heatCapacity: 1000, // J/(kg·K)
        corrosionResistance: 5,
        fragmentationPattern: 'natural'
    },
    { ...environments[0], reflectiveSurfaces: true }, // Indoor with reflection (Standard Atmosphere modified)
    [
        { material: 'human', thickness: 30, density: 1.0, distance: 5, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'reinforced concrete', thickness: 30, density: 2.4, distance: 10, crossSection: 10, criticalDamageThreshold: 70 }
    ]
);

runComprehensiveScenario(
    'Swarm of Explosive-Laden Micro-Drones',
    3, // meters (per drone)
    { type: 'sphere', dimensions: { radius: 2.5 } }, // cm (approximate for 50 g charge)
    0.05, // kg per drone
    true,
    {
        material: 'Plastic',
        density: 1.2, // g/cm³
        thickness: 0.1, // cm
        yieldStrength: 30, // MPa
        ultimateStrength: 50, // MPa
        elasticModulus: 2, // GPa
        hardness: 20, // Brinell
        elongation: 20, // %
        thermalConductivity: 0.2, // W/(m·K)
        meltingPoint: 373, // K
        fractureToughness: 2, // MPa·m^(1/2)
        poissonRatio: 0.4,
        heatCapacity: 1000, // J/(kg·K)
        corrosionResistance: 6,
        fragmentationPattern: 'natural'
    },
    environments[0], // Standard atmosphere, open area
    [
        { material: 'human', thickness: 30, density: 1.0, distance: 3, crossSection: 0.75, criticalDamageThreshold: 35 }
    ]
);

runComprehensiveScenario(
    'Explosive-Laced Paint in Art Installation',
    10, // meters
    { type: 'cylinder', dimensions: { radius: 15, height: 170 } }, // cm (sculpture dimensions)
    3, // kg
    false,
    {
        material: 'Plaster',
        density: 1.5, // g/cm³
        thickness: 2, // cm
        yieldStrength: 5, // MPa
        ultimateStrength: 10, // MPa
        elasticModulus: 5, // GPa
        hardness: 15, // Brinell
        elongation: 1, // %
        thermalConductivity: 0.5, // W/(m·K)
        meltingPoint: 373, // K (decomposition)
        fractureToughness: 0.5, // MPa·m^(1/2)
        poissonRatio: 0.2,
        heatCapacity: 800, // J/(kg·K)
        corrosionResistance: 3,
        fragmentationPattern: 'natural'
    },
    { ...environments[0], reflectiveSurfaces: true }, // Indoor with reflection (Standard Atmosphere modified)
    [
        { material: 'human', thickness: 30, density: 1.0, distance: 5, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'artifact', thickness: 1, density: 1.5, distance: 10, crossSection: 1, criticalDamageThreshold: 50 }
    ]
);

runComprehensiveScenario(
    'Subterranean Acoustic Explosive in Subway System',
    20, // meters
    { type: 'cylinder', dimensions: { radius: 15, height: 170 } }, // cm
    10, // kg
    false,
    {
        material: 'Steel (mild)',
        density: 7.85, // g/cm³
        thickness: 1.0, // cm
        yieldStrength: 250, // MPa
        ultimateStrength: 420, // MPa
        elasticModulus: 200, // GPa
        hardness: 130, // Brinell
        elongation: 25, // %
        thermalConductivity: 45, // W/(m·K)
        meltingPoint: 1723, // K
        fractureToughness: 50, // MPa·m^(1/2)
        poissonRatio: 0.3,
        heatCapacity: 465, // J/(kg·K)
        corrosionResistance: 4,
        fragmentationPattern: 'natural'
    },
    { ...environments[0], reflectiveSurfaces: true }, // Tunnel environment (Standard Atmosphere modified)
    [
        { material: 'human', thickness: 30, density: 1.0, distance: 10, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'reinforced concrete', thickness: 30, density: 2.4, distance: 20, crossSection: 10, criticalDamageThreshold: 70 }
    ]
);

runComprehensiveScenario(
    'Bioluminescent Explosive in Marine Environment',
    15, // meters
    { type: 'sphere', dimensions: { radius: 15 } }, // cm
    4, // kg
    true,
    {
        material: 'Gelatinous',
        density: 1.0, // g/cm³
        thickness: 0.5, // cm
        yieldStrength: 0.1, // MPa
        ultimateStrength: 0.2, // MPa
        elasticModulus: 0.01, // GPa
        hardness: 1, // Brinell
        elongation: 100, // %
        thermalConductivity: 0.5, // W/(m·K)
        meltingPoint: 300, // K
        fractureToughness: 0.1, // MPa·m^(1/2)
        poissonRatio: 0.45,
        heatCapacity: 1500, // J/(kg·K)
        corrosionResistance: 2,
        fragmentationPattern: 'natural'
    },
    environments[3], // Underwater environment
    [
        { material: 'human', thickness: 30, density: 1.0, distance: 5, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'steel', thickness: 2, density: 7.85, distance: 15, crossSection: 10, criticalDamageThreshold: 50 }
    ]
);

runComprehensiveScenario(
    'IED with ANFO in Transit Hub',
    1, // meters
    { type: 'cube', dimensions: { side: 100 } }, // cm
    1000, // kg
    false,
    {
        material: 'Steel (mild)',
        density: 7.85,
        thickness: 0.5,
        yieldStrength: 250,
        ultimateStrength: 420,
        elasticModulus: 200,
        hardness: 130,
        elongation: 25,
        thermalConductivity: 45,
        meltingPoint: 1723,
        fractureToughness: 50,
        poissonRatio: 0.3,
        heatCapacity: 465,
        corrosionResistance: 4,
        fragmentationPattern: 'natural'
    },
    { ...environments[0], reflectiveSurfaces: true }, // Standard atmosphere with reflection
    [
        { material: 'human', thickness: 30, density: 1.0, distance: 5, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'human', thickness: 30, density: 1.0, distance: 10, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'human', thickness: 30, density: 1.0, distance: 15, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'human', thickness: 30, density: 1.0, distance: 20, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'human', thickness: 30, density: 1.0, distance: 25, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'human', thickness: 30, density: 1.0, distance: 30, crossSection: 0.75, criticalDamageThreshold: 35 },
        { material: 'reinforced concrete', thickness: 30, density: 2.4, distance: 10, crossSection: 10, criticalDamageThreshold: 70 }
    ] // Human and concrete targets
);