diff --git a/cxx/fitpack/CMakeLists.txt b/cxx/fitpack/CMakeLists.txt
index b347d94..5820206 100644
--- a/cxx/fitpack/CMakeLists.txt
+++ b/cxx/fitpack/CMakeLists.txt
@@ -25,5 +25,6 @@ FortranCInterface_HEADER(FortranCInterface.h
 )
 
 add_library(fitpack STATIC ${src_files} fitpack.h)
+set(FITPACK_INCLUDE_DIR ${CMAKE_CURRENT_BINARY_DIR})
 
 
diff --git a/cxx/fitpack/fitpack.h b/cxx/fitpack/fitpack.h
index 7cec587..cf62cde 100644
--- a/cxx/fitpack/fitpack.h
+++ b/cxx/fitpack/fitpack.h
@@ -1,6 +1,10 @@
 #pragma once
 
+#include <expected>
+#include <limits>
+#include <memory>
 #include <ranges>
+#include <vector>
 
 namespace mcc::fitpack
 {
@@ -125,10 +129,242 @@ void sphere(int* iopt,
 }
 
 
-template <std::ranges::input_range TethaT, std::ranges::input_range PhiT>
-int fitpack_sphere(const TethaT& tetha, const PhiT& phi)
+
+template <std::ranges::contiguous_range TethaT,
+          std::ranges::contiguous_range PhiT,
+          std::ranges::contiguous_range FuncT,
+          typename WeightT,
+          std::ranges::contiguous_range TKnotT,
+          std::ranges::contiguous_range PKnotT,
+          std::ranges::contiguous_range CoeffT>
+int fitpack_sphere_smooth(const TethaT& tetha,
+                          const PhiT& phi,
+                          const FuncT& func,
+                          const WeightT& weight,
+                          double s_par,
+                          int& ntest,
+                          int& npest,
+                          TKnotT& tetha_knots,
+                          PKnotT& phi_knots,
+                          CoeffT& coeffs,
+                          double& resi2_sum,
+                          double eps = std::numeric_limits<double>::epsilon())
+    requires((std::ranges::contiguous_range<WeightT> || std::convertible_to<WeightT, double>) &&
+             std::ranges::output_range<CoeffT, double>)
 {
+    static_assert(std::same_as<std::ranges::range_value_t<TethaT>, double> &&
+                      std::same_as<std::ranges::range_value_t<PhiT>, double> &&
+                      std::same_as<std::ranges::range_value_t<FuncT>, double> &&
+                      std::same_as<std::ranges::range_value_t<TKnotT>, double> &&
+                      std::same_as<std::ranges::range_value_t<PKnotT>, double>,
+                  "Input ranges elements type must be double!");
+
+    if constexpr (std::ranges::contiguous_range<WeightT>) {
+        static_assert(std::same_as<std::ranges::range_value_t<WeightT>, double>,
+                      "Input ranges elements type must be double!");
+    }
+
+    auto m = std::min(std::min(std::ranges::size(tetha), std::ranges::size(phi)), std::ranges::size(func));
+    if (m < 2) {
+        return 10;
+    }
+
+    int res = 0;
+
+    /*   do checking here to avoid possibly unnecessary memory allocations  */
+
+    // number of knots must be >= 8 for the qubic b-splines
+    if (ntest < 8 || npest < 8) {
+        return 10;
+    }
+
+    if constexpr (std::ranges::contiguous_range<WeightT>) {
+        m = std::min(m, std::ranges::size(weight));
+        if (m < 2) {
+            return 10;
+        }
+    }
+
+    if (s_par <= 0.0) {
+        return 10;
+    }
+
+    if (eps <= 0.0 || eps >= 1.0) {
+        return 10;
+    }
+
+    const auto nt = std::ranges::size(tetha_knots);
+    const auto np = std::ranges::size(phi_knots);
+
+    if (nt < ntest) {
+        std::ranges::fill_n(std::back_inserter(tetha_knots), ntest - nt, 0.0);
+    }
+    if (np < npest) {
+        std::ranges::fill_n(std::back_inserter(phi_knots), npest - np, 0.0);
+    }
+
+
+    // compute working arrays sizes according to sphere.f
+    const int u = ntest - 7;
+    const int v = npest - 7;
+    const int uv = u * v;
+    const int v2 = v * v;
+
+    const int lwrk1 = 185 + 52 * v + 10 * u + 14 * uv + 8 * (u - 1) * v2 + 8 * m;
+    const int lwrk2 = 48 + 21 * v + 7 * uv + 4 * (u - 1) * v2;
+    const int kwrk = m + uv;
+
+    std::unique_ptr<double> wrk1{new double(lwrk1)};
+    std::unique_ptr<double> wrk2{new double(lwrk2)};
+    std::unique_ptr<int> iwrk{new int(kwrk)};
+
+    const int n_coeffs = (ntest - 4) * (npest - 4);
+    if (std::ranges::size(coeffs) < n_coeffs) {  // resize
+        std::ranges::fill_n(std::back_inserter(coeffs), n_coeffs - std::ranges::size(coeffs), 0.0);
+    }
+
+    const int iopt = 0;
+
+    if constexpr (std::ranges::contiguous_range<WeightT>) {
+        res = sphere(&iopt, &m, std::ranges::data(tetha), std::ranges::data(phi), std::ranges::data(func),
+                     std::ranges::data(weight), &s_par, &ntest, &npest, &eps, &ntest, std::ranges::data(tetha_knots),
+                     &npest, std::ranges::data(phi_knots), std::ranges::data(coeffs), &resi2_sum, wrk1.get(), &lwrk1,
+                     wrk2.get(), &lwrk2, iwrk.get(), &kwrk, &res);
+    } else {
+        std::vector<double> weight_vec(m, weight);
+
+        res = sphere(&iopt, &m, std::ranges::data(tetha), std::ranges::data(phi), std::ranges::data(func),
+                     weight_vec.data(), &s_par, &ntest, &npest, &eps, &ntest, std::ranges::data(tetha_knots), &npest,
+                     std::ranges::data(phi_knots), std::ranges::data(coeffs), &resi2_sum, wrk1.get(), &lwrk1,
+                     wrk2.get(), &lwrk2, iwrk.get(), &kwrk, &res);
+    }
+
+
+    return res;
 }
 
+/*
+    least-squares fitting (iopt=-1)
+*/
+template <std::ranges::contiguous_range TethaT,
+          std::ranges::contiguous_range PhiT,
+          std::ranges::contiguous_range FuncT,
+          typename WeightT,
+          std::ranges::contiguous_range TKnotT,
+          std::ranges::contiguous_range PKnotT,
+          std::ranges::contiguous_range CoeffT>
+int fitpack_sphere_fit(const TethaT& tetha,
+                       const PhiT& phi,
+                       const FuncT& func,
+                       const WeightT& weight,
+                       const TKnotT& tetha_knots,
+                       const PKnotT& phi_knots,
+                       CoeffT& coeffs,
+                       double& resi2_sum,
+                       double eps = std::numeric_limits<double>::epsilon())
+{
+    const double s = 100.0;
+    int nt = std::ranges::size(tetha_knots);
+    int np = std::ranges::size(phi_knots);
+
+    return fitpack_sphere_smooth(tetha, phi, func, weight, s, nt, np, tetha_knots, phi_knots, coeffs, resi2_sum, eps);
+}
+
+/*
+
+template <std::ranges::contiguous_range TethaT,
+          std::ranges::contiguous_range PhiT,
+          std::ranges::contiguous_range FuncT,
+          typename WeightT,
+          std::ranges::contiguous_range TKnotT,
+          std::ranges::contiguous_range PKnotT,
+          std::ranges::contiguous_range CoeffT>
+int fitpack_sphere_fit(const TethaT& tetha,
+                       const PhiT& phi,
+                       const FuncT& func,
+                       const WeightT& weight,
+                       const TKnotT& tetha_knots,
+                       const PKnotT& phi_knots,
+                       CoeffT& coeffs,
+                       double& resi2_sum,
+                       double eps = std::numeric_limits<double>::epsilon())
+    requires((std::ranges::contiguous_range<WeightT> || std::convertible_to<WeightT, double>) &&
+             std::ranges::output_range<CoeffT, double>)
+{
+    static_assert(std::same_as<std::ranges::range_value_t<TethaT>, double> &&
+                      std::same_as<std::ranges::range_value_t<PhiT>, double> &&
+                      std::same_as<std::ranges::range_value_t<FuncT>, double> &&
+                      std::same_as<std::ranges::range_value_t<TKnotT>, double> &&
+                      std::same_as<std::ranges::range_value_t<PKnotT>, double>,
+                  "Input ranges elements type must be double!");
+
+    if constexpr (std::ranges::contiguous_range<WeightT>) {
+        static_assert(std::same_as<std::ranges::range_value_t<WeightT>, double>,
+                      "Input ranges elements type must be double!");
+    }
+
+    auto m = std::min(std::min(std::ranges::size(tetha), std::ranges::size(phi)), std::ranges::size(func));
+    if (m < 2) {
+        return 10;
+    }
+
+    int res = 0;
+
+    const auto nt = std::ranges::size(tetha_knots);
+    const auto np = std::ranges::size(phi_knots);
+
+    // number of knots must be >= 8 for the qubic b-splines
+    if (std::min(nt, np) < 8) {
+        return 10;
+    }
+
+    if constexpr (std::ranges::contiguous_range<WeightT>) {
+        m = std::min(m, std::ranges::size(weight));
+        if (m < 2) {
+            return 10;
+        }
+    }
+
+    // compute working arrays sizes according to sphere.f
+    const int u = nt - 7;
+    const int v = np - 7;
+    const int uv = u * v;
+    const int v2 = v * v;
+
+    const int lwrk1 = 185 + 52 * v + 10 * u + 14 * uv + 8 * (u - 1) * v2 + 8 * m;
+    const int lwrk2 = 48 + 21 * v + 7 * uv + 4 * (u - 1) * v2;
+    const int kwrk = m + uv;
+
+    std::unique_ptr<double> wrk1{new double(lwrk1)};
+    std::unique_ptr<double> wrk2{new double(lwrk2)};
+    std::unique_ptr<int> iwrk{new int(kwrk)};
+
+    const int n_coeffs = (nt - 4) * (np - 4);
+    if (std::ranges::size(coeffs) < n_coeffs) {  // resize
+        std::ranges::fill_n(std::back_inserter(coeffs), n_coeffs - std::ranges::size(coeffs), 0.0);
+    }
+
+    const int iopt = -1;
+    const double s = 100.0;
+
+    if constexpr (std::ranges::contiguous_range<WeightT>) {
+        res = sphere(&iopt, &m, std::ranges::data(tetha), std::ranges::data(phi), std::ranges::data(func),
+                     std::ranges::data(weight), &s, &nt, &np, &eps, &nt, std::ranges::data(tetha_knots), &np,
+                     std::ranges::data(phi_knots), std::ranges::data(coeffs), &resi2_sum, wrk1.get(), &lwrk1,
+                     wrk2.get(), &lwrk2, iwrk.get(), &kwrk, &res);
+    } else {
+        std::vector<double> weight_vec(m, weight);
+
+        res = sphere(&iopt, &m, std::ranges::data(tetha), std::ranges::data(phi), std::ranges::data(func),
+                     weight_vec.data(), &s, &nt, &np, &eps, &nt, std::ranges::data(tetha_knots), &np,
+                     std::ranges::data(phi_knots), std::ranges::data(coeffs), &resi2_sum, wrk1.get(), &lwrk1,
+                     wrk2.get(), &lwrk2, iwrk.get(), &kwrk, &res);
+    }
+
+
+    return res;
+}
+
+*/
 
 }  // namespace mcc::fitpack