Browse Source

Testing linear analisys

lintest
Michael Uleysky 4 months ago
parent
commit
aaafdf28e7
  1. 2
      actions/actionuv.h
  2. 52
      include/uvdata.h

2
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("div", "", "Velocity divergence", "(" + velunit + ")/" + distunit, "");
if(!err.Exist() && !headwrited) err = fw.AddVariable("rot", "", "Velocity rotor", "(" + 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("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("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.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); 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, "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, "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, "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(data, "ke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j); }, it);
if(!err.Exist()) if(!err.Exist())
err = fw.WriteVariable( err = fw.WriteVariable(

52
include/uvdata.h

@ -293,6 +293,58 @@ template<class OneVarData> class UVData: public UVDataDims<OneVarData, internal:
return sn * sn + ss * ss - w * w; 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 auto StablePoints(size_t ix, size_t iy) const
{ {
std::vector<struct B::StPoint> points; std::vector<struct B::StPoint> points;

Loading…
Cancel
Save