You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
226 lines
7.7 KiB
226 lines
7.7 KiB
1 year ago
|
#pragma once
|
||
|
#include "geohelpers.h"
|
||
|
#include <vector>
|
||
|
|
||
|
using michlib::M_PI;
|
||
|
using michlib::real;
|
||
|
using michlib::Sqrt;
|
||
|
|
||
|
enum class Metric
|
||
|
{
|
||
|
EUCLID,
|
||
|
SPHERIC
|
||
|
};
|
||
|
|
||
|
class UVData
|
||
|
{
|
||
|
static constexpr real fillval = 1.0e10;
|
||
|
real x0 = 0.0, y0 = 0.0;
|
||
|
size_t nx = 0, ny = 0;
|
||
|
real xstep = 0.0, ystep = 0.0;
|
||
|
|
||
|
std::vector<real> u, v;
|
||
|
Metric metric;
|
||
|
|
||
|
real D(real lon1, real lat1, real lon2, real lat2) const
|
||
|
{
|
||
|
switch(metric)
|
||
|
{
|
||
|
case(Metric::EUCLID): return michlib::Hypot(lon2 - lon1, lat2 - lat1);
|
||
|
// Spheric distance in km
|
||
|
case(Metric::SPHERIC): return michlib::GCD(lon1 * M_PI / 180.0, lat1 * M_PI / 180.0, lon2 * M_PI / 180.0, lat2 * M_PI / 180.0) * 6371.0;
|
||
|
}
|
||
|
return 0.0;
|
||
|
}
|
||
|
|
||
|
// For derivatives
|
||
|
real Ud(size_t ix, size_t iy) const { return IsCoast(ix, iy) ? 0.0 : U(ix, iy); }
|
||
|
real Vd(size_t ix, size_t iy) const { return IsCoast(ix, iy) ? 0.0 : V(ix, iy); }
|
||
|
|
||
|
real dUdX(size_t ix, size_t iy) const
|
||
|
{
|
||
|
if(IsCoast(ix, iy)) return fillval;
|
||
|
if(ix == 0) return (Ud(1, iy) - Ud(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy));
|
||
|
if(ix == nx - 1) return (Ud(nx - 1, iy) - Ud(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy));
|
||
|
return 0.5 * ((Ud(ix + 1, iy) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
|
||
|
(Ud(ix, iy) - Ud(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
|
||
|
}
|
||
|
real dVdX(size_t ix, size_t iy) const
|
||
|
{
|
||
|
if(IsCoast(ix, iy)) return fillval;
|
||
|
if(ix == 0) return (Vd(1, iy) - Vd(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy));
|
||
|
if(ix == nx - 1) return (Vd(nx - 1, iy) - Vd(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy));
|
||
|
return 0.5 * ((Vd(ix + 1, iy) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
|
||
|
(Vd(ix, iy) - Vd(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
|
||
|
}
|
||
|
real dUdY(size_t ix, size_t iy) const
|
||
|
{
|
||
|
if(IsCoast(ix, iy)) return fillval;
|
||
|
if(iy == 0) return (Ud(ix, 1) - Ud(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1));
|
||
|
if(iy == ny - 1) return (Ud(ix, ny - 1) - Ud(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1));
|
||
|
return 0.5 * ((Ud(ix, iy + 1) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
|
||
|
(Ud(ix, iy) - Ud(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
|
||
|
}
|
||
|
real dVdY(size_t ix, size_t iy) const
|
||
|
{
|
||
|
if(IsCoast(ix, iy)) return fillval;
|
||
|
if(iy == 0) return (Vd(ix, 1) - Vd(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1));
|
||
|
if(iy == ny - 1) return (Vd(ix, ny - 1) - Vd(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1));
|
||
|
return 0.5 * ((Vd(ix, iy + 1) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
|
||
|
(Vd(ix, iy) - Vd(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
|
||
|
}
|
||
|
|
||
|
// For stationary points
|
||
|
real LonR(real x, [[maybe_unused]] real y) const { return x0 + x * xstep; }
|
||
|
real LatR([[maybe_unused]] real x, real y) const { return y0 + y * ystep; }
|
||
|
|
||
|
public:
|
||
|
enum STPOINT
|
||
|
{
|
||
|
SADDLE = 0,
|
||
|
SACICFOCUS,
|
||
|
SKNOT,
|
||
|
UACICFOCUS,
|
||
|
UKNOT,
|
||
|
SCICFOCUS,
|
||
|
UCICFOCUS,
|
||
|
NOPOINT
|
||
|
};
|
||
|
|
||
|
struct StPoint
|
||
|
{
|
||
|
real x = 0.0, y = 0.0;
|
||
|
STPOINT type = NOPOINT;
|
||
|
};
|
||
|
|
||
|
UVData() = default;
|
||
|
UVData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_, Metric m_ = Metric::SPHERIC):
|
||
|
x0(x0_), y0(y0_), nx(nx_), ny(ny_), xstep(xs_), ystep(ys_), u(nx_ * ny_), v(nx_ * ny_), metric(m_)
|
||
|
{
|
||
|
}
|
||
|
|
||
|
const real& U(size_t i) const { return u[i]; }
|
||
|
const real& V(size_t i) const { return v[i]; }
|
||
|
const real& U(size_t ix, size_t iy) const { return U(iy * nx + ix); }
|
||
|
const real& V(size_t ix, size_t iy) const { return V(iy * nx + ix); }
|
||
|
|
||
|
real& U(size_t i) { return u[i]; }
|
||
|
real& V(size_t i) { return v[i]; }
|
||
|
real& U(size_t ix, size_t iy) { return U(iy * nx + ix); }
|
||
|
real& V(size_t ix, size_t iy) { return V(iy * nx + ix); }
|
||
|
|
||
|
size_t N() const { return u.size(); }
|
||
|
size_t Nx() const { return nx; }
|
||
|
size_t Ny() const { return ny; }
|
||
|
|
||
|
real Lon(size_t ix, [[maybe_unused]] size_t iy) const { return x0 + ix * xstep; }
|
||
|
real Lat([[maybe_unused]] size_t ix, size_t iy) const { return y0 + iy * ystep; }
|
||
|
|
||
|
real Lon(size_t i) const { return Lon(i % nx, i / nx); }
|
||
|
real Lat(size_t i) const { return Lat(i % nx, i / nx); }
|
||
|
|
||
|
explicit operator bool() const { return N() != 0; }
|
||
|
|
||
|
bool IsCoast(size_t i) const { return U(i) == fillval || V(i) == fillval; }
|
||
|
bool IsCoast(size_t ix, size_t iy) const { return U(ix, iy) == fillval || V(ix, iy) == fillval; }
|
||
|
|
||
|
static real Fillval() { return fillval; }
|
||
|
|
||
|
real Div(size_t i) const { return Div(i % nx, i / nx); }
|
||
|
real Rot(size_t i) const { return Rot(i % nx, i / nx); }
|
||
|
real Div(size_t ix, size_t iy) const
|
||
|
{
|
||
|
if(!*this) return 0.0;
|
||
|
real ux = dUdX(ix, iy);
|
||
|
real vy = dVdY(ix, iy);
|
||
|
return (ux == fillval || vy == fillval) ? fillval : (ux + vy);
|
||
|
}
|
||
|
real Rot(size_t ix, size_t iy) const
|
||
|
{
|
||
|
if(!*this) return 0.0;
|
||
|
real vx = dVdX(ix, iy);
|
||
|
real uy = dUdY(ix, iy);
|
||
|
return (vx == fillval || uy == fillval) ? fillval : (vx - uy);
|
||
|
}
|
||
|
// Okubo-Weiss parameter
|
||
|
real OW(size_t i) const { return OW(i % nx, i / nx); }
|
||
|
real OW(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 = vx - uy;
|
||
|
return sn * sn + ss * ss - w * w;
|
||
|
}
|
||
|
|
||
|
auto StablePoints(size_t ix, size_t iy) const
|
||
|
{
|
||
|
std::vector<struct StPoint> points;
|
||
|
if(!*this) return points;
|
||
|
if(ix >= Nx() - 1 || iy >= Ny() - 1) return points;
|
||
|
if(IsCoast(ix, iy) || IsCoast(ix + 1, iy) || IsCoast(ix, iy + 1) || IsCoast(ix + 1, iy + 1)) return points;
|
||
|
|
||
|
real au = U(ix, iy) - U(ix, iy + 1) - U(ix + 1, iy) + U(ix + 1, iy + 1);
|
||
|
real bu = U(ix, iy + 1) - U(ix, iy);
|
||
|
real cu = U(ix + 1, iy) - U(ix, iy);
|
||
|
real du = U(ix, iy);
|
||
|
|
||
|
real av = V(ix, iy) - V(ix, iy + 1) - V(ix + 1, iy) + V(ix + 1, iy + 1);
|
||
|
real bv = V(ix, iy + 1) - V(ix, iy);
|
||
|
real cv = V(ix + 1, iy) - V(ix, iy);
|
||
|
real dv = V(ix, iy);
|
||
|
|
||
|
real D = au * bv - av * bu;
|
||
|
real a = (au * av * cu - au * au * cv) / D;
|
||
|
real b = (au * bv * cu - au * bu * cv + au * av * du - au * au * dv) / D;
|
||
|
real c = (au * bv * du - au * bu * dv) / D;
|
||
|
real des = b * b - 4.0 * a * c;
|
||
|
if(des < 0) return points;
|
||
|
des = Sqrt(des);
|
||
|
|
||
|
real p = av * cu - au * cv;
|
||
|
real q = av * bu - au * bv;
|
||
|
real r = av * du - au * dv;
|
||
|
real x1 = 0.5 * (-b + des) / a;
|
||
|
real x2 = 0.5 * (-b - des) / a;
|
||
|
real y1 = -(p * x1 + r) / q;
|
||
|
real y2 = -(p * x2 + r) / q;
|
||
|
|
||
|
auto PointType = [au = au, av = av, bu = bu, bv = bv, cu = cu, cv = cv](real x, real y) -> auto
|
||
|
{
|
||
|
real a = 1.0;
|
||
|
real b = -y * au - x * av - bv - cu;
|
||
|
real c = y * (au * bv - av * bu) + x * (av * cu - au * cv) + bv * cu - bu * cv;
|
||
|
real des = b * b - 4.0 * a * c;
|
||
|
if(des > 0.0)
|
||
|
{
|
||
|
if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0 && c > 0.0)
|
||
|
return UKNOT;
|
||
|
else if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0)
|
||
|
return SADDLE;
|
||
|
else
|
||
|
return SKNOT;
|
||
|
}
|
||
|
else
|
||
|
{
|
||
|
bool acyclcrit = (av * y + cv > 0.0);
|
||
|
if(b < 0.0)
|
||
|
return acyclcrit ? SACICFOCUS : SCICFOCUS;
|
||
|
else
|
||
|
return acyclcrit ? UACICFOCUS : UCICFOCUS;
|
||
|
}
|
||
|
};
|
||
|
|
||
|
if(x1 >= 0.0 && x1 < 1.0 && y1 >= 0.0 && y1 < 1.0) points.emplace_back(LonR(x1 + ix, y1 + iy), LatR(x1 + ix, y1 + iy), PointType(x1, y1));
|
||
|
if(x2 >= 0.0 && x2 < 1.0 && y2 >= 0.0 && y2 < 1.0) points.emplace_back(LonR(x2 + ix, y2 + iy), LatR(x2 + ix, y2 + iy), PointType(x2, y2));
|
||
|
return points;
|
||
|
}
|
||
|
};
|