Fix PureResolve

src/Chemistry/Mixtures/Implements/HomogeneousMixture.cs

@@ -564,6 +564,54 @@ public class HomogeneousMixture
     /// 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>
@@ -632,10 +680,17 @@ public class HomogeneousMixture
         HashSet<Compound> toAbsorb = new ();
         foreach (Compound c in hOther.Compounds)
         {
-            double maxRes = Resolvability(c);
+            double eigenRes = EigenResolvability(c);
+            double maxRes = Resolvability(eigenRes);
             double oRes = hOther.Resolvability(c);
-            if(oRes < 0 || maxRes < 0 || maxRes < oRes)
+            if(oRes < 0 || maxRes < oRes)
                 continue;
+            if (eigenRes.IsEqualApprox(0))
+            {
+                toAbsorb.Add(c);
+                continue;
+            }
+
             double maxResAmount = Amount;
             foreach (Compound ct in Compounds)
                 if (ct.IsometricTo(c))
@@ -652,7 +707,12 @@ public class HomogeneousMixture
             toAbsorb.Add(cSplit);
         }
         foreach (Compound c in toAbsorb)
+        {
+            double v = c.Volume;
+            hOther.Volume -= v;
+            Volume += v;
             AddCompound(c);
+        }
         if(hOther.Amount.IsEqualApprox(0))
             HeterogeneousMixture.RemoveLayer(hOther);
     }