using System;
using UnityEngine;

// EarthCoreCalculations.cs
// This static class provides specific calculations related to the Earth's core dynamics,
// implementing the terms of the Navier-Stokes, Magnetohydrodynamic (MHD), Energy, and Compositional equations.
// These methods are designed to be called for individual grid points within a discretized core.

public static class EarthCoreCalculations
{
    // --- Existing Moment of Inertia and Density Calculation Methods ---
    // (Assuming these are already present and correct, as per your previous files)

    // public static double CalculateSolidSphereMomentOfInertia(double mass_kg, double radius_m) { ... }
    // public static double CalculateSphericalShellMomentOfInertia(double mass_kg, double innerRadius_m, double outerRadius_m) { ... }
    // public static double CalculateAverageDensity(double mass_kg, double volume_m3) { ... }


    // --- Core Dynamics Equation Terms ---
    // These methods represent the components of the Navier-Stokes and MHD equations.
    // They will be called by CoreDynamicsSimulator for each grid point.
    // The coefficients and relationships within these methods will be empirically determined by you.

    /// <summary>
    /// Calculates the local density of the outer core fluid, accounting for
    /// thermal and compositional variations. This is a simplified equation of state.
    /// </summary>
    /// <param name="referenceDensity">Reference density (EarthConstantsAndCalculations.OUTER_CORE_REFERENCE_DENSITY_KG_M3).</param>
    /// <param name="referenceTemperature">A reference temperature for the core (K).</param>
    /// <param name="currentTemperature">Current local temperature (K).</param>
    /// <param name="thermalExpCoeff">Thermal expansion coefficient (EarthConstantsAndCalculations.OUTER_CORE_THERMAL_EXPANSION_COEFFICIENT_PER_K).</param>
    /// <param name="compositionalAnomaly">Local compositional anomaly (e.g., concentration of light elements).</param>
    /// <param name="compositionalExpCoeff">Compositional expansion coefficient (EarthConstantsAndCalculations.OUTER_CORE_COMPOSITIONAL_EXPANSION_COEFFICIENT).</param>
    /// <returns>Local density in kg/m^3.</returns>
    public static double CalculateLocalOuterCoreDensity(
        double referenceDensity,
        double referenceTemperature, // e.g., average core temperature
        double currentTemperature,
        double thermalExpCoeff,
        double compositionalAnomaly, // e.g., deviation from average light element concentration
        double compositionalExpCoeff)
    {
        // Simplified linear equation of state for density variations: rho = rho_0 * [1 - alpha * (T - T_0) + beta * (C - C_0)]
        // Where alpha is thermalExpCoeff, beta is compositionalExpCoeff
        // NEED EMPIRICAL DATA for the exact form and coefficients for Earth's core.
        double density = referenceDensity * (1.0 - thermalExpCoeff * (currentTemperature - referenceTemperature) + compositionalExpCoeff * compositionalAnomaly);
        return density;
    }

    /// <summary>
    /// Calculates the Coriolis force term (2 * rho * Omega x v) for a fluid element.
    /// Part of the Navier-Stokes equation.
    /// </summary>
    /// <param name="localDensity_kg_m3">Local density of the fluid element.</param>
    /// <param name="earthAngularVelocity_rad_s">Earth's angular velocity vector (Omega).</param>
    /// <param name="fluidVelocity_m_s">Local fluid velocity vector (v).</param>
    /// <returns>Coriolis force vector (N/m^3, force per unit volume).</returns>
    public static Vector3d CalculateCoriolisForceTerm(
        double localDensity_kg_m3,
        Vector3d earthAngularVelocity_rad_s,
        Vector3d fluidVelocity_m_s)
    {
        // F_coriolis = 2 * rho * (Omega x v)
        return 2.0 * localDensity_kg_m3 * earthAngularVelocity_rad_s.Cross(fluidVelocity_m_s);
    }

    /// <summary>
    /// Calculates the Buoyancy force term ((rho_0 - rho) * g) for a fluid element.
    /// Part of the Navier-Stokes equation.
    /// </summary>
    /// <param name="localDensity_kg_m3">Local density of the fluid element.</param>
    /// <param name="referenceDensity_kg_m3">Reference density for buoyancy (EarthConstantsAndCalculations.OUTER_CORE_REFERENCE_DENSITY_KG_M3).</param>
    /// <param name="gravityVector_m_s2">Local gravitational acceleration vector (g).</param>
    /// <returns>Buoyancy force vector (N/m^3, force per unit volume).</returns>
    public static Vector3d CalculateBuoyancyForceTerm(
        double localDensity_kg_m3,
        double referenceDensity_kg_m3,
        Vector3d gravityVector_m_s2)
    {
        // Buoyancy force = (rho_0 - rho) * g
        // If localDensity_kg_m3 is less than referenceDensity_kg_m3, this force will be opposite to gravity (upwards).
        // NEED EMPIRICAL DATA for the reference density and precise buoyancy formulation.
        return (referenceDensity_kg_m3 - localDensity_kg_m3) * gravityVector_m_s2;
    }

    /// <summary>
    /// Calculates the Lorentz force term ((1/mu) * (curl B) x B) for a fluid element.
    /// Part of the Navier-Stokes equation and the induction equation.
    /// </summary>
    /// <param name="currentMagneticField_T">Local magnetic field vector (B).</param>
    /// <param name="magneticPermeability_H_m">Magnetic permeability (mu).</param>
    /// <param name="curlB">The curl of the magnetic field (nabla x B), calculated from neighbors.</param>
    /// <returns>Lorentz force vector (N/m^3, force per unit volume).</returns>
    public static Vector3d CalculateLorentzForceTerm(
        Vector3d currentMagneticField_T,
        double magneticPermeability_H_m,
        Vector3d curlB) // This parameter is passed from the simulator
    {
        // F_Lorentz = (1/mu) * (curl B) x B
        if (magneticPermeability_H_m == 0) return Vector3d.zero; // Avoid division by zero
        return (1.0 / magneticPermeability_H_m) * curlB.Cross(currentMagneticField_T);
    }

    /// <summary>
    /// Calculates the Viscous force term (nu * nabla^2 v) for a fluid element.
    /// Part of the Navier-Stokes equation.
    /// </summary>
    /// <param name="dynamicViscosity_Pa_s">Dynamic viscosity (nu).</param>
    /// <param name="laplacianV">The Laplacian of the velocity vector (nabla^2 v), calculated from neighbors.</param>
    /// <returns>Viscous force vector (N/m^3, force per unit volume).</returns>
    public static Vector3d CalculateViscousForceTerm(
        double dynamicViscosity_Pa_s,
        Vector3d laplacianV) // This parameter is passed from the simulator
    {
        // F_viscous = nu * nabla^2 v
        return dynamicViscosity_Pa_s * laplacianV;
    }

    /// <summary>
    /// Calculates the pressure gradient term (-nabla P) from the sum of other forces.
    /// In a full Navier-Stokes solver, pressure is often solved implicitly or through a Poisson equation.
    /// For an empirical model, this might represent the net force that *would* be balanced by a pressure gradient.
    /// </summary>
    /// <param name="inertialTerm">rho * (dv/dt + (v.grad)v) - Placeholder for future dv/dt and advection terms.</param>
    /// <param name="coriolisTerm">2 * rho * Omega x v</param>
    /// <param name="buoyancyTerm">(rho_0 - rho) * g</param>
    /// <param name="lorentzTerm">(1/mu) * (curl B) x B</param>
    /// <param name="viscousTerm">nu * nabla^2 v</param>
    /// <returns>The conceptual pressure gradient vector (Pa/m).</returns>
    public static Vector3d CalculatePressureGradientFromForces(
        Vector3d inertialTerm, // Placeholder for future dv/dt and advection terms
        Vector3d coriolisTerm,
        Vector3d buoyancyTerm,
        Vector3d lorentzTerm,
        Vector3d viscousTerm)
    {
        // -∇P = inertialTerm + coriolisTerm + buoyancyTerm + lorentzTerm + viscousTerm
        // This is the sum of all forces (per unit volume) that the pressure gradient must balance.
        // NEED EMPIRICAL DATA for the relative importance and exact balance of these terms.
        return inertialTerm + coriolisTerm + buoyancyTerm + lorentzTerm + viscousTerm;
    }


    /// <summary>
    /// Calculates the advection term (curl(v x B)) for the magnetic induction equation.
    /// This describes how the magnetic field is "advected" by fluid flow.
    /// </summary>
    /// <param name="vCrossB_curl">The curl of (v x B), calculated from neighbors.</param>
    /// <returns>The advection term for magnetic field change (T/s).</returns>
    public static Vector3d CalculateMagneticAdvectionTerm(Vector3d vCrossB_curl)
    {
        return vCrossB_curl; // dB/dt_advection = curl(v x B)
    }

    /// <summary>
    /// Calculates the diffusion term ((1/(mu*sigma)) * nabla^2 B) for the magnetic induction equation.
    /// This describes how the magnetic field diffuses through the fluid.
    /// </summary>
    /// <param name="magneticDiffusivity">Magnetic diffusivity (eta = 1 / (mu * sigma)).</param>
    /// <param name="laplacianB">The Laplacian of the magnetic field (nabla^2 B), calculated from neighbors.</param>
    /// <returns>The diffusion term for magnetic field change (T/s).</returns>
    public static Vector3d CalculateMagneticDiffusionTerm(
        double magneticDiffusivity,
        Vector3d laplacianB) // This parameter is passed from the simulator
    {
        // dB/dt_diffusion = eta * nabla^2 B
        return magneticDiffusivity * laplacianB;
    }

    // --- NEW: Thermal Dynamics Terms (Energy Equation) ---

    /// <summary>
    /// Calculates the thermal advection term (-v . nabla T) for the energy equation.
    /// Describes heat transport by fluid flow.
    /// </summary>
    /// <param name="fluidVelocity_m_s">Local fluid velocity vector (v).</param>
    /// <param name="gradTemperature">Gradient of temperature (nabla T), calculated from neighbors.</param>
    /// <returns>Change in temperature due to advection (K/s).</returns>
    public static double CalculateThermalAdvectionTerm(Vector3d fluidVelocity_m_s, Vector3d gradTemperature)
    {
        // -v . nabla T
        return -fluidVelocity_m_s.Dot(gradTemperature);
    }

    /// <summary>
    /// Calculates the thermal conduction term (kappa * nabla^2 T) for the energy equation.
    /// Describes heat diffusion through the material.
    /// </summary>
    /// <param name="thermalDiffusivity_m2_per_s">Thermal diffusivity (kappa).</param>
    /// <param name="laplacianTemperature">Laplacian of temperature (nabla^2 T), calculated from neighbors.</param>
    /// <returns>Change in temperature due to conduction (K/s).</returns>
    public static double CalculateThermalConductionTerm(double thermalDiffusivity_m2_per_s, double laplacianTemperature)
    {
        // kappa * nabla^2 T
        return thermalDiffusivity_m2_per_s * laplacianTemperature;
    }

    /// <summary>
    /// Calculates the heat source term (Q / (rho * Cp)) for the energy equation.
    /// Accounts for internal heat generation (e.g., radiogenic decay, latent heat).
    /// </summary>
    /// <param name="heatProduction_W_per_m3">Heat generated per unit volume (Watts/m^3).</param>
    /// <param name="localDensity_kg_m3">Local density (rho).</param>
    /// <param name="specificHeatCapacity_J_per_kg_K">Specific heat capacity (Cp).</param>
    /// <returns>Change in temperature due to heat sources (K/s).</returns>
    public static double CalculateHeatSourceTerm(double heatProduction_W_per_m3, double localDensity_kg_m3, double specificHeatCapacity_J_per_kg_K)
    {
        if (localDensity_kg_m3 == 0 || specificHeatCapacity_J_per_kg_K == 0) return 0;
        // Q / (rho * Cp)
        return heatProduction_W_per_m3 / (localDensity_kg_m3 * specificHeatCapacity_J_per_kg_K);
    }

    // --- NEW: Compositional Dynamics Terms ---

    /// <summary>
    /// Calculates the compositional advection term (-v . nabla C) for the compositional equation.
    /// Describes transport of compositional anomalies by fluid flow.
    /// </summary>
    /// <param name="fluidVelocity_m_s">Local fluid velocity vector (v).</param>
    /// <param name="gradCompositionalAnomaly">Gradient of compositional anomaly (nabla C), calculated from neighbors.</param>
    /// <returns>Change in compositional anomaly due to advection (dimensionless/s).</returns>
    public static double CalculateCompositionalAdvectionTerm(Vector3d fluidVelocity_m_s, Vector3d gradCompositionalAnomaly)
    {
        // -v . nabla C
        return -fluidVelocity_m_s.Dot(gradCompositionalAnomaly);
    }

    /// <summary>
    /// Calculates the compositional diffusion term (D * nabla^2 C) for the compositional equation.
    /// Describes diffusion of compositional anomalies.
    /// </summary>
    /// <param name="compositionalDiffusivity_m2_per_s">Compositional diffusivity (D).</param>
    /// <param name="laplacianCompositionalAnomaly">Laplacian of compositional anomaly (nabla^2 C), calculated from neighbors.</param>
    /// <returns>Change in compositional anomaly due to diffusion (dimensionless/s).</returns>
    public static double CalculateCompositionalDiffusionTerm(double compositionalDiffusivity_m2_per_s, double laplacianCompositionalAnomaly)
    {
        // D * nabla^2 C
        return compositionalDiffusivity_m2_per_s * laplacianCompositionalAnomaly;
    }
}

Planetary Instability Model PIM - Copyright (C) 2025 James Pacha