From aaafdf28e73cb3c646bc2036b5cf4dbcba50181e Mon Sep 17 00:00:00 2001
From: Michael Uleysky <uleysky@poi.dvo.ru>
Date: Wed, 3 Jul 2024 14:33:03 +1000
Subject: [PATCH] Testing linear analisys

---
 actions/actionuv.h |  2 ++
 include/uvdata.h   | 52 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 54 insertions(+)

diff --git a/actions/actionuv.h b/actions/actionuv.h
index ebe4d41..8f27322 100644
--- a/actions/actionuv.h
+++ b/actions/actionuv.h
@@ -237,6 +237,7 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
     if(!err.Exist() && !headwrited) err = fw.AddVariable("div", "", "Velocity divergence", "(" + velunit + ")/" + distunit, "");
     if(!err.Exist() && !headwrited) err = fw.AddVariable("rot", "", "Velocity rotor", "(" + velunit + ")/" + distunit, "");
     if(!err.Exist() && !headwrited) err = fw.AddVariable("ow", "", "Okubo-Weiss parameter", "(" + velunit + ")2/" + distunit + "2", "");
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("discr", "", "", "(" + velunit + ")2/" + distunit + "2", "");
     if(!err.Exist() && !headwrited) err = fw.AddVariable("ke", "", "Squared velocity module, u^2+v^2", "(" + velunit + ")2", "");
     if(!err.Exist() && !headwrited) err = fw.AddVariable("eke", "", "Squared velocity dispersion aka eddy kinetic energy, <u^2+v^2>-<u>^2-<v>^2", "(" + velunit + ")2", "");
     if(!err.Exist() && !headwrited) err = fw.WriteGrid(data);
@@ -246,6 +247,7 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
     if(!err.Exist()) err = fw.WriteVariable(data, "div", [&data = std::as_const(data)](size_t i, size_t j) { return data.Div(i, j); }, it);
     if(!err.Exist()) err = fw.WriteVariable(data, "rot", [&data = std::as_const(data)](size_t i, size_t j) { return data.Rot(i, j); }, it);
     if(!err.Exist()) err = fw.WriteVariable(data, "ow", [&data = std::as_const(data)](size_t i, size_t j) { return data.OW(i, j); }, it);
+    if(!err.Exist()) err = fw.WriteVariable(data, "discr", [&data = std::as_const(data)](size_t i, size_t j) { return data.Discr(i, j); }, it);
     if(!err.Exist()) err = fw.WriteVariable(data, "ke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j); }, it);
     if(!err.Exist())
      err = fw.WriteVariable(
diff --git a/include/uvdata.h b/include/uvdata.h
index d98e2c8..a3f0fe8 100644
--- a/include/uvdata.h
+++ b/include/uvdata.h
@@ -293,6 +293,58 @@ template<class OneVarData> class UVData: public UVDataDims<OneVarData, internal:
   return sn * sn + ss * ss - w * w;
  }
 
+ // Discriminant of characteristic equation
+ real Discr(size_t i) const { return Discr(i % nx, i / nx); }
+ real Discr(size_t ix, size_t iy) const
+ {
+  if(!*this) return 0.0;
+  real ux = dUdX(ix, iy);
+  real uy = dUdY(ix, iy);
+  real vx = dVdX(ix, iy);
+  real vy = dVdY(ix, iy);
+
+  if(ux == Fillval() || uy == Fillval() || vx == Fillval() || vy == Fillval()) return Fillval();
+
+  real sn = ux - vy;
+  real ss = vx + uy;
+  real w  = michlib::Abs(uy - vx);
+
+  if(sn * sn + 4.0 * uy * vx >= 0)
+   return michlib::Acos(w / Sqrt(ss * ss + sn * sn)) * 2.0 / M_PI;
+  else
+  {
+   real L = (ux + vy) * 0.5;
+   real W = michlib::Sqrt(-(sn * sn + 4.0 * uy * vx)) * 0.5;
+   michlib::message(W);
+
+   real C2 = uy / W;
+   real C4 = (vy - L) / W;
+
+   real A = 2.0 * C4;
+   real B = C4 * C4 + C2 * C2 - 1.0;
+
+   real t1 = michlib::Atan2(-B, A) * 0.5;
+   real t2 = michlib::Atan2(B, -A) * 0.5;
+
+   real d1, d2;
+   {
+    real tsn  = michlib::Sin(t1);
+    real tcn  = michlib::Cos(t1);
+    real temp = tcn + C4 * tsn;
+    d1        = C2 * C2 * tsn * tsn + temp * temp;
+   }
+   {
+    real tsn  = michlib::Sin(t2);
+    real tcn  = michlib::Cos(t2);
+    real temp = tcn + C4 * tsn;
+    d2        = C2 * C2 * tsn * tsn + temp * temp;
+   }
+
+   real mind = std::min(d1, d2), maxd = std::max(d1, d2);
+   return -michlib::Sqrt(mind / maxd);
+  }
+ }
+
  auto StablePoints(size_t ix, size_t iy) const
  {
   std::vector<struct B::StPoint> points;