using System.Numerics;
using Skeleton.Algebra.Extensions;
using Skeleton.Algebra.ScalarFieldStructure;
using Skeleton.DataStructure;
using Skeleton.DataStructure.Link;
using Skeleton.DataStructure.Packs;
using Skeleton.Utils;
using Skeleton.Utils.Helpers;
using Skeleton.Utils.InverseSampling;
using VirtualChemistry.Chemistry.Bonds.Implements;
using VirtualChemistry.Chemistry.Compounds.Implements;
using VirtualChemistry.Chemistry.Mixtures.Resolver;
using VirtualChemistry.Constants;
using VirtualChemistry.DataStructure.Packs;

namespace VirtualChemistry.Chemistry.Mixtures.Implements;
/// <summary>
/// mixture of many well mixed compounds that can not be seperated by psy methods
/// </summary>
public class HomogeneousMixture
{
    /// <summary>
    /// Constructor by default
    /// </summary>
    public HomogeneousMixture()
    {
        LogStateMatrix = new CacheItem<LieSU3>(x =>
            Compounds.Select(compound => compound.LogStateMatrix.GetFrom(x) * compound.Concentration)
                .DefaultIfEmpty(LieSU3.Zero)
                .SpLieSum()
        );
        StateMatrix = new(x => LogStateMatrix.GetFrom(x).Exp());
        CompressionElasticity = new(x =>
            {
                DiagR3 spectrum = new (0d, 2d * Math.PI, 4d * Math.PI);
                return (LogStateMatrix.GetFrom(x).CayleyNegative() * StateMatrix.GetFrom(x))
                    .ConjugateOn(spectrum)
                    .Rayleigh(LogStateMatrix.GetFrom(x).ColumnAverage);
            }
        );
        CompressionBias = new(x =>
            {
                DiagR3 w1 = new (1d / 11d, 7d / 11d, 10d / 11d);
                DiagR3 w2 = new (2d / 17d, 8d / 17d, 16d / 17d);
                double p1 = StateMatrix.GetFrom(x)
                    .ConjugateOn(w1)
                    .Rayleigh(StateMatrix.GetFrom(x).Hermitian().ColumnAverage);
                double p2 = StateMatrix.GetFrom(x)
                    .Hermitian()
                    .ConjugateOn(w2)
                    .Rayleigh(StateMatrix.GetFrom(x).ColumnAverageBivariant);
                double x1 = InverseSampling.StandardNormal(p1) / 4d - 2d * Math.PI;
                double x2 = InverseSampling.StandardNormal(p2) / 3d + 2d * Math.PI;
                return Math.Atan(Math.Abs(x1) > Math.Abs(x2) ? x1 : x2) * 2d;
            }
        );
        FreeDensity = new(x =>
            {
                double avg = (Math.Exp(-4d * Math.PI) + Math.Exp(Math.PI)) / 2d;
                DiagR3 spectrum = new DiagR3(Math.Exp(-4d * Math.PI), avg, Math.Exp(Math.PI));
                double tempMod = Math.Exp(-Phase);
                return tempMod * LogStateMatrix.GetFrom(x)
                    .CayleyPositive()
                    .ConjugateOn(spectrum)
                    .Rayleigh(StateMatrix.GetFrom(x).ColumnAverageBivariant);
            }
        );
        HeatConductivity = new(x => 
            StateMatrix.GetFrom(x).ConjugateOn(ChemistryConstant.EnergyConductivitySpectrum)
            .Rayleigh(LogStateMatrix.GetFrom(x) * StateMatrix.GetFrom(x).ColumnAverageBivariant)
        );
        ResetAllCache();
    }

    /// <summary>
    /// singleton used to represent no mixture
    /// </summary>
    public static readonly HomogeneousMixture Null = new();

