timur/mcc
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

...

Browse Source
This commit is contained in:
Timur A. Fatkhullin
2026-03-29 22:18:46 +03:00
parent 8116658e0d
commit a4b3a32d94
2 changed files with 194 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

213
include/mcc/mcc_bsplines.h
View File

@@ -7,6 +7,8 @@
#include <FortranCInterface.h> #include <FortranCInterface.h>
extern "C" { extern "C" {
/* fitting on XY-grid */
void surfit(int* iopt, void surfit(int* iopt,
int* m, int* m,
double* x, double* x,
@@ -38,6 +40,7 @@ void surfit(int* iopt,
int* kwrk, int* kwrk,
int* ier); int* ier);
/* XY-grid */
void bispev(double* tx, void bispev(double* tx,
int* nx, int* nx,
double* ty, double* ty,
@@ -56,6 +59,7 @@ void bispev(double* tx,
int* kwrk, int* kwrk,
int* ier); int* ier);
/* XY-grid */
void parder(double* tx, void parder(double* tx,
int* nx, int* nx,
double* ty, double* ty,
@@ -76,7 +80,7 @@ void parder(double* tx,
int* kwrk, int* kwrk,
int* ier); int* ier);
/* fitting on sphere */
void sphere(int* iopt, void sphere(int* iopt,
int* m, int* m,
double* teta, double* teta,
@@ -100,8 +104,44 @@ void sphere(int* iopt,
int* iwrk, int* iwrk,
int* kwrk, int* kwrk,
int* ier); int* ier);
}
/* calculation for set of points (not grid) */
void bispeu(double* tx,
int* nx,
double* ty,
int* ny,
double* c,
int* kx,
int* ky,
double* x,
double* y,
double* z,
int* m,
double* wrk,
int* lwrk,
int* ier);
/* calculation for set of points (not grid) */
void pardeu(double* tx,
int* nx,
double* ty,
int* ny,
double* c,
int* kx,
int* ky,
int* nux,
int* nuy,
double* x,
double* y,
double* z,
int* m,
double* wrk,
int* lwrk,
int* iwrk,
int* kwrk,
int* ier);
}
namespace mcc::bsplines namespace mcc::bsplines
@@ -275,14 +315,14 @@ template <std::ranges::contiguous_range TXT,
std::ranges::contiguous_range YT, std::ranges::contiguous_range YT,
std::ranges::contiguous_range CoeffT, std::ranges::contiguous_range CoeffT,
std::ranges::contiguous_range FuncT> std::ranges::contiguous_range FuncT>
int fitpack_eval_spl2d(const TXT& tx, int fitpack_eval_spl2d_grid(const TXT& tx,
const TYT& ty, const TYT& ty,
const CoeffT& coeffs, const CoeffT& coeffs,
const XT& x, const XT& x,
const YT& y, const YT& y,
FuncT& func, FuncT& func,
int kx = 3, int kx = 3,
int ky = 3) int ky = 3)
{ {
static_assert(std::same_as<std::ranges::range_value_t<TXT>, double> && static_assert(std::same_as<std::ranges::range_value_t<TXT>, double> &&
std::same_as<std::ranges::range_value_t<TYT>, double> && std::same_as<std::ranges::range_value_t<TYT>, double> &&
@@ -334,6 +374,70 @@ int fitpack_eval_spl2d(const TXT& tx,
} }
/* for set of points */
template <std::ranges::contiguous_range TXT,
std::ranges::contiguous_range TYT,
std::ranges::contiguous_range XT,
std::ranges::contiguous_range YT,
std::ranges::contiguous_range CoeffT,
std::ranges::contiguous_range FuncT>
int fitpack_eval_spl2d(const TXT& tx,
const TYT& ty,
const CoeffT& coeffs,
const XT& x,
const YT& y,
FuncT& func,
int kx = 3,
int ky = 3)
{
static_assert(std::same_as<std::ranges::range_value_t<TXT>, double> &&
std::same_as<std::ranges::range_value_t<TYT>, double> &&
std::same_as<std::ranges::range_value_t<XT>, double> &&
std::same_as<std::ranges::range_value_t<YT>, double> &&
std::same_as<std::ranges::range_value_t<CoeffT>, double> &&
std::same_as<std::ranges::range_value_t<FuncT>, double>,
"Input ranges elements type must be double!");
if (kx < 0 || ky < 0) {
return 10;
}
if (std::ranges::size(x) != std::ranges::size(y)) {
return 10;
}
int m = std::ranges::size(x);
int ier = 0;
int ntx = std::ranges::size(tx);
int nty = std::ranges::size(ty);
auto n_coeffs = (ntx - kx - 1) * (nty - ky - 1);
if (std::ranges::size(coeffs) < n_coeffs) {
return 10;
}
if (std::ranges::size(func) < m) {
std::ranges::fill_n(std::back_inserter(func), m - std::ranges::size(func), 0.0);
}
auto tx_ptr = const_cast<double*>(std::ranges::data(tx));
auto ty_ptr = const_cast<double*>(std::ranges::data(ty));
auto coeffs_ptr = const_cast<double*>(std::ranges::data(coeffs));
auto x_ptr = const_cast<double*>(std::ranges::data(x));
auto y_ptr = const_cast<double*>(std::ranges::data(y));
int lwrk = kx + ky + 2;
std::vector<double> wrk(lwrk);
bispeu(tx_ptr, &ntx, ty_ptr, &nty, coeffs_ptr, &kx, &ky, x_ptr, y_ptr, std::ranges::data(func), &m, wrk.data(),
&lwrk, &ier);
return ier;
}
// scalar x,y version // scalar x,y version
template <std::ranges::contiguous_range TXT, template <std::ranges::contiguous_range TXT,
std::ranges::contiguous_range TYT, std::ranges::contiguous_range TYT,
@@ -375,16 +479,16 @@ template <std::ranges::contiguous_range TXT,
std::ranges::contiguous_range YT, std::ranges::contiguous_range YT,
std::ranges::contiguous_range CoeffT, std::ranges::contiguous_range CoeffT,
std::ranges::contiguous_range PderT> std::ranges::contiguous_range PderT>
int fitpack_parder_spl2d(const TXT& tx, int fitpack_parder_spl2d_grid(const TXT& tx,
const TYT& ty, const TYT& ty,
const CoeffT& coeffs, const CoeffT& coeffs,
const XT& x, const XT& x,
const YT& y, const YT& y,
PderT& pder, PderT& pder,
int dx, // partial derivatives order along X int dx, // partial derivatives order along X
int dy, // partial derivatives order along Y int dy, // partial derivatives order along Y
int kx = 3, int kx = 3,
int ky = 3) int ky = 3)
{ {
static_assert(std::same_as<std::ranges::range_value_t<TXT>, double> && static_assert(std::same_as<std::ranges::range_value_t<TXT>, double> &&
std::same_as<std::ranges::range_value_t<TYT>, double> && std::same_as<std::ranges::range_value_t<TYT>, double> &&
@@ -438,6 +542,75 @@ int fitpack_parder_spl2d(const TXT& tx,
return ier; return ier;
} }
template <std::ranges::contiguous_range TXT,
std::ranges::contiguous_range TYT,
std::ranges::contiguous_range XT,
std::ranges::contiguous_range YT,
std::ranges::contiguous_range CoeffT,
std::ranges::contiguous_range PderT>
int fitpack_parder_spl2d(const TXT& tx,
const TYT& ty,
const CoeffT& coeffs,
const XT& x,
const YT& y,
PderT& pder,
int dx, // partial derivatives order along X
int dy, // partial derivatives order along Y
int kx = 3,
int ky = 3)
{
static_assert(std::same_as<std::ranges::range_value_t<TXT>, double> &&
std::same_as<std::ranges::range_value_t<TYT>, double> &&
std::same_as<std::ranges::range_value_t<XT>, double> &&
std::same_as<std::ranges::range_value_t<YT>, double> &&
std::same_as<std::ranges::range_value_t<CoeffT>, double> &&
std::same_as<std::ranges::range_value_t<PderT>, double>,
"Input ranges elements type must be double!");
if (kx < 0 || ky < 0 || dx < 0 || dy < 0) {
return 10;
}
int ntx = std::ranges::size(tx);
int nty = std::ranges::size(ty);
auto n_coeffs = (ntx - kx - 1) * (nty - ky - 1);
if (std::ranges::size(coeffs) < n_coeffs) {
return 10;
}
int m = std::ranges::size(x);
if (std::ranges::size(x) != std::ranges::size(y)) {
return 10;
}
if (std::ranges::size(pder) < m) {
std::ranges::fill_n(std::back_inserter(pder), m - std::ranges::size(pder), 0.0);
}
auto tx_ptr = const_cast<double*>(std::ranges::data(tx));
auto ty_ptr = const_cast<double*>(std::ranges::data(ty));
auto coeffs_ptr = const_cast<double*>(std::ranges::data(coeffs));
auto x_ptr = const_cast<double*>(std::ranges::data(x));
auto y_ptr = const_cast<double*>(std::ranges::data(y));
int ier = 0;
int lwrk = m * (kx + 1 - dx) + m * (ky + 1 - dy) + (ntx - kx - 1) * (nty - ky - 1);
std::vector<double> wrk(lwrk);
int kwrk = 2 * m;
std::vector<int> iwrk(kwrk);
pardeu(tx_ptr, &ntx, ty_ptr, &nty, coeffs_ptr, &kx, &ky, &dx, &dy, x_ptr, y_ptr, std::ranges::data(pder), &m,
wrk.data(), &lwrk, iwrk.data(), &kwrk, &ier);
return ier;
}
// scalar version // scalar version
template <std::ranges::contiguous_range TXT, template <std::ranges::contiguous_range TXT,
std::ranges::contiguous_range TYT, std::ranges::contiguous_range TYT,

2
include/mcc/mcc_pcm_construct.h
View File

@@ -375,7 +375,7 @@ protected:
return result; return result;
} }
result.bspline_fit_err = bsplines::fitpack_sphere_fit(theta_hw, _table.hw_colon, _table.colon_res, 1.0, ty, tx, result.bspline_fit_err = bsplines::fitpack_sphere_fit(theta_hw, _table.hw_colon, _table.colat_res, 1.0, ty, tx,
pcm_data.bspline.coeffsY, resi2y); pcm_data.bspline.coeffsY, resi2y);
if (result.bspline_fit_err > 0) { if (result.bspline_fit_err > 0) {
result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT;