From 42a125e6961a37695539f49cd395d944863ed699 Mon Sep 17 00:00:00 2001 From: Andres Morfin Veytia <78442543+amorfinv@users.noreply.github.com> Date: Wed, 29 Jul 2026 22:09:18 +0200 Subject: [PATCH] refactor: remove cpp tools --- minisky/tools/__init__.py | 15 +- minisky/tools/src_cpp/cgeo.cpp | 698 ---------------------------- minisky/tools/src_cpp/geo.hpp | 150 ------ minisky/tools/src_cpp/pyhelpers.hpp | 98 ---- minisky/tools/src_cpp/setup.py | 10 - 5 files changed, 2 insertions(+), 969 deletions(-) delete mode 100644 minisky/tools/src_cpp/cgeo.cpp delete mode 100644 minisky/tools/src_cpp/geo.hpp delete mode 100644 minisky/tools/src_cpp/pyhelpers.hpp delete mode 100644 minisky/tools/src_cpp/setup.py diff --git a/minisky/tools/__init__.py b/minisky/tools/__init__.py index c7de0f2..059db78 100644 --- a/minisky/tools/__init__.py +++ b/minisky/tools/__init__.py @@ -1,21 +1,10 @@ """Aeronautics and geodesy tool library of MiniSky. Bundles the utility modules used throughout the simulator: unit -conversions and the ISA atmosphere (aero), geodesy functions (geo, or the -compiled cgeo variant when available), +conversions and the ISA atmosphere (aero), geodesy functions (geo), text/value converters (convert), named area shapes and inside-tests (areafilter), the navigation database (navdata), and position-text parsing (position). """ -from typing import TYPE_CHECKING - -if TYPE_CHECKING: - from . import geo as geo -else: - try: - from . import cgeo as geo # type: ignore[import-not-found] - except ImportError: - from . import geo - -from . import aero, areafilter, convert, navdata, position # noqa: E402 +from . import aero, areafilter, convert, geo, navdata, position diff --git a/minisky/tools/src_cpp/cgeo.cpp b/minisky/tools/src_cpp/cgeo.cpp deleted file mode 100644 index e63a9e3..0000000 --- a/minisky/tools/src_cpp/cgeo.cpp +++ /dev/null @@ -1,698 +0,0 @@ -#define NPY_NO_DEPRECATED_API NPY_1_10_API_VERSION -#include "Python.h" -#include "numpy/arrayobject.h" -#include "geo.hpp" -#include -#include -#define DEG2RAD 0.017453292519943295 -#define RAD2DEG 57.29577951308232 -#define M2NM 0.0005399568034557236 -#define NM2M 1852.0 - -static PyObject* cgeo_rwgs84(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL; - double lat; - if (!PyArg_ParseTuple(args, "O", &arg1)) - return NULL; - - // Check if arg is an array - if (PyArray_Check(arg1)) { - PyArrayObject *arr1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - int nd = PyArray_NDIM(arr1); - npy_intp size = PyArray_SIZE(arr1); - npy_intp* shape = PyArray_DIMS(arr1); - - PyObject* radii = PyArray_SimpleNew(nd, shape, NPY_DOUBLE); - - double* pLatd = (double*)PyArray_DATA(arr1); - double* pR = (double*)PyArray_DATA((PyArrayObject*)radii); - - while (--size >= 0) { - lat = DEG2RAD * *pLatd; - *pR = rwgs84(sin(lat), cos(lat)); - ++pLatd; ++pR; - } - - Py_DECREF(arr1); - return radii; - } else { - // arg is a scalar - lat = DEG2RAD * PyFloat_AsDouble(arg1); - return Py_BuildValue("d", rwgs84(sin(lat), cos(lat))); - } -}; - -static PyObject* cgeo_qdrdist(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - qdr_d_in ll1, ll2; - - // First check if lat1/lon1 are arrays - if (PyArray_Check(arg1) && PyArray_Check(arg2)) { - PyArrayObject* arr1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - PyArrayObject* arr2 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat1 = (double*)PyArray_DATA(arr1), - *plon1 = (double*)PyArray_DATA(arr2); - npy_intp size = PyArray_SIZE(arr1); - - // Create return vectors - PyObject* vqdr = PyArray_SimpleNew(1, &size, NPY_DOUBLE); - PyObject* vdst = PyArray_SimpleNew(1, &size, NPY_DOUBLE); - double *pqdr = (double*)PyArray_DATA((PyArrayObject*)vqdr); - double *pdst = (double*)PyArray_DATA((PyArrayObject*)vdst); - - // Check if lat2/lon2 are also arrays - if (PyArray_Check(arg3) && PyArray_Check(arg4)) { - PyArrayObject* arr3 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - PyArrayObject* arr4 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat2 = (double*)PyArray_DATA(arr3), - *plon2 = (double*)PyArray_DATA(arr4); - while (--size >= 0) { - ll1.init(DEG2RAD * *plat1, DEG2RAD * *plon1); - ll2.init(DEG2RAD * *plat2, DEG2RAD * *plon2); - *pqdr = RAD2DEG * qdr(ll1, ll2); - *pdst = M2NM * dist(ll1, ll2); - ++plat1; ++plon1; ++plat2; ++plon2; ++pqdr; ++pdst; - } - Py_DECREF(arr3); - Py_DECREF(arr4); - } else { - // lat2/lon2 are scalars - ll2.init(DEG2RAD * PyFloat_AsDouble(arg3), DEG2RAD * PyFloat_AsDouble(arg4)); - while (--size >= 0) { - ll1.init(DEG2RAD * *plat1, DEG2RAD * *plon1); - *pqdr = RAD2DEG * qdr(ll1, ll2); - *pdst = M2NM * dist(ll1, ll2); - ++plat1; ++plon1; ++pqdr; ++pdst; - } - } - Py_DECREF(arr1); - Py_DECREF(arr2); - return Py_BuildValue("NN", vqdr, vdst); - } - // Arguments should be all scalars - ll1.init(DEG2RAD * PyFloat_AsDouble(arg1), DEG2RAD * PyFloat_AsDouble(arg2)); - ll2.init(DEG2RAD * PyFloat_AsDouble(arg3), DEG2RAD * PyFloat_AsDouble(arg4)); - return Py_BuildValue("dd", RAD2DEG * qdr(ll1, ll2), M2NM * dist(ll1, ll2)); -}; - -static PyObject* cgeo_qdrdist_matrix(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - PyArrayObject *lat1 = NULL, *lon1 = NULL, *lat2 = NULL, *lon2 = NULL; - if (!PyArg_ParseTuple(args, "OO|OO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - lat1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lon1 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lat2 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lon2 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - if (lat1 == NULL || lon1 == NULL) return NULL; - - double *plat1 = (double*)PyArray_DATA(lat1), - *plon1 = (double*)PyArray_DATA(lon1); - - double *plat2 = (lat2 == NULL ? plat1 : (double*)PyArray_DATA(lat2)), - *plon2 = (lon2 == NULL ? plon1 : (double*)PyArray_DATA(lon2)); - - // Determine sizes - npy_intp size = PyArray_SIZE(lat1); - - int i = 0, j = 0; - - // Create ll2 data for efficient nested loop - std::vector ll2(size); - std::vector::iterator pll2 = ll2.begin(); - while (i < size) { - pll2->init(DEG2RAD * *plat2, DEG2RAD * *plon2); - ++i; ++plat2; ++plon2; ++pll2; - } - - // Create output matrices - npy_intp shape[] = {size, size}; - PyObject* vqdr = PyArray_SimpleNew(2, shape, NPY_DOUBLE); - PyObject* vdst = PyArray_SimpleNew(2, shape, NPY_DOUBLE); - - // Nested loop to calculate qdr and dist matrices - i = 0; - double *pqdr = (double*)PyArray_DATA((PyArrayObject*)vqdr); - double *pdst = (double*)PyArray_DATA((PyArrayObject*)vdst); - - qdr_d_in ll1; - while (i < size) { - ll1.init(DEG2RAD * *plat1, DEG2RAD * *plon1); - pll2 = ll2.begin(); - while (j < size) { - if (i == j) { - *pqdr = 0.0; - *pdst = 0.0; - } else { - *pqdr = RAD2DEG * qdr(ll1, *pll2); - *pdst = M2NM * dist(ll1, *pll2); - } - ++j; ++pll2; ++pqdr; ++pdst; - } - ++i; ++plat1; ++plon1; - j = 0; - } - //} - Py_DECREF(lat1); - Py_DECREF(lon1); - Py_XDECREF(lat2);// Py_XDECREF checks for NULL - Py_XDECREF(lon2); - - return Py_BuildValue("NN", vqdr, vdst); -}; - -static PyObject* cgeo_latlondist(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - qdr_d_in ll1, ll2; - - // First check if lat1/lon1 are arrays - if (PyArray_Check(arg1) && PyArray_Check(arg2)) { - PyArrayObject* arr1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - PyArrayObject* arr2 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat1 = (double*)PyArray_DATA(arr1), - *plon1 = (double*)PyArray_DATA(arr2); - npy_intp size = PyArray_SIZE(arr1); - - // Create return vector - PyObject* dst = PyArray_SimpleNew(1, &size, NPY_DOUBLE); - double *pdst = (double*)PyArray_DATA((PyArrayObject*)dst); - - // Check if lat2/lon2 are also arrays - if (PyArray_Check(arg3) && PyArray_Check(arg4)) { - PyArrayObject* arr3 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - PyArrayObject* arr4 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat2 = (double*)PyArray_DATA(arr3), - *plon2 = (double*)PyArray_DATA(arr4); - while (--size >= 0) { - ll1.init(DEG2RAD * *plat1, DEG2RAD * *plon1); - ll2.init(DEG2RAD * *plat2, DEG2RAD * *plon2); - *pdst = M2NM * dist(ll1, ll2); - ++plat1; ++plon1; ++plat2; ++plon2; ++pdst; - } - Py_DECREF(arr3); - Py_DECREF(arr4); - } else { - // lat2/lon2 are scalars - ll2.init(DEG2RAD * PyFloat_AsDouble(arg3), DEG2RAD * PyFloat_AsDouble(arg4)); - while (--size >= 0) { - ll1.init(DEG2RAD * *plat1, DEG2RAD * *plon1); - *pdst = M2NM * dist(ll1, ll2); - ++plat1; ++plon1; ++pdst; - } - } - Py_DECREF(arr1); - Py_DECREF(arr2); - return dst; - } - // Arguments should be all scalars - ll1.init(DEG2RAD * PyFloat_AsDouble(arg1), DEG2RAD * PyFloat_AsDouble(arg2)); - ll2.init(DEG2RAD * PyFloat_AsDouble(arg3), DEG2RAD * PyFloat_AsDouble(arg4)); - return Py_BuildValue("d", M2NM * dist(ll1, ll2)); -}; - -static PyObject* cgeo_latlondist_matrix(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - PyArrayObject *lat1 = NULL, *lon1 = NULL, *lat2 = NULL, *lon2 = NULL; - if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - lat1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lon1 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lat2 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lon2 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - if (lat1 == NULL || lon1 == NULL) return NULL; - - double *plat1 = (double*)PyArray_DATA(lat1), - *plon1 = (double*)PyArray_DATA(lon1); - - double *plat2 = (lat2 == NULL ? plat1 : (double*)PyArray_DATA(lat2)), - *plon2 = (lon2 == NULL ? plon1 : (double*)PyArray_DATA(lon2)); - - bool equal_latlon_arrays = (plat1 == plat2); - - // Determine sizes - npy_intp size = PyArray_SIZE(lat1); - - int i = 0, j = 0; - - // Create ll2 data for efficient nested loop - std::vector ll2(size); - std::vector::iterator pll2 = ll2.begin(); - while (i < size) { - pll2->init(DEG2RAD * *plat2, DEG2RAD * *plon2); - ++i; ++plat2; ++plon2; ++pll2; - } - - // Create output matrices - npy_intp shape[] = {size, size}; - PyObject* dst = PyArray_SimpleNew(2, shape, NPY_DOUBLE); - - // Nested loop to calculate dist matrix - i = 0; - double *pdst = (double*)PyArray_DATA((PyArrayObject*)dst); - if (equal_latlon_arrays) { - double *pdst_T = pdst; - std::vector::iterator pll1 = ll2.begin(); - pll2 = ll2.begin(); - while (i < size) { - while (j < size) { - if (i == j) { - *pdst = 0.0; - } else { - *pdst = *pdst_T = M2NM * dist(*pll1, *pll2); - } - ++j; ++pll2; ++pdst; - pdst_T += size; - } - ++i; ++pll1; - pdst += i; - pdst_T = pdst; - j = i; - pll2 = ll2.begin() + j; - } - } else { - qdr_d_in ll1; - while (i < size) { - ll1.init(DEG2RAD * *plat1, DEG2RAD * *plon1); - pll2 = ll2.begin(); - while (j < size) { - if (i == j) { - *pdst = 0.0; - } else { - *pdst = M2NM * dist(ll1, *pll2); - } - ++j; ++pll2; ++pdst; - } - ++i; ++plat1; ++plon1; - j = 0; - } - } - Py_DECREF(lat1); - Py_DECREF(lon1); - Py_XDECREF(lat2);// Py_XDECREF checks for NULL - Py_XDECREF(lon2); - - return Py_BuildValue("N", dst); -}; - -static PyObject* cgeo_wgsg(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL; - if (!PyArg_ParseTuple(args, "O", &arg1)) - return NULL; - - // Check if arg is an array - if (PyArray_Check(arg1)) { - PyArrayObject *arr1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - int nd = PyArray_NDIM(arr1); - npy_intp size = PyArray_SIZE(arr1); - npy_intp* shape = PyArray_DIMS(arr1); - PyObject* g = PyArray_SimpleNew(nd, shape, NPY_DOUBLE); - double* pLatd = (double*)PyArray_DATA(arr1); - double* pg = (double*)PyArray_DATA((PyArrayObject*)g); - - while (--size >= 0) { - *pg = wgsg(DEG2RAD * *pLatd); - ++pLatd; ++pg; - } - - Py_DECREF(arr1); - return g; - } else { - // arg is a scalar - return Py_BuildValue("d", wgsg(DEG2RAD * PyFloat_AsDouble(arg1))); - } -}; - -static PyObject* cgeo_qdrpos(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - // Check if args are arrays - if (PyArray_Check(arg1) && PyArray_Check(arg2) && - PyArray_Check(arg3) && PyArray_Check(arg4)) { - PyArrayObject *arr1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY), - *arr2 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY), - *arr3 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY), - *arr4 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat1 = (double*)PyArray_DATA((PyArrayObject*)arr1), - *plon1 = (double*)PyArray_DATA((PyArrayObject*)arr2), - *pqdr = (double*)PyArray_DATA((PyArrayObject*)arr3), - *pdst = (double*)PyArray_DATA((PyArrayObject*)arr4); - - // Determine sizes - npy_intp size = PyArray_SIZE(arr1); - - // Create output matrices - PyObject *lat2 = PyArray_SimpleNew(1, &size, NPY_DOUBLE), - *lon2 = PyArray_SimpleNew(1, &size, NPY_DOUBLE); - double *plat2 = (double*)PyArray_DATA((PyArrayObject*)lat2), - *plon2 = (double*)PyArray_DATA((PyArrayObject*)lon2); - - // Calculate the output vectors - pos newpos; - while (--size >= 0) { - newpos = qdrpos(DEG2RAD * *plat1, DEG2RAD * *plon1, DEG2RAD * *pqdr, NM2M * *pdst); - *plat2 = RAD2DEG * newpos.lat; - *plon2 = RAD2DEG * newpos.lon; - ++plat1; ++plon1; ++pdst; ++pqdr; ++plat2; ++plon2; - } - - Py_DECREF(arr1); - Py_DECREF(arr2); - Py_DECREF(arr3); - Py_DECREF(arr4); - - return Py_BuildValue("NN", lat2, lon2); - } else { - // Args should be scalars - pos newpos = qdrpos(DEG2RAD * PyFloat_AsDouble(arg1), DEG2RAD * PyFloat_AsDouble(arg2), - DEG2RAD * PyFloat_AsDouble(arg3), NM2M * PyFloat_AsDouble(arg4)); - return Py_BuildValue("dd", RAD2DEG * newpos.lat, RAD2DEG * newpos.lon); - } -}; - -static PyObject* cgeo_kwikdist(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - // First check if lat1/lon1 are arrays - if (PyArray_Check(arg1) && PyArray_Check(arg2)) { - PyArrayObject* arr1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - PyArrayObject* arr2 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat1 = (double*)PyArray_DATA(arr1), - *plon1 = (double*)PyArray_DATA(arr2); - npy_intp size = PyArray_SIZE(arr1); - - // Create return vector - PyObject* dst = PyArray_SimpleNew(1, &size, NPY_DOUBLE); - double *pdst = (double*)PyArray_DATA((PyArrayObject*)dst); - - // Check if lat2/lon2 are also arrays - if (PyArray_Check(arg3) && PyArray_Check(arg4)) { - PyArrayObject* arr3 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - PyArrayObject* arr4 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat2 = (double*)PyArray_DATA(arr3), - *plon2 = (double*)PyArray_DATA(arr4); - while (--size >= 0) { - *pdst = M2NM * kwikdist(kwik_in( - DEG2RAD * *plat1, DEG2RAD * *plon1, - DEG2RAD * *plat2, DEG2RAD * *plon2)); - ++plat1; ++plon1; ++plat2; ++plon2; ++pdst; - } - Py_DECREF(arr3); - Py_DECREF(arr4); - } else { - // lat2/lon2 are scalars - double lat2 = DEG2RAD * PyFloat_AsDouble(arg3), - lon2 = DEG2RAD * PyFloat_AsDouble(arg4); - while (--size >= 0) { - *pdst = *pdst = M2NM * kwikdist(kwik_in( - DEG2RAD * *plat1, DEG2RAD * *plon1, lat2, lon2)); - ++plat1; ++plon1; ++pdst; - } - } - Py_DECREF(arr1); - Py_DECREF(arr2); - return dst; - } - // Arguments should be all scalars - return Py_BuildValue("d", M2NM * kwikdist( - kwik_in(DEG2RAD * PyFloat_AsDouble(arg1), - DEG2RAD * PyFloat_AsDouble(arg2), - DEG2RAD * PyFloat_AsDouble(arg3), - DEG2RAD * PyFloat_AsDouble(arg4)))); -}; - -static PyObject* cgeo_kwikdist_matrix(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - PyArrayObject *lat1 = NULL, *lon1 = NULL, *lat2 = NULL, *lon2 = NULL; - if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - lat1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lon1 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lat2 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lon2 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - if (lat1 == NULL || lon1 == NULL) return NULL; - - double *plat1 = (double*)PyArray_DATA(lat1), - *plon1 = (double*)PyArray_DATA(lon1); - - double *plat2 = (lat2 == NULL ? plat1 : (double*)PyArray_DATA(lat2)), - *plon2 = (lon2 == NULL ? plon1 : (double*)PyArray_DATA(lon2)); - - bool equal_latlon_arrays = (plat1 == plat2); - - // Determine sizes - npy_intp size = PyArray_SIZE(lat1); - - // Create output matrices - npy_intp shape[] = {size, size}; - PyObject* dst = PyArray_SimpleNew(2, shape, NPY_DOUBLE); - double *pdst = (double*)PyArray_DATA((PyArrayObject*)dst); - // Nested loop to calculate dist matrix - int i = 0, j = 0; - if (equal_latlon_arrays) { - double *pdst_T = pdst; - while (i < size) { - while (j < size) { - if (i == j) { - *pdst = 0.0; - } else { - *pdst = *pdst_T = M2NM * kwikdist( - kwik_in(DEG2RAD * *plat1, DEG2RAD * *plon1, DEG2RAD * *plat2, DEG2RAD * *plon2)); - } - ++j; ++plat2; ++plon2; ++pdst; - pdst_T += size; - } - ++i; ++plat1; ++plon1; - pdst += i; - pdst_T = pdst; - j = i; - plat2 = (double*)PyArray_DATA(lat2) + j; - plon2 = (double*)PyArray_DATA(lon2) + j; - } - } else { - while (i < size) { - while (j < size) { - if (i == j) { - *pdst = 0.0; - } else { - *pdst = M2NM * kwikdist( - kwik_in(DEG2RAD * *plat1, DEG2RAD * *plon1, DEG2RAD * *plat2, DEG2RAD * *plon2)); - } - ++j; ++plat2; ++plon2; ++pdst; - } - ++i; ++plat1; ++plon1; - j = 0; - plat2 = (double*)PyArray_DATA(lat2); - plon2 = (double*)PyArray_DATA(lon2); - } - } - Py_DECREF(lat1); - Py_DECREF(lon1); - Py_XDECREF(lat2);// Py_XDECREF checks for NULL - Py_XDECREF(lon2); - - return Py_BuildValue("N", dst); -}; - -static PyObject* cgeo_kwikqdrdist(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - // First check if lat1/lon1 are arrays - if (PyArray_Check(arg1) && PyArray_Check(arg2)) { - PyArrayObject* arr1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - PyArrayObject* arr2 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat1 = (double*)PyArray_DATA(arr1), - *plon1 = (double*)PyArray_DATA(arr2); - npy_intp size = PyArray_SIZE(arr1); - - // Create return vectors - PyObject *qdr = PyArray_SimpleNew(1, &size, NPY_DOUBLE), - *dst = PyArray_SimpleNew(1, &size, NPY_DOUBLE); - double *pqdr = (double*)PyArray_DATA((PyArrayObject*)qdr), - *pdst = (double*)PyArray_DATA((PyArrayObject*)dst); - - // Check if lat2/lon2 are also arrays - if (PyArray_Check(arg3) && PyArray_Check(arg4)) { - PyArrayObject* arr3 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - PyArrayObject* arr4 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - double *plat2 = (double*)PyArray_DATA(arr3), - *plon2 = (double*)PyArray_DATA(arr4); - while (--size >= 0) { - kwik_in in(DEG2RAD * *plat1, DEG2RAD * *plon1, DEG2RAD * *plat2, DEG2RAD * *plon2); - *pqdr = RAD2DEG * kwikqdr(in); - *pdst = M2NM * kwikdist(in); - ++plat1; ++plon1; ++plat2; ++plon2; ++pqdr; ++pdst; - } - Py_DECREF(arr3); - Py_DECREF(arr4); - } else { - // lat2/lon2 are scalars - double lat2 = DEG2RAD * PyFloat_AsDouble(arg3), - lon2 = DEG2RAD * PyFloat_AsDouble(arg4); - while (--size >= 0) { - kwik_in in(DEG2RAD * *plat1, DEG2RAD * *plon1, lat2, lon2); - *pqdr = RAD2DEG * kwikqdr(in); - *pdst = M2NM * kwikdist(in); - ++plat1; ++plon1; ++pqdr; ++pdst; - } - } - Py_DECREF(arr1); - Py_DECREF(arr2); - return dst; - } - // Arguments should be all scalars - kwik_in in(DEG2RAD * PyFloat_AsDouble(arg1), - DEG2RAD * PyFloat_AsDouble(arg2), - DEG2RAD * PyFloat_AsDouble(arg3), - DEG2RAD * PyFloat_AsDouble(arg4)); - return Py_BuildValue("dd", RAD2DEG * kwikqdr(in), M2NM * kwikdist(in)); -}; - -static PyObject* cgeo_kwikqdrdist_matrix(PyObject* self, PyObject* args) -{ - PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL, *arg4 = NULL; - PyArrayObject *lat1 = NULL, *lon1 = NULL, *lat2 = NULL, *lon2 = NULL; - if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4)) - return NULL; - - lat1 = (PyArrayObject*)PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lon1 = (PyArrayObject*)PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lat2 = (PyArrayObject*)PyArray_FROM_OTF(arg3, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - lon2 = (PyArrayObject*)PyArray_FROM_OTF(arg4, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); - if (lat1 == NULL || lon1 == NULL) return NULL; - - double *plat1 = (double*)PyArray_DATA(lat1), - *plon1 = (double*)PyArray_DATA(lon1); - - double *plat2 = (lat2 == NULL ? plat1 : (double*)PyArray_DATA(lat2)), - *plon2 = (lon2 == NULL ? plon1 : (double*)PyArray_DATA(lon2)); - - bool equal_latlon_arrays = (plat1 == plat2); - - // Determine sizes - npy_intp size = PyArray_SIZE(lat1); - - // Create output matrices - npy_intp shape[] = {size, size}; - PyObject *qdr = PyArray_SimpleNew(2, shape, NPY_DOUBLE), - *dst = PyArray_SimpleNew(2, shape, NPY_DOUBLE); - double *pqdr = (double*)PyArray_DATA((PyArrayObject*)qdr), - *pdst = (double*)PyArray_DATA((PyArrayObject*)dst); - // Nested loop to calculate dist matrix - int i = 0, j = 0; - if (equal_latlon_arrays) { - double *pqdr_T = pqdr, - *pdst_T = pdst; - while (i < size) { - while (j < size) { - if (i == j) { - *pqdr = 0.0; - *pdst = 0.0; - } else { - kwik_in in(DEG2RAD * *plat1, DEG2RAD * *plon1, DEG2RAD * *plat2, DEG2RAD * *plon2); - *pqdr = *pqdr_T = RAD2DEG * kwikqdr(in); - *pdst = *pdst_T = M2NM * kwikdist(in); - } - ++j; ++plat2; ++plon2; ++pqdr; ++pdst; - pqdr_T += size; pdst_T += size; - } - ++i; ++plat1; ++plon1; - pqdr += i; pdst += i; - pqdr_T = pqdr; pdst_T = pdst; - j = i; - plat2 = (double*)PyArray_DATA(lat2) + j; - plon2 = (double*)PyArray_DATA(lon2) + j; - } - } else { - while (i < size) { - while (j < size) { - if (i == j) { - *pqdr = 0.0; - *pdst = 0.0; - } else { - kwik_in in(DEG2RAD * *plat1, DEG2RAD * *plon1, DEG2RAD * *plat2, DEG2RAD * *plon2); - *pqdr = RAD2DEG * kwikqdr(in); - *pdst = M2NM * kwikdist(in); - } - ++j; ++plat2; ++plon2; ++pqdr; ++pdst; - } - ++i; ++plat1; ++plon1; - j = 0; - plat2 = (double*)PyArray_DATA(lat2); - plon2 = (double*)PyArray_DATA(lon2); - } - } - Py_DECREF(lat1); - Py_DECREF(lon1); - Py_XDECREF(lat2);// Py_XDECREF checks for NULL - Py_XDECREF(lon2); - - return Py_BuildValue("NN", qdr, dst); -}; - -static struct PyMethodDef methods[] = { - {"rwgs84", cgeo_rwgs84, METH_VARARGS, "Get local earth radius using WGS'84 spec."}, - {"rwgs84_matrix", cgeo_rwgs84, METH_VARARGS, "Get local earth radius using WGS'84 spec (for vectors)."}, - {"qdrdist", cgeo_qdrdist, METH_VARARGS, "Calculate bearing and distance between lat1+lon1 and lat2+lon2"}, - {"qdrdist_matrix", cgeo_qdrdist_matrix, METH_VARARGS, "Calculate bearing and distance matrices between vectors lat1+lon1/lat2+lon2"}, - {"latlondist", cgeo_latlondist, METH_VARARGS, "Calculate distance between lat1+lon1 and lat2+lon2"}, - {"latlondist_matrix", cgeo_latlondist_matrix, METH_VARARGS, "Calculate distance matrix between vectors lat1+lon1/lat2+lon2"}, - {"wgsg", cgeo_wgsg, METH_VARARGS, "Gravity acceleration at a given latitude according to WGS'84"}, - {"qdrpos", cgeo_qdrpos, METH_VARARGS, "Calculate position from reference position, bearing and distance"}, - {"kwikdist", cgeo_kwikdist, METH_VARARGS, "Quick and dirty dist [nm]"}, - {"kwikdist_matrix", cgeo_kwikdist_matrix, METH_VARARGS, "Quick and dirty dist [nm] (for vectors)"}, - {"kwikqdrdist", cgeo_kwikqdrdist, METH_VARARGS, "Quick and dirty dist [nm] and bearing [deg]"}, - {"kwikqdrdist_matrix", cgeo_kwikqdrdist_matrix, METH_VARARGS, "Quick and dirty dist [nm] and bearing [deg] (for vectors)"}, - {NULL, NULL, 0, NULL} -}; - -#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */ -#define PyMODINIT_FUNC void -#endif - -#if PY_MAJOR_VERSION >= 3 -static struct PyModuleDef cgeodef = -{ - PyModuleDef_HEAD_INIT, - "cgeo", /* name of module */ - "", /* module documentation, may be NULL */ - -1, /* size of per-interpreter state of the module, or -1 if the module keeps state in global variables. */ - methods -}; - -PyMODINIT_FUNC PyInit_cgeo(void) -{ - import_array(); - return PyModule_Create(&cgeodef); -}; -#else -PyMODINIT_FUNC initcgeo() -{ - Py_InitModule("cgeo", methods); - import_array(); -}; -#endif diff --git a/minisky/tools/src_cpp/geo.hpp b/minisky/tools/src_cpp/geo.hpp deleted file mode 100644 index 19d73ae..0000000 --- a/minisky/tools/src_cpp/geo.hpp +++ /dev/null @@ -1,150 +0,0 @@ -#include - -// Earth major and minor axes -static const double a = 6378137.0, // [m] Major semi-axis WGS-84 - b = 6356752.314245, // [m] Minor semi-axis WGS-84 - re = 6371000.0, // [m] average earth radius - a2 = a * a, - b2 = b * b, - a4 = a2 * a2, - b4 = b2 * b2; - -// Return the distance from the Earth's center to a point on the spheroid -// surface at geodetic latitude, lat (radians). Pass cos(lat) and sin(lat) -inline double rwgs84(const double& sinlat, const double& coslat) -{ - // See https://en.wikipedia.org/wiki/Earth_radius#Geocentric_radius - double sinlat2 = sinlat * sinlat, - coslat2 = coslat * coslat; - - // Calculate and return the radius in meters: - return sqrt((a4 * coslat2 + b4 * sinlat2) / - (a2 * coslat2 + b2 * sinlat2)); -} - -struct qdr_d_in { - double lat, lon, sinlat, coslat; - void init(const double& lat, const double& lon) { - this->lat = lat; this->lon = lon; - this->sinlat = sin(lat); this->coslat = cos(lat); - }}; - -// return great-circle distance between two points. -// The implementation uses the Haversine formula: -// a = sin²(Δφ/2) + cos φ1 * cos φ2 * sin²(Δλ/2) -// c = 2 * atan2( sqrt(a), sqrt(1−a) ) -// d = R * c -// where -// φ is latitude (radians) -// λ is longitude -// R is Earth’s radius -// see http://www.movable-type.co.uk/scripts/latlong.html#ortho-dist -inline double dist(const qdr_d_in& ll1, const qdr_d_in& ll2) -{ - double sindlat2 = sin(0.5 * (ll2.lat - ll1.lat)), - sindlon2 = sin(0.5 * (ll2.lon - ll1.lon)); - - double r; - if (ll1.lat * ll2.lat >= 0.0) { - r = rwgs84(sin(0.5 * (ll1.lat + ll2.lat)), cos(0.5 * (ll1.lat + ll2.lat))); - } else { - r = 0.5 * ( fabs(ll1.lat) * (rwgs84(ll1.sinlat, ll1.coslat) + a) + - fabs(ll2.lat) * (rwgs84(ll2.sinlat, ll2.coslat) + a)) - / ( fabs(ll1.lat) + fabs(ll2.lat)); - } - - double root = sindlat2 * sindlat2 + ll1.coslat * ll2.coslat * sindlon2 * sindlon2; - return 2.0 * r * atan2(sqrt(root), sqrt(1.0 - root)); -} - -// initial bearing (radians) initial heading for a great-circle route from -// point ll1, the start point, to point ll2, the end point. -// θ = atan2( sin Δλ * cos φ2 , cos φ1 * sin φ2 − sin φ1 * cos φ2 * cos Δλ ) -// where -// φ1 is the latitude (radians) of the start point -// λ1 is the longitude (radians) of the start point -// φ2 is the latitude (radians) of the end point -// λ2 is the longitude (radians) of the end point -// Δλ is the difference in longitude -// see http://www.movable-type.co.uk/scripts/latlong.html#bearing -// see http://williams.best.vwh.net/avform.htm -inline double qdr(const qdr_d_in& ll1, const qdr_d_in& ll2) -{ - return atan2(sin(ll2.lon - ll1.lon) * ll2.coslat, - ll2.sinlat * ll1.coslat - ll1.sinlat * ll2.coslat * cos(ll2.lon - ll1.lon)); -} - -inline double wgsg(const double& lat) -{ - static const double geq = 9.7803; // m/s2 g at equator - static const double e2 = 6.694e-3; // eccentricity - static const double k = 0.001932; // derived from flattening f, 1/f = 298.257223563 - - double sinlat = sin(lat); - return geq * (1.0 + k * sinlat * sinlat) / sqrt(1.0 - e2 * sinlat * sinlat); -} - -// position on Earth, lat, lon (radians) -struct pos {double lat, lon;}; - -// Return the new position starting from point (lon1, lat1) -// with initial bearing qdr, for a distance dist -// φ2 = asin( sin φ1 * cos δ + cos φ1 * sin δ * cos θ ) -// λ2 = λ1 + atan2( sin θ * sin δ * cos φ1, cos δ − sin φ1 * sin φ2 ) -// where -// φ is latitude, -// λ is longitude, -// θ is the bearing (clockwise from north), -// δ is the angular distance d/R -// d is the distance travelled -// R is Earth’s radius -// see http://www.movable-type.co.uk/scripts/latlong.html#destPoint -inline pos qdrpos(const double& lat1, const double& lon1, const double& qdr, const double& dist) -{ - // Calculate new position - double sinlat = sin(lat1), - coslat = cos(lat1); - double R = rwgs84(sinlat, coslat); - double sdr = sin(dist / R), - cdr = cos(dist / R); - pos newpos; - newpos.lat = asin(sinlat * cdr + coslat * sdr * cos(qdr)); - - newpos.lon = lon1 + atan2(sin(qdr) * sdr * coslat, - cdr - sinlat * sin(newpos.lat)); - return newpos; -} - -inline pos kwikpos(const double& lat1, const double& lon1, const double& qdr, const double& dist) -{ - pos newpos; - newpos.lat = lat1 + cos(qdr) * dist / re; - newpos.lon = lon1 + sin(qdr) * dist / re / cos(lat1); - return newpos; -} - -struct kwik_in {double dlat, dlon, cavelat; - kwik_in(const double& lat1, const double& lon1, const double& lat2, const double& lon2) : - dlat(lat2 - lat1), dlon(lon2 - lon1), cavelat(cos(0.5 * (lat1 + lat2))) {}; -}; - -// quick distance calculation (using equirectangular distance approximation) -inline double kwikdist(const kwik_in& in) -{ - double dangle = sqrt(in.dlat * in.dlat + in.dlon * in.dlon * in.cavelat * in.cavelat); - return re * dangle; -} - -// struct posvec {double dx, dy; posvec(const double& dx, const double& dy): dx(dx), dy(dy) {}}; -// inline posvec kwikposvec(const double &lat1, const double &lon1, const double &lat2, const double &lon2) -// { -// return posvec( -// re * (lon1 - lon2) * cos(0.5 * (lat1 + lat2)), -// re * (lat1 - lat2) -// ); -// }; - -inline double kwikqdr(const kwik_in &in) -{ - return fmod(atan2(in.dlon * in.cavelat, in.dlat), 360.0); -} \ No newline at end of file diff --git a/minisky/tools/src_cpp/pyhelpers.hpp b/minisky/tools/src_cpp/pyhelpers.hpp deleted file mode 100644 index 906c5ea..0000000 --- a/minisky/tools/src_cpp/pyhelpers.hpp +++ /dev/null @@ -1,98 +0,0 @@ -#define NPY_NO_DEPRECATED_API NPY_1_10_API_VERSION -#include -#include "structmember.h" -#include "numpy/arrayobject.h" -#include - -template int atype(); -template<> int atype() {return NPY_DOUBLE;}; -template<> int atype() {return NPY_BOOL;}; - -struct PyAttr { - PyObject* attr; - PyAttr(PyObject* attr, bool refowned = false) : attr(NULL) { - if (refowned) { - this->attr = attr; - } - } - PyAttr(PyObject* parent, const char* name) : - attr(PyObject_GetAttrString(parent, name)) {} - PyAttr(const PyAttr& parent, const char* name) : - attr(PyObject_GetAttrString(parent.attr, name)) {} - ~PyAttr() {Py_XDECREF(attr);} -}; - -template -struct PyArrayAttr: public PyAttr { - T *ptr, *ptr_start; - PyArrayObject* arr; - PyArrayAttr(PyObject* parent, const char* name) : - PyAttr(parent, name) {init(attr);} - PyArrayAttr(const PyAttr& parent, const char* name) : - PyAttr(parent, name) {init(attr);} - PyArrayAttr(PyObject* attr) : PyAttr(attr) { - init(attr); - } - - PyArrayAttr(const int length) : PyAttr(NULL) { - int nd = 1; - npy_intp dims[] = {length}; - arr = (PyArrayObject*)PyArray_SimpleNew(nd, dims, T_ENUM); - if (arr != NULL) { - ptr_start = ptr = (T*)PyArray_DATA(arr); - } - } - - void init(PyObject* obj) { - if (obj != NULL) { - arr = (PyArrayObject*)PyArray_FROM_OTF(obj, atype(), NPY_ARRAY_IN_ARRAY); - if (arr != NULL) { - ptr_start = ptr = (T*)PyArray_DATA(arr); - } else { - std::cout << "Couldn't create array from python object" << std::endl; - } - } - } - ~PyArrayAttr() - { - Py_XDECREF(arr); - } - - operator bool() const {return (arr != NULL);} - npy_intp size() const {return PyArray_SIZE(arr);} -}; - -typedef PyArrayAttr PyDoubleArrayAttr; -typedef PyArrayAttr PyBoolArrayAttr; - -struct PyListAttr: public PyAttr { - PyListAttr(int size=0) : PyAttr(PyList_New(size), true) {} - PyListAttr(PyObject* attr) : PyAttr(attr) {} - PyListAttr(PyObject* parent, const char* name) : PyAttr(parent, name) {} - PyListAttr(const PyAttr& parent, const char* name) : PyAttr(parent, name) {} - PyObject* operator[](Py_ssize_t idx) const {return PyList_GetItem(attr, idx);} - inline int setItem(const Py_ssize_t& idx, const int& item) {return PyList_SetItem(attr, idx, PyLong_FromLong(item));} - inline int setItem(const Py_ssize_t& idx, const double& item) {return PyList_SetItem(attr, idx, PyFloat_FromDouble(item));} - inline int setItem(const Py_ssize_t& idx, PyObject* item) {return PyList_SetItem(attr, idx, item);} - inline int append(const int& item) { - PyObject* o = PyLong_FromLong(item); - int index = PyList_Append(attr, o); - Py_DECREF(o); - return index;} - inline int append(const double& item) { - PyObject* o = PyFloat_FromDouble(item); - int index = PyList_Append(attr, o); - Py_DECREF(o); - return index;} - inline int append(PyObject* item) {return PyList_Append(attr, item);} -}; - -double GetAttrDouble(PyObject* parent, const char* name) { - PyAttr a(parent, name); - return PyFloat_AsDouble(a.attr); -}; - -int GetAttrInt(PyObject* parent, const char* name) { - PyAttr a(parent, name); - return PyLong_AsLong(a.attr); -}; diff --git a/minisky/tools/src_cpp/setup.py b/minisky/tools/src_cpp/setup.py deleted file mode 100644 index a0ffb16..0000000 --- a/minisky/tools/src_cpp/setup.py +++ /dev/null @@ -1,10 +0,0 @@ -#!/usr/bin/env python -# -*- coding: UTF-8 -*- - -from distutils.core import Extension, setup - -import numpy as np - -ext_modules = [Extension("cgeo", sources=["cgeo.cpp"])] - -setup(name="cgeo", version="1.0", include_dirs=[np.get_include()], ext_modules=ext_modules)