    private IEnumerable<BaseBond> ConnectedIndex()
    {
        HashSet<BaseBond> bs = new();
        foreach (Compound c in Compounds)
        {
            if(c.Amount.IsEqualApprox(0))
                continue;


            foreach (BaseBond b in c.Bonds)
            {
                if (
                    b.Connected &&
                    b.Energy > b.AntiBondingEnergy &&
                    b.ConnectedBond.Energy > b.ConnectedBond.AntiBondingEnergy &&
                    !bs.Contains(b.ConnectedBond)
                )
                {
                    bs.Add(b);
                }
            }
        }

        return bs;
        //Compounds.Where(c => !(0.5 * c.Amount).IsEqualApprox(0)).SelectMany(x => x.ConnectedIndex());
    }
    public void RPC()
    {

        int s1 = 0;
        int s2 = 0;
        int s3 = 0;
        int s4 = 0;
        foreach (Compound c in Compounds)
        {
            foreach (BaseBond b in c.Bonds)
            {
                switch (b.BondingGroup)
                {
                    case 0:
                        s1++;
                        break;
                    case 1:
                        s2++;
                        break;
                    case 2:
                        s3++;
                        break;
                    default:
                        s4++;
                        break;
                }
            }
        }
        Console.WriteLine($"{s1} ----  {s2} ---- {s3} --- {s4}");
    }

    private BondGroups ReactionGroup()
    {
        HashSet<BaseBond> group1 = new();
        HashSet<BaseBond> group2 = new();
        HashSet<BaseBond> group3 = new();
        HashSet<BaseBond> group4 = new();
        foreach (Compound c in Compounds)
        {
            if ((c.Amount * 0.125).IsEqualApprox(0) || c.Atoms.Count > 15) 
                continue;
            foreach(BaseBond b in c.Bonds)
            {
                if(b.Connected || !(b.Energy > b.BondingEnergy))
                    continue;
                switch (b.BondingGroup)
                {
                    case 0:
                        group1.Add(b);
                        break;
                    case 1:
                        group2.Add(b);
                        break;
                    case 2:
                        group3.Add(b);
                        break;
                    case 3:
                        group4.Add(b);
                        break;
                }
            }
        }
        return new(group1, group2, group3, group4);
    }

    /// <summary>
    /// elasticity range: 0, 4pi
    /// elasticity controls the speed k drops from e^pi to e^-4pi
    /// elasticity -> 0; k -> constant
    /// b -> 4pi; k -> e^pi
    /// b -> -4pi; k -> e^-4pi
    /// </summary>
    public CacheItem<double> CompressionElasticity { get; }

    /// <summary>
    /// bias range: -4pi, 4pi
    /// bias controls the temp of triple point
    /// higher b -> higher melting point and higher boiling point
    /// </summary>
    public CacheItem<double> CompressionBias { get; }

    /// <summary>
    /// phase -> 1 implies material is more like gas
    /// phase -> -1 implies material is more like solid
    /// 
    /// </summary>
    public double Phase => Compounds.Sum(c => c.Concentration * c.Phase);

    /// <summary>
    /// stiffness range: e^-4pi, e^pi
    /// higher stiffness -> harder to compress -> solid like
    /// lower stiffness -> easier to compress -> gas like
    /// Stiffness K, Density D, Free Density D0, Temperature T
    /// P = K(1/D - 1/D0)
    /// 
    /// Temperature range: -pi, pi
    /// </summary>
    public double CompressionStiffness => Math.Exp(Math.PI / 2d * (5d * Phase - 3d));

    /// <summary>
    /// e^-4pi ... e^pi
    /// </summary>
    public CacheItem<double> FreeDensity { get; }

    /// <summary>
    /// volume that under 0 compression
    /// </summary>
    public double FreeVolume => Amount / FreeDensity.Get;

    /// <summary>
    /// P = K (1/D - 1/D0)
    /// </summary>
    public double Pressure => (Volume - FreeVolume) * CompressionStiffness;

    /// <summary>
    /// Volume of the mixture
    /// </summary>
    public double Volume { get; set; }

    /// <summary>
    /// Heat transfer ability
    /// </summary>
    public CacheItem<double> HeatConductivity { get; }

    /// <summary>
    /// energy / amount = temperature
    /// </summary>
    public double Temperature
    {
        get => Energy / Amount;
        set => Energy = value * Amount;
    }

    /// <summary>
    /// combust rate
    /// </summary>
    public double CombustRate
    {
        get
        {
            H3 w = StateMatrix.Get!.ConjugateOn(new DiagR3(0, Math.PI, 2 * Math.PI));
            return w.Rayleigh(
                LogStateMatrix.Get!.ColumnAverageBivariant +
                StateMatrix.Get.Hermitian().ColumnAverageBivariant
            );
        }
    }

    /// <summary>
    /// total Energy the mixture holds 
    /// </summary>
    public double Energy
    {
        get => Compounds.Sum(c => c.Energy);
        set
        {
            foreach (Compound c in Compounds)
                c.Energy = value * c.Concentration;
        }
    }
    /// <summary>
    /// reset all cache
    /// </summary>
    public void ResetAllCache()
    {
        LayerCache = LinkNode<HomogeneousMixture>.Null;
    }
    /// <summary>
    /// the heterogeneous mixture that holds this homogeneous mixture
    /// </summary>
    public HeterogeneousMixture HeterogeneousMixture { get; set; } = HeterogeneousMixture.Null;

    /// <summary>
    /// Compounds of this mixture
    /// </summary>
    public HashSet<Compound> Compounds { get; set; } = new HashSet<Compound>();
    /// <summary>
    /// Compounds whose ratio is greater than 1%
    /// </summary>
    public IEnumerable<Compound> MainCompounds => Compounds.Where(c => c.Concentration > 0.001d);
    private LinkNode<HomogeneousMixture> LayerCache { get; set; } = LinkNode<HomogeneousMixture>.Null;
    /// <summary>
    /// 
    /// </summary>
    public LinkNode<HomogeneousMixture> Layer
    {
        get => LayerCache == LinkNode<HomogeneousMixture>.Null
            ? HeterogeneousMixture.LayerOrder.Find(this)
            : LayerCache;
        set => LayerCache = value;
    }

    /// <summary>
    /// State matrix of this mixture
    /// </summary>
    public CacheItem<SU3> StateMatrix { get; }

    /// <summary>
    /// Log of the state matrix
    /// </summary>
    public CacheItem<LieSU3> LogStateMatrix { get; }
    
    /// <summary>
    /// Split the homogeneous mixture into another heterogeneous mixture by ratio or amount;
    /// </summary>
    /// <param name="amount"> The amount or ratio</param>
    /// <param name="into">another heterogeneous mixture</param>
    /// <param name="byRatio">determine if the amount variable is ratio or amount</param>
    /// <returns></returns>
    public HomogeneousMixture Split(double amount, HeterogeneousMixture into, bool byRatio = false)
    {
        double ebs = Energy;
        double vbs = Volume;
        double ratio = byRatio ? amount : amount / Amount;
        HomogeneousMixture res = new HomogeneousMixture();
        foreach (Compound c in Compounds)
        {
            Compound split = c.Split(amount, byRatio);
            if(split != Compound.Null)
                res.AddCompound(split);
        }
        res.Energy = ebs*ratio;
        res.Volume = vbs*ratio;
        Energy = ebs * (1 - ratio);
        Volume = vbs * (1 - ratio);
        into.AddLayer(res);
        return res;
    }
    /// <summary>
    /// Remove a compound from the homogeneous mixture
    /// </summary>
    /// <param name="compound"></param>
    public void RemoveCompound(Compound compound)
    {
        double deltaV = compound.Volume;
        Compounds.Remove(compound);
        compound.HomogeneousMixture = Null;
        Volume -= deltaV;
    }

    /// <summary>
    /// Add compound into the homogeneous mixture
    /// </summary>
    /// <param name="compound">compound to add</param>
    /// <param name="resetCache"></param>
    /// <param name="skipCombine">add without combine to isomorphic compound</param>
    public void AddCompound(Compound compound, bool resetCache = true, bool skipCombine=false)
    {
        if (compound.HomogeneousMixture == this)
            return;
        bool added = false;
        if(!skipCombine)
        {
            foreach (Compound sCompound in Compounds.Where(x => x.Expression == compound.Expression))
            {
                if (sCompound.IsometricTo(compound))
                {
                    sCompound.Amount += compound.Amount;
                    compound.Amount = 0d;
                    sCompound.Energy += compound.Energy;
                    added = true;
                    break;
                }
            }
        }

        compound.HomogeneousMixture = added ? Null : this;
        if (!added)
            Compounds.Add(compound);
        if(resetCache)
            ResetAllCache();
    }
    /// <summary>
    /// Amount of this mixture 
    /// </summary>
    public double Amount
    {
        get => Compounds
            .Select(c => c.Amount)
            .DefaultIfEmpty(0)
            .Sum();
        set
        {
            Dictionary<Compound, double> res = new Dictionary<Compound, double>();
            foreach (Compound c in Compounds)
                res[c] = value * c.Concentration;
            foreach (Compound c in Compounds)
                c.Amount = res[c];
        }
    }

    /// <summary>
    /// ratio of this mixture in the heterogeneous mixture
    /// </summary>
    public double Ratio =>
        Amount / (HeterogeneousMixture == HeterogeneousMixture.Null ? Amount : HeterogeneousMixture.Amount);
    
    
    /// <summary>
    /// density of the mixture, amount of molecular per volume
    /// </summary>
    public double Density => Amount / Volume;
    
    /// <summary>
    /// save mixture into string
    /// </summary>
    /// <returns></returns>
    public string Dump() =>
        $"`ENV?ENG{Energy.ExactDoubleString()}?V{Volume.ExactDoubleString()}`%ENV`COMPS{String.Join(HomogeneousMixtureResolver.CompoundsSplit, Compounds.Select(c => c.Dump()))}`%COMPS";

    /// <summary>
    /// update volume
    /// </summary>
    public void VolumeUpdate()
    {
        double sVolume = Volume + HeterogeneousMixture.ForceFromContainer / CompressionStiffness;
        Volume = Math.Max(sVolume, 1E-6);
    }

    /// <summary>
    /// exchange energy to neighbor mixtures 
    /// </summary>
    public void HeatExchange()
    {
        if (Layer != Layer.Parent.First)
        {
            HomogeneousMixture mPrevious = Layer.Previous.Value;
            double balanceTemp = (Energy + mPrevious.Energy) / (Amount + mPrevious.Amount);
            double dTemp = balanceTemp - Temperature;
            double conductivity = Math.Sqrt(HeatConductivity.Get * mPrevious.HeatConductivity.Get);
            double transferTemp = dTemp * conductivity;
            Temperature += transferTemp;
            mPrevious.Temperature -= transferTemp;
        }
        if (Layer != Layer.Parent.Last)
        {
            HomogeneousMixture mNext = Layer.Next.Value;
            double balanceTemp = (Energy + mNext.Energy) / (Amount + mNext.Amount);
            double dTemp = balanceTemp - Temperature;
            double conductivity = Math.Sqrt(HeatConductivity.Get * mNext.HeatConductivity.Get);
            double transferTemp = dTemp * conductivity;
            Temperature += transferTemp;
            mNext.Temperature -= transferTemp;
        }

        double envBalanceTemp = HeterogeneousMixture.Container.EnvironmentTemperature;
        double envDTemp = envBalanceTemp - Temperature;
        double envConductivity = HeatConductivity.Get;
        Temperature += envDTemp * envConductivity;
    }
    private void ReactConnecting(HashSet<BaseBond> g1, HashSet<BaseBond> g2)
    {
        double maxCd = -1;
        BaseBond xb1 = BaseBond.Null;
        BaseBond xb2 = BaseBond.Null;
        foreach (BaseBond b1 in g1)
        {
            foreach (BaseBond b2 in g2)
            {
                H3 cd = new(Complex.ImaginaryOne / 2 * (b1.LogState.Get + b2.LogState.Get));
                double cdX = cd.Rayleigh((b1.LogState.Get + b2.LogState.Get).ColumnAverage);
                if (cdX > maxCd)
                {
                    maxCd = cdX;
                    xb1 = b1;
                    xb2 = b2;
                }
            }
        }
        if (xb1 == BaseBond.Null)
            return;
        if (xb1 == xb2)
        {
            PlainPack2<Compound, BaseBond> s = xb1.Compound.SplitWithBondMap(0.5, xb2, true);
            xb2 = s.Item2;
        }

        if (xb1.Compound.Amount.IsSignificantlyGreaterThen(xb2.Compound.Amount))
        {
            double amountDiff = xb1.Compound.Amount - xb2.Compound.Amount;
            Compound oth = xb1.Compound.Split(amountDiff);
            xb1.HomogeneousMixture.AddCompound(oth);
        }

        if (xb1.Compound.Amount.IsSignificantlyLessThen(xb2.Compound.Amount))
        {
            double amountDiff = xb2.Compound.Amount - xb1.Compound.Amount;
            Compound oth = xb2.Compound.Split(amountDiff);
            xb2.HomogeneousMixture.AddCompound(oth);
        }
        xb1.Connect(xb2);
    }

    private int ReactionType()
    {
        Complex d = StateMatrix.Get.Trace();
        return d.Phase switch
        {
            > 7 * Math.PI / 4 or < Math.PI / 4 => 0,
            > Math.PI / 4 and < 3 * Math.PI / 4 => 1,
            > 3 * Math.PI / 4 and < 5 * Math.PI / 4 => 2,
            _ => 3
        };
    }

    /// <summary>
    /// react within the homogeneous mixture
    /// </summary>
    public void React()
    {
        BaseBond[] connected = ConnectedIndex().ToArray();
        foreach (BaseBond bond in connected)
        {
            if (bond.Compound.Amount.IsEqualApprox(0))
                continue;
            if (
                bond.Energy > bond.AntiBondingEnergy &&
                bond.ConnectedBond.Energy > bond.ConnectedBond.AntiBondingEnergy &&
                bond.DisconnectCondition &&
                bond.ConnectedBond.DisconnectCondition
            )
            {
                bond.Disconnect();
            }
            
        }
        for (int i = 1; i <= 3; i++)
        {
            BondGroups bg = ReactionGroup();
            int type = ReactionType();
            if(type == 0)
            {
                ReactConnecting(bg.Group1, bg.Group3);
                ReactConnecting(bg.Group2, bg.Group4);
            }
            else if (type == 1)
            {
                ReactConnecting(bg.Group1, bg.Group1);
                ReactConnecting(bg.Group3, bg.Group3);
            }
            else if (type == 2)
            {
                ReactConnecting(bg.Group2, bg.Group2);
                ReactConnecting(bg.Group3, bg.Group3);
            }
            else
            {
                ReactConnecting(bg.Group1, bg.Group4);
                ReactConnecting(bg.Group2, bg.Group3);
            }
        }

    }
    /// <summary>
    /// higher level mixture can resolve lower level mixture
    /// </summary>
    public double ResolveLevel => Compounds.Sum(c => c.Concentration * c.ResolveLevel);
    /// <summary>
    /// eigen resolvability of other mixture
    /// </summary>
    /// <param name="oth"></param>
    /// <returns></returns>
    public double EigenResolvability(HomogeneousMixture oth)
    {
        LieSU3 s = oth.LogStateMatrix.Get ^ LogStateMatrix.Get;
        C3 a = new(1, 0.5, 1);
        C3 b = new(0.5, 1, 0.5);
        Complex w = a * s * b;
        ComplexFieldStructure.Structure.Fix(w);
        double res = w.Phase;
        if (res > Math.PI)
            res = 2 * Math.PI - res;
        return res;
    }
    
    /// <summary>
    /// eigen resolvability from -pi to pi
    /// </summary>
    /// <param name="c"></param>
    /// <returns></returns>
    public double EigenResolvability(Compound c)
    {
        LieSU3 s = c.LogStateMatrix.Get ^ LogStateMatrix.Get;
        C3 a = new(1, 0.5, 1);
        C3 b = new(0.5, 1, 0.5);
        Complex w = a * s * b;
        w = ComplexFieldStructure.Structure.Fix(w);
        double res = w.Phase;
        if (res > Math.PI)
            res = 2 * Math.PI - res;
        return res;
    }
    /// <summary>
    /// calculate resolvability from eigen resolvability
    /// </summary>
    /// <param name="eigenRes"></param>
    /// <returns></returns>
    public double Resolvability(double eigenRes)
    {
        double res = -Math.Tan((Math.Sin(eigenRes) * Math.PI - Math.PI) / 2);
        if (Math.Abs(res) < 1E-6)
            return res;
        return res;
    }
    
    /// <summary>
    /// ability to resolve other compound
    /// </summary>
    /// <param name="c"></param>
    /// <returns></returns>
    public double Resolvability(Compound c)
    {
        LieSU3 s = c.LogStateMatrix.Get ^ LogStateMatrix.Get;
        C3 a = new(1, 0.5, 1);
        C3 b = new(0.5, 1, 0.5);
        Complex w = a * s * b;
        w = ComplexFieldStructure.Structure.Fix(w);
        double res = -Math.Tan((Math.Sin(w.Phase) * Math.PI - Math.PI) / 2);
        if (Math.Abs(res) < 1E-6)
            return 0;
        return res;
    }
    /// <summary>
    ///  a tensor of order 3 abstracted from this mixture
    /// </summary>
    public C3 Ita => Compounds
        .Select(c => c.Concentration * c.StateMatrix.Get.ColumnAverageBivariant)
        .Aggregate((a, b) => a + b);
    /// <summary>
    /// Red component  
    /// </summary>
    public Byte ColorRed => (Byte)
    (
        StateMatrix.Get
            .ConjugateOn(ChemistryConstant.OpticalFilter.RedDefinite)
            .Rayleigh(Ita) * 255
    );
    /// <summary>
    /// Green component
    /// </summary>
    public Byte ColorGreen => (Byte)
    (
        StateMatrix.Get
            .ConjugateOn(ChemistryConstant.OpticalFilter.GreenDefinite)
            .Rayleigh(Ita) * 255
    );
    /// <summary>
    /// Blue component
    /// </summary>
    public Byte ColorBlue => (Byte)
    (
        StateMatrix.Get
            .ConjugateOn(ChemistryConstant.OpticalFilter.BlueDefinite)
            .Rayleigh(Ita) * 255
    );
    /// <summary>
    /// Transparent component
    /// </summary>
    public Byte ColorTransparent => (Byte)
    (
        StateMatrix.Get
            .ConjugateOn(ChemistryConstant.OpticalFilter.OpacityDefinite)
            .Rayleigh(Ita) * 255
    );

    public double OsmoticLevel => StateMatrix.Get
        .ConjugateOn(new DiagR3(Phase, Math.PI, -Phase))
        .Rayleigh(LogStateMatrix.Get.ColumnAverage);
    /// <summary>
    /// Should not be used by any methods except ResolveNext and ResolvePrev
    /// </summary>
    private void PureResolve(HomogeneousMixture hOther)
    {
        HashSet< (Compound, double)> toAbsorb = new ();
        if (hOther.OsmoticLevel > OsmoticLevel)
            return;
        foreach (Compound c in hOther.Compounds)
        {
            double eigenRes = EigenResolvability(c);
            double maxRes = Resolvability(eigenRes);
            double oRes = hOther.Resolvability(c);
            if(oRes < 0 || maxRes < oRes)
                continue;
            if (eigenRes.IsEqualApprox(0))
            {
                toAbsorb.Add((c, c.Amount));
                continue;
            }
            double maxResAmount = Amount;
            foreach (Compound ct in Compounds)
                if (ct.IsometricTo(c))
                    maxResAmount -= ct.Amount;
            maxResAmount *= maxRes;
            if(maxResAmount < 0)
                continue;
            double splitAmount = Math.Min(maxResAmount, c.Amount);
            if(splitAmount.IsEqualApprox(0))
                continue;
            toAbsorb.Add((c, splitAmount));
        }
        foreach ((Compound c, double a) in toAbsorb)
        {
            Compound sp = c.Split(a);
            double v = sp.Volume;
            hOther.Volume -= v;
            Volume += v;
            /*if (c.Amount.IsEqualApprox(0))
                hOther.RemoveCompound(c);*/
            AddCompound(sp);
        }
        if(hOther.Amount.IsEqualApprox(0))
            HeterogeneousMixture.RemoveLayer(hOther);
    }
    private void ResolveNext()
    {
        if (HeterogeneousMixture == HeterogeneousMixture.Null || Layer.Next.IsTail)
            return;
        PureResolve(Layer.Next.Value);
    }
    private void ResolvePrevious()
    {
        if (HeterogeneousMixture == HeterogeneousMixture.Null || Layer.Previous.IsHead)
            return;
        PureResolve(Layer.Previous.Value);
    }
    
    private bool Resolve()
    {
        ResolveNext();
        ResolvePrevious();
        return false;
    }

    private bool Precipitate()
    {
        HashSet<(Compound, double)> toDeposit = new();
        foreach (Compound c in Compounds)
        {
            double maxRes = Resolvability(c);
            if (maxRes < 0 || maxRes > c.Concentration)
                continue;
            double maxResAmount = maxRes * (Amount - c.Amount);
            if(maxResAmount.IsEqualApprox(0))
                continue;
            double aDiff = c.Amount - maxResAmount;
            if (aDiff.IsSignificantlyGreaterThen(0) )
                toDeposit.Add((c, aDiff));
        }

        foreach ((Compound c, double a) in toDeposit)
        {
            Compound oth = c.Split(a);
            if (oth.FreeDensity.Get > FreeDensity.Get)
            {
                if (Layer.Next.IsEnding)
                {
                    HomogeneousMixture m = HeterogeneousMixture.AddLayer();
                    m.Layer.MoveAfter(Layer);
                    m.AddCompound(oth);
                    m.Volume = m.FreeVolume;
                    return true;
                }
                Layer.Next.Value.AddCompound(oth);
            }
            else
            {
                if (Layer.Previous.IsEnding)
                {
                    HomogeneousMixture m = HeterogeneousMixture.AddLayer();
                    m.Layer.MoveBefore(Layer);
                    m.AddCompound(oth);
                    m.Volume = m.FreeVolume;
                    return true;
                }
                Layer.Previous.Value.AddCompound(oth);
            }
        }
        return false;
    }

    /// <summary>
    /// 
    /// </summary>
    /// <returns>if any new layer created or deleted</returns>
    public bool PrecipitateAndResolve()
    {
        bool a = Resolve();
        bool b = Precipitate();
        return a || b;
    }

    /// <summary>
    /// combine identical compounds in the mixture
    /// </summary>
    public void Degenerate()
    {
        HashSet<Compound> toRemove = new();
        foreach (IGrouping<string, Compound> g in Compounds.GroupBy(c => c.Expression))
        {
            HashSet<Compound> classified = new();
            foreach (Compound c in g)
            {
                classified.Add(c);
                foreach (Compound cx in g.Where(x => !classified.Contains(x)))
                {
                    if (!c.IsometricTo(cx))
                        continue;
                    classified.Add(cx);
                    c.Amount += cx.Amount;
                    cx.Amount = 0;
                    toRemove.Add(cx);
                }
            }
        }

        foreach (Compound c in toRemove)
        {
            c.HomogeneousMixture = Null;
            Compounds.Remove(c);
        }
    }

}