diff --git a/cpp/example/CMakeLists.txt b/cpp/example/CMakeLists.txt
index 51914e686..a32b53a5c 100644
--- a/cpp/example/CMakeLists.txt
+++ b/cpp/example/CMakeLists.txt
@@ -4,11 +4,14 @@ include_directories(SYSTEM ${Sopt_INCLUDE_DIRS})
 include_directories("${PROJECT_SOURCE_DIR}/cpp")
 
 add_example(gridding LIBRARIES libpurify NOTEST)
+add_example(rm_gridding_example LIBRARIES libpurify NOTEST)
 add_example(image_wproj_chirp LIBRARIES libpurify NOTEST)
+add_example(image_rm_window LIBRARIES libpurify NOTEST)
 add_example(compare_wprojection LIBRARIES libpurify NOTEST)
 add_example(generate_vis_data LIBRARIES libpurify NOTEST)
 
 add_example(plot_wkernel LIBRARIES libpurify NOTEST)
+add_example(plot_rmkernel LIBRARIES libpurify NOTEST)
 add_example(wavelet_decomposition LIBRARIES libpurify NOTEST)
 add_example(histogram_equalisation LIBRARIES libpurify NOTEST)
 
diff --git a/cpp/example/image_rm_window.cc b/cpp/example/image_rm_window.cc
new file mode 100644
index 000000000..d0ab3c49e
--- /dev/null
+++ b/cpp/example/image_rm_window.cc
@@ -0,0 +1,65 @@
+#include "purify/config.h"
+#include "purify/types.h"
+#include "purify/directories.h"
+#include "purify/logging.h"
+#include "purify/operators.h"
+#include "purify/pfitsio.h"
+#include "purify/rm_operators.h"
+#include "purify/utilities.h"
+#include <sopt/power_method.h>
+
+#include <boost/math/special_functions/sinc.hpp>
+
+using namespace purify;
+using namespace purify::notinstalled;
+
+int main(int nargs, char const **args) {
+#define ARGS_MACRO(NAME, ARGN, VALUE, TYPE) \
+  auto const NAME =                         \
+      static_cast<TYPE>((nargs > ARGN) ? std::stod(static_cast<std::string>(args[ARGN])) : VALUE);
+
+  ARGS_MACRO(width, 1, 10, t_real)
+  ARGS_MACRO(imsize, 2, 256, t_uint)
+  ARGS_MACRO(cell, 3, 2400, t_real)
+  ARGS_MACRO(Jl_max, 4, 0, t_uint)
+#undef ARGS_MACRO
+  purify::logging::initialize();
+  purify::logging::set_level("debug");
+  // Gridding example
+  auto const oversample_ratio = 2;
+  auto const power_iters = 100;
+  auto const power_tol = 1e-5;
+  auto const Ju = 4;
+  t_int const imsizex = imsize;
+  const std::string suffix = "_" + std::to_string(static_cast<int>(width)) + "_" +
+                             std::to_string(static_cast<int>(imsize)) + "_" +
+                             std::to_string(static_cast<int>(cell));
+
+  auto const kernel = "kb";
+
+  t_uint const number_of_pixels = imsizex;
+  t_uint const number_of_vis = std::floor(number_of_pixels * 2);
+  // Generating random uv(w) coverage
+  t_real const sigma_m = constant::pi / 3;
+
+  auto header =
+      pfitsio::header_params("", " ", 1, 0, 0, stokes::I, cell, cell, 0, 1, 0, 0, 0, 0, 0);
+  const auto measure_opw = measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+      Vector<t_real>::Zero(1), Vector<t_real>::Constant(1, width), Vector<t_complex>::Ones(1),
+      imsizex, cell, oversample_ratio, kernels::kernel_from_string.at(kernel), Ju, Jl_max, 1e-6,
+      1e-6);
+  Vector<t_complex> const measurements =
+      (measure_opw->adjoint() * Vector<t_complex>::Constant(1, 1)).conjugate();
+  t_uint max_pos = 0;
+  measurements.array().real().maxCoeff(&max_pos);
+  Image<t_complex> out = measurements * std::sqrt(imsizex * oversample_ratio);
+  header.fits_name = "rm_window_real" + suffix + ".fits";
+  pfitsio::write2d(out.real(), header);
+  header.fits_name = "rm_window_imag" + suffix + ".fits";
+  pfitsio::write2d(out.imag(), header);
+  Vector<t_complex> true_window = Vector<t_complex>::Zero(imsize);
+  for (t_int i = 0; i < imsize; i++)
+    true_window(i) = boost::math::sinc_pi(width * std::abs(i - (imsize) / 2.) * cell);
+  header.fits_name = "true_window_real" + suffix + ".fits";
+  pfitsio::write2d(true_window.real(), header);
+}
diff --git a/cpp/example/plot_rmkernel.cc b/cpp/example/plot_rmkernel.cc
new file mode 100644
index 000000000..d236569ed
--- /dev/null
+++ b/cpp/example/plot_rmkernel.cc
@@ -0,0 +1,96 @@
+
+#include "purify/config.h"
+
+#include <iostream>
+#include "catch.hpp"
+#include "purify/directories.h"
+
+#include "purify/types.h"
+#include "purify/logging.h"
+
+#include "purify/kernels.h"
+
+#include "purify/directories.h"
+#include "purify/pfitsio.h"
+#include "purify/rm_kernel_integration.h"
+
+using namespace purify;
+using namespace purify::notinstalled;
+
+using namespace purify;
+
+int main(int nargs, char const **args) {
+  auto const output_name = (nargs > 2) ? static_cast<std::string>(args[1]) : "kernel";
+#define ARGS_MACRO(NAME, ARGN, VALUE, TYPE) \
+  auto const NAME =                         \
+      static_cast<TYPE>((nargs > ARGN) ? std::stod(static_cast<std::string>(args[ARGN])) : VALUE);
+
+  ARGS_MACRO(dl2, 2, 10, t_real)
+  ARGS_MACRO(absolute_error, 3, 1e-2, t_real)
+  ARGS_MACRO(imsize, 4, 1024, t_uint)
+  ARGS_MACRO(cell, 5, 20, t_real)
+
+#undef ARGS_MACRO
+  purify::logging::initialize();
+  purify::logging::set_level("debug");
+  t_uint const J = 4;
+  t_int const Jl = 30;
+  const t_real oversample_ratio = 2;
+  PURIFY_LOW_LOG("image size: {}", imsize);
+  PURIFY_LOW_LOG("cell size: {}", cell);
+  const t_real du = constant::pi / (oversample_ratio * cell * imsize);
+  PURIFY_LOW_LOG("1/du: {}", 1. / du);
+  auto const normu = [=](const t_real l) -> t_real { return kernels::ft_kaiser_bessel(l, J); };
+
+  const t_uint max_evaluations = 1e8;
+  const t_real relative_error = 0;
+  t_uint e;
+  const t_real norm = std::abs(projection_kernels::exact_rm_projection_integration(
+      0, 0, du, oversample_ratio, normu, 1e9, 0, 1e-8, integration::method::p, e));
+  auto const ftkernelu = [=](const t_real l) -> t_real {
+    return kernels::ft_kaiser_bessel(l, J) / norm;
+  };
+  t_uint total = 0;
+  t_uint rtotal = 0;
+  t_int const Ju = 20;
+  t_real const upsample = 5;
+  t_int const Ju_max = std::floor(Ju * upsample);
+  const t_int Jl_max = 100;
+  Vector<t_complex> kernel = Vector<t_complex>::Zero(Ju_max * Jl_max);
+  Vector<t_real> evals = Vector<t_real>::Zero(Ju_max * Jl_max);
+  Vector<t_real> support = Vector<t_real>::Zero(Jl_max);
+  PURIFY_LOW_LOG("dlambda^2: {}", dl2);
+  PURIFY_LOW_LOG("Kernel size: {} x {}", Ju_max, Jl_max);
+  PURIFY_LOW_LOG("FoV (rad/m^2): {}", imsize * cell);
+  PURIFY_LOW_LOG("du: {}", du);
+#pragma omp parallel for collapse(2)
+  for (int i = 0; i < Jl_max; i++) {
+    for (int j = 0; j < Ju_max; j++) {
+      t_uint evaluations = 0;
+      t_uint revaluations = 0;
+      t_uint e;
+      if (j / upsample < (dl2 * i / Jl_max / du + J) / 2) {
+        kernel(j * Jl_max + i) = projection_kernels::exact_rm_projection_integration(
+            j / upsample, dl2 * i / Jl_max, du, oversample_ratio, ftkernelu, max_evaluations,
+            absolute_error, 0, integration::method::p, evaluations);
+        evals(j * Jl_max + i) = evaluations;
+#pragma omp critical(print)
+        {
+          PURIFY_LOW_LOG("u : {}, dlambda^2 : {}", j / upsample, dl2 * i / Jl_max);
+          PURIFY_LOW_LOG("Kernel value: {}", kernel(j * Jl_max + i));
+          PURIFY_LOW_LOG("Evaluations: {}", evaluations);
+        }
+        if (kernel(j * Jl_max + i) != kernel(j * Jl_max + i))
+          throw std::runtime_error("Kernel value is a nan.");
+        total += evaluations;
+      }
+    }
+  }
+  PURIFY_LOW_LOG("Total Evaluations: {}", total);
+  PURIFY_LOW_LOG("FoV (rad/m^2): {}", imsize * cell);
+  PURIFY_LOW_LOG("du: {}", du);
+  pfitsio::write2d(kernel.real(), Jl_max, Ju_max, output_name + "_real.fits");
+  pfitsio::write2d(kernel.imag(), Jl_max, Ju_max, output_name + "_imag.fits");
+  pfitsio::write2d(kernel.cwiseAbs(), Jl_max, Ju_max, output_name + ".fits");
+  pfitsio::write2d(evals, Jl_max, Ju_max, output_name + "_evals.fits");
+}
diff --git a/cpp/example/rm_gridding_example.cc b/cpp/example/rm_gridding_example.cc
new file mode 100644
index 000000000..0b5c89383
--- /dev/null
+++ b/cpp/example/rm_gridding_example.cc
@@ -0,0 +1,145 @@
+#include "purify/config.h"
+#include "purify/types.h"
+#include "purify/logging.h"
+#include "purify/pfitsio.h"
+#include "purify/rm_operators.h"
+#include "purify/utilities.h"
+#include <sopt/power_method.h>
+
+#include <sopt/credible_region.h>
+#include <sopt/imaging_padmm.h>
+#include <sopt/power_method.h>
+#include <sopt/relative_variation.h>
+#include <sopt/utilities.h>
+#include <sopt/wavelets.h>
+#include <sopt/wavelets/sara.h>
+
+using namespace purify;
+
+namespace {
+void run_admm(const Vector<t_complex> &y,
+              const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> &measure_op,
+              const sopt::LinearTransform<Vector<t_complex>> &Psi, const t_real epsilon,
+              const std::string &prefix, const t_real oversample_ratio, const t_real channels,
+              const t_real imsizex, const t_real measure_op_norm) {
+  const Vector<t_complex> dmap = (measure_op->adjoint() * y) / channels * measure_op_norm *
+                                 measure_op_norm * oversample_ratio * imsizex;  // in jy/beam
+  pfitsio::write2d(dmap.real(), prefix + "_rm_dmap_real.fits");
+  pfitsio::write2d(dmap.imag(), prefix + "_rm_dmap_imag.fits");
+
+  auto const padmm =
+      sopt::algorithm::ImagingProximalADMM<t_complex>(y)
+          .itermax(100000)
+          .gamma((Psi.adjoint() * (measure_op->adjoint() * y).eval()).cwiseAbs().maxCoeff() * 1e-3)
+          .relative_variation(1e-7)
+          .l2ball_proximal_epsilon(epsilon)
+          .tight_frame(false)
+          .l1_proximal_tolerance(1e-3)
+          .l1_proximal_nu(1.)
+          .l1_proximal_itermax(50)
+          .l1_proximal_positivity_constraint(false)
+          .l1_proximal_real_constraint(false)
+          .residual_convergence(epsilon)
+          .lagrange_update_scale(0.9)
+          .nu(1e0)
+          .Psi(Psi)
+          .Phi(*measure_op);
+
+  auto const diagnostic = padmm();
+  Image<t_complex> image = Image<t_complex>::Map(diagnostic.x.data(), 1, imsizex);
+  pfitsio::write2d(image.real(), prefix + "_sol_real.fits");
+  pfitsio::write2d(image.imag(), prefix + "_sol_imag.fits");
+  Vector<t_complex> residuals = measure_op->adjoint() * (y - ((*measure_op) * diagnostic.x)) /
+                                channels * measure_op_norm * measure_op_norm * oversample_ratio *
+                                imsizex;
+  Image<t_complex> residual_image = Image<t_complex>::Map(residuals.data(), 1, imsizex);
+  pfitsio::write2d(residual_image.real(), prefix + "_res_real.fits");
+  pfitsio::write2d(residual_image.imag(), prefix + "_res_imag.fits");
+}
+}  // namespace
+
+int main(int nargs, char const **args) {
+  purify::logging::initialize();
+  purify::logging::set_level("debug");
+  sopt::logging::initialize();
+  sopt::logging::set_level("debug");
+  /* POSSUM
+    const t_int channels = 1100;
+    const t_real start_c = 700e6;
+    const t_real end_c = 1800e6;
+    const t_real snr = 30;
+    const t_real dnu = 1e6;  // Hz
+    auto const oversample_ratio = 2;
+    auto const power_iters = 100;
+    auto const power_tol = 1e-4;
+    auto const Ju = 4;
+    auto const imsizex = 8192 * 2;
+  */
+
+  const t_int channels = 768;
+  const t_real start_c = 200.32e6;
+  const t_real end_c = 231.04e6;
+  const t_real snr = 30;
+  const t_real dnu = 40e3;  // Hz
+  auto const oversample_ratio = 2;
+  auto const power_iters = 100;
+  auto const power_tol = 1e-4;
+  auto const Ju = 4;
+  auto const imsizex = 8192 * 2;
+
+  auto const kernel = "kb";
+  const t_int Jl_max = 100;
+  const Vector<t_real> freq = Vector<t_real>::LinSpaced(channels, start_c, end_c);
+  const Vector<t_real> lambda2 = ((constant::c / (freq.array() - dnu / 2.)).square() +
+                                  (constant::c / (freq.array() + dnu / 2.)).square()) /
+                                 2;
+  const Vector<t_real> widths = ((constant::c / (freq.array() - dnu / 2.)).square() -
+                                 (constant::c / (freq.array() + dnu / 2.)).square());
+
+  const t_real cell = constant::pi / (2 * lambda2.maxCoeff());  // rad/m^2
+  const auto measure_op_stuff = sopt::algorithm::normalise_operator<Vector<t_complex>>(
+      measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+          lambda2, widths, Vector<t_complex>::Ones(lambda2.size()), imsizex, cell, oversample_ratio,
+          kernels::kernel_from_string.at(kernel), Ju, Jl_max, 1e-6, 1e-6),
+      power_iters, power_tol, Vector<t_complex>::Random(imsizex).eval());
+  const auto measure_op = std::get<2>(measure_op_stuff);
+  const t_real measure_op_norm = std::get<0>(measure_op_stuff);
+  const t_real scale = std::sqrt(imsizex * oversample_ratio);
+
+  const auto measure_op_no_correction_stuff =
+      sopt::algorithm::normalise_operator<Vector<t_complex>>(
+          measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+              lambda2, Vector<t_real>::Zero(lambda2.size()),
+              Vector<t_complex>::Ones(lambda2.size()), imsizex, cell, oversample_ratio,
+              kernels::kernel_from_string.at(kernel), Ju, Jl_max, 1e-6, 1e-6),
+          power_iters, power_tol, Vector<t_complex>::Random(imsizex).eval());
+  const auto measure_op_no_correction = std::get<2>(measure_op_no_correction_stuff);
+  const t_real measure_op_norm_no_correction = std::get<0>(measure_op_no_correction_stuff);
+
+  Vector<t_complex> input = Vector<t_complex>::Zero(imsizex);
+  for (t_int i = 0; i < 15; i++)
+    input((i + 1) * 1024 - 1) =
+        std::exp(-2. * t_complex(0., 1.) * constant::pi * static_cast<t_real>(i) / 8.);
+  pfitsio::write2d(input.real(), "rm_input_real.fits");
+  pfitsio::write2d(input.imag(), "rm_input_imag.fits");
+
+  Vector<t_complex> y = (*measure_op * input);
+  t_real const sigma = utilities::SNR_to_standard_deviation(y, snr);
+
+  // adding noise to visibilities
+  y = utilities::add_noise(y, 0., sigma);
+  auto const epsilon = utilities::calculate_l2_radius(y.size(), sigma);
+
+  sopt::wavelets::SARA const sara{std::make_tuple("Dirac", 3u)};
+  //, std::make_tuple("DB1", 3u), std::make_tuple("DB2", 3u),
+  //     std::make_tuple("DB3", 3u),   std::make_tuple("DB4", 3u), std::make_tuple("DB5", 3u),
+  //     std::make_tuple("DB6", 3u),   std::make_tuple("DB7", 3u),
+  //  std::make_tuple("DB8", 3u)};
+
+  auto const Psi = sopt::linear_transform<t_complex>(sara, imsizex, 1);
+
+  run_admm(y, measure_op, Psi, epsilon, "corrected", oversample_ratio, channels, imsizex,
+           measure_op_norm);
+  run_admm(y, measure_op_no_correction, Psi, epsilon, "not_corrected", oversample_ratio, channels,
+           imsizex, measure_op_norm_no_correction);
+}
diff --git a/cpp/purify/CMakeLists.txt b/cpp/purify/CMakeLists.txt
index a21a95efb..f5255487b 100644
--- a/cpp/purify/CMakeLists.txt
+++ b/cpp/purify/CMakeLists.txt
@@ -28,12 +28,14 @@ set(HEADERS
   wide_field_utilities.h
   wkernel_integration.h
   wproj_operators.h
+  rm_operators.h
+  rm_kernel_integration.h
   "${PROJECT_BINARY_DIR}/include/purify/config.h")
 
 set(SOURCES utilities.cc pfitsio.cc
   kernels.cc wproj_utilities.cc operators.cc uvfits.cc yaml-parser.cc
   read_measurements.cc distribute.cc integration.cc wide_field_utilities.cc wkernel_integration.cc 
-  wproj_operators.cc)
+  wproj_operators.cc rm_operators.cc rm_kernel_integration.cc)
 
 if(TARGET casacore::ms)
   list(APPEND SOURCES casacore.cc)
diff --git a/cpp/purify/rm_kernel_integration.cc b/cpp/purify/rm_kernel_integration.cc
new file mode 100644
index 000000000..f1f3ee16e
--- /dev/null
+++ b/cpp/purify/rm_kernel_integration.cc
@@ -0,0 +1,30 @@
+#include "rm_kernel_integration.h"
+#include <boost/math/special_functions/sinc.hpp>
+
+namespace purify {
+
+namespace projection_kernels {
+t_complex fourier_rm_kernel(const t_real x, const t_real dlambda2, const t_real u,
+                            const t_real du) {
+  return std::exp(t_complex(0., -2 * constant::pi * u * x)) *
+         ((dlambda2 > 0) ? boost::math::sinc_pi(dlambda2 * x * constant::pi / du) : 1.);
+}
+
+t_complex exact_rm_projection_integration(const t_real u, const t_real dlambda2, const t_real du,
+                                          const t_real oversample_ratio,
+                                          const std::function<t_complex(t_real)> &ftkernelu,
+                                          const t_uint max_evaluations, const t_real absolute_error,
+                                          const t_real relative_error,
+                                          const integration::method method, t_uint &evaluations) {
+  auto const func = [&](const Vector<t_real> &x) -> t_complex {
+    evaluations++;
+    return ftkernelu(x(0)) * fourier_rm_kernel(x(0), dlambda2, u, du);
+  };
+  Vector<t_real> xmax = Vector<t_real>::Zero(1);
+  xmax(0) = oversample_ratio / 2.;
+  const Vector<t_real> xmin = -xmax;
+  return integration::integrate(xmin, xmax, func, integration::norm_type::paired, absolute_error,
+                                relative_error, max_evaluations, method);
+}
+}  // namespace projection_kernels
+}  // namespace purify
diff --git a/cpp/purify/rm_kernel_integration.h b/cpp/purify/rm_kernel_integration.h
new file mode 100644
index 000000000..07bfc3598
--- /dev/null
+++ b/cpp/purify/rm_kernel_integration.h
@@ -0,0 +1,23 @@
+#ifndef PURIFY_RM_KERNEL_INTEGRATION_H
+#define PURIFY_RM_KERNEL_INTEGRATION_H
+#include "purify/config.h"
+#include "purify/types.h"
+#include "purify/logging.h"
+
+#include "purify/integration.h"
+namespace purify {
+
+namespace projection_kernels {
+//! fourier integrand with sinc function
+t_complex fourier_rm_kernel(const t_real x, const t_real dlambda2, const t_real u, const t_real du);
+//! integrate fourier integrand to get rm gridding kernel
+t_complex exact_rm_projection_integration(const t_real u, const t_real dlambda2, const t_real du,
+                                          const t_real oversample_ratio,
+                                          const std::function<t_complex(t_real)> &ftkernelu,
+                                          const t_uint max_evaluations, const t_real absolute_error,
+                                          const t_real relative_error,
+                                          const integration::method method, t_uint &evaluations);
+}  // namespace projection_kernels
+// namespace projection_kernels
+}  // namespace purify
+#endif
diff --git a/cpp/purify/rm_operators.cc b/cpp/purify/rm_operators.cc
new file mode 100644
index 000000000..3260bc6d3
--- /dev/null
+++ b/cpp/purify/rm_operators.cc
@@ -0,0 +1,101 @@
+#include "purify/rm_operators.h"
+#include "purify/rm_kernel_integration.h"
+namespace purify {
+
+namespace rm_details {
+t_real pixel_to_lambda2(const t_real cell, const t_uint imsize, const t_real oversample_ratio) {
+  return constant::pi / (oversample_ratio * cell * imsize);
+}
+//! Construct gridding matrix
+Sparse<t_complex> init_gridding_matrix_1d(const Vector<t_real> &u, const Vector<t_complex> &weights,
+                                          const t_uint imsizex_, const t_real oversample_ratio,
+                                          const std::function<t_real(t_real)> kernelu,
+                                          const t_uint Ju) {
+  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
+  const t_uint rows = u.size();
+  const t_uint cols = ftsizeu_;
+
+  Sparse<t_complex> interpolation_matrix(rows, cols);
+  interpolation_matrix.reserve(Vector<t_int>::Constant(rows, Ju));
+
+  const t_complex I(0, 1);
+  const t_int ju_max = std::min(Ju, ftsizeu_);
+  // If I collapse this for loop there is a crash when using MPI... Sparse<>::insert() doesn't work
+  // right
+#pragma omp parallel for
+  for (t_int m = 0; m < rows; ++m) {
+    for (t_int ju = 1; ju < ju_max + 1; ++ju) {
+      const t_real k_u = std::floor(u(m) - ju_max * 0.5);
+      const t_uint index = utilities::mod(k_u + ju, ftsizeu_);
+      interpolation_matrix.insert(m, index) = std::exp(-2 * constant::pi * I * (k_u + ju) * 0.5) *
+                                              kernelu(u(m) - (k_u + ju)) * weights(m);
+    }
+  }
+  return interpolation_matrix;
+}
+Sparse<t_complex> init_gridding_matrix_1d(const Vector<t_real> &u, const Vector<t_real> &widths,
+                                          const Vector<t_complex> &weights, const t_uint imsizex_,
+                                          const t_real cell, const t_real oversample_ratio,
+                                          const std::function<t_real(t_real)> ftkernelu,
+                                          const t_uint Ju, const t_uint J_max,
+                                          const t_real rel_error, const t_real abs_error) {
+  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
+  const t_int rows = u.size();
+  const t_int cols = ftsizeu_;
+
+  const t_real du = rm_details::pixel_to_lambda2(cell, imsizex_, oversample_ratio);
+
+  // count gridding coefficients with variable support size
+  Vector<t_int> total_coeffs = Vector<t_int>::Zero(widths.size());
+  for (int i = 0; i < widths.size(); i++) {
+    const t_int Ju_max = std::min<t_int>(std::floor(Ju + widths(i) / du), J_max);
+    total_coeffs(i) = Ju_max;
+  }
+  const t_real num_of_coeffs = total_coeffs.array().cast<t_real>().sum();
+  PURIFY_HIGH_LOG("Max channel width in wavelengths squared: {}", widths.maxCoeff());
+  PURIFY_HIGH_LOG("Using {} rows (coefficients per a row {}), and memory of {} MegaBytes",
+                  total_coeffs.size(), total_coeffs.array().cast<t_real>().mean(),
+                  16. * num_of_coeffs / std::pow(10., 6));
+  Sparse<t_complex> interpolation_matrix(rows, cols);
+  try {
+    interpolation_matrix.reserve(total_coeffs);
+  } catch (std::bad_alloc e) {
+    throw std::runtime_error(
+        "Not enough memory for coefficients, choose upper limit on support size Jw.");
+  }
+  interpolation_matrix.reserve(Vector<t_int>::Constant(rows, Ju));
+
+  const t_complex I(0, 1);
+  const t_uint max_evaluations = 1e9;
+#pragma omp parallel for
+  for (t_int m = 0; m < rows; ++m) {
+    const t_int ju_max = std::min<t_int>(std::floor(Ju + widths(m) / du), J_max);
+    for (t_int ju = 1; ju < ju_max + 1; ++ju) {
+      const t_real k_u = std::floor(u(m) - ju_max * 0.5);
+      const t_uint index = utilities::mod(k_u + ju, ftsizeu_);
+      t_uint evaluations = 0;
+      interpolation_matrix.insert(m, index) =
+          std::exp(-2 * constant::pi * I * (k_u + ju) * 0.5) *
+          projection_kernels::exact_rm_projection_integration(
+              u(m) - (k_u + ju), widths(m), du, oversample_ratio, ftkernelu, max_evaluations,
+              abs_error, rel_error, integration::method::p, evaluations) *
+          weights(m);
+    }
+  }
+  return interpolation_matrix;
+}
+
+Image<t_complex> init_correction1d(const t_real oversample_ratio, const t_uint imsizex_,
+                                   const std::function<t_real(t_real)> ftkernelu) {
+  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
+  const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5);
+
+  Array<t_real> range;
+  range.setLinSpaced(ftsizeu_, 0.5, ftsizeu_ - 0.5);
+  return (1e0 / range.segment(x_start, imsizex_).unaryExpr([&](t_real x) {
+           return ftkernelu(x / static_cast<t_real>(ftsizeu_) - 0.5);
+         })) *
+         t_complex(1., 0.);
+}
+}  // namespace rm_details
+}  // namespace purify
diff --git a/cpp/purify/rm_operators.h b/cpp/purify/rm_operators.h
new file mode 100644
index 000000000..1115b2338
--- /dev/null
+++ b/cpp/purify/rm_operators.h
@@ -0,0 +1,269 @@
+#ifndef PURIFY_RM_OPERATOR_H
+#define PURIFY_RM_OPERATOR_H
+
+#include "purify/config.h"
+#include "purify/types.h"
+
+#include "purify/operators.h"
+
+namespace purify {
+
+namespace rm_details {
+
+//! Construct gridding matrix
+Sparse<t_complex> init_gridding_matrix_1d(const Vector<t_real> &u, const Vector<t_complex> &weights,
+                                          const t_uint imsizex_, const t_real oversample_ratio,
+                                          const std::function<t_real(t_real)> kernelu,
+                                          const t_uint Ju);
+//! Construct gridding matrix with fde
+Sparse<t_complex> init_gridding_matrix_1d(const Vector<t_real> &u, const Vector<t_real> &widths,
+                                          const Vector<t_complex> &weights, const t_uint imsizex_,
+                                          const t_real cell, const t_real oversample_ratio,
+                                          const std::function<t_real(t_real)> ftkernelu,
+                                          const t_uint Ju, const t_uint J_max,
+                                          const t_real rel_error, const t_real abs_error);
+
+//! Given the Fourier transform of a gridding kernel, creates the scaling image for gridding
+//! correction.
+Image<t_complex> init_correction1d(const t_real oversample_ratio, const t_uint imsizex_,
+                                   const std::function<t_real(t_real)> ftkernelu);
+//! gives size of lambda^2 pixel given faraday cell size (in rad/m^2) and image size
+t_real pixel_to_lambda2(const t_real cell, const t_uint imsize, const t_real oversample_ratio);
+}  // namespace rm_details
+
+namespace operators {
+
+//! constructs lambdas that apply degridding matrix with adjoint
+template <class T, class... ARGS>
+std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T>> init_gridding_matrix_1d(
+    ARGS &&... args) {
+  const std::shared_ptr<const Sparse<t_complex>> interpolation_matrix =
+      std::make_shared<const Sparse<t_complex>>(
+          rm_details::init_gridding_matrix_1d(std::forward<ARGS>(args)...));
+  const std::shared_ptr<const Sparse<t_complex>> adjoint =
+      std::make_shared<const Sparse<t_complex>>(interpolation_matrix->adjoint());
+
+  return std::make_tuple(
+      [=](T &output, const T &input) {
+        output = utilities::sparse_multiply_matrix(*interpolation_matrix, input);
+      },
+      [=](T &output, const T &input) {
+        output = utilities::sparse_multiply_matrix(*adjoint, input);
+      });
+}
+
+//! Construsts zero padding operator
+template <class T>
+std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T>> init_zero_padding_1d(
+    const Image<typename T::Scalar> &S, const t_real oversample_ratio) {
+  const t_uint imsizex_ = S.size();
+  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
+  const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5);
+  auto direct = [=](T &output, const T &x) {
+    assert(x.size() == imsizex_);
+    output = Vector<t_complex>::Zero(ftsizeu_);
+#pragma omp parallel for
+    for (t_uint j = 0; j < imsizex_; j++) output(x_start + j) = S(j) * x(j);
+  };
+  auto indirect = [=](T &output, const T &x) {
+    assert(x.size() == ftsizeu_);
+    output = T::Zero(imsizex_);
+#pragma omp parallel for
+    for (t_uint j = 0; j < imsizex_; j++) output(j) = std::conj(S(j)) * x(x_start + j);
+  };
+  return std::make_tuple(direct, indirect);
+}  // namespace operators
+
+//! Construsts FFT operator
+template <class T>
+std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T>> init_FFT_1d(
+    const t_uint imsizex_, const t_real oversample_factor_,
+    const fftw_plan fftw_plan_flag_ = fftw_plan::measure) {
+  t_int const ftsizeu_ = std::floor(imsizex_ * oversample_factor_);
+  t_int plan_flag = (FFTW_MEASURE | FFTW_PRESERVE_INPUT);
+  switch (fftw_plan_flag_) {
+  case (fftw_plan::measure):
+    plan_flag = (FFTW_MEASURE | FFTW_PRESERVE_INPUT);
+    break;
+  case (fftw_plan::estimate):
+    plan_flag = (FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
+    break;
+  }
+
+#ifdef PURIFY_OPENMP_FFTW
+  PURIFY_LOW_LOG("Using OpenMP threading with FFTW.");
+  fftw_init_threads();
+#endif
+  Vector<typename T::Scalar> src = Vector<t_complex>::Zero(ftsizeu_);
+  Vector<typename T::Scalar> dst = Vector<t_complex>::Zero(ftsizeu_);
+  // creating plans
+  const auto del = [](fftw_plan_s *plan) { fftw_destroy_plan(plan); };
+  // fftw plan with threads needs to be used before each fftw_plan is created
+#ifdef PURIFY_OPENMP_FFTW
+  fftw_plan_with_nthreads(omp_get_max_threads());
+#endif
+  const std::shared_ptr<fftw_plan_s> m_plan_forward(
+      fftw_plan_dft_1d(ftsizeu_, reinterpret_cast<fftw_complex *>(src.data()),
+                       reinterpret_cast<fftw_complex *>(dst.data()), FFTW_FORWARD, plan_flag),
+      del);
+  // fftw plan with threads needs to be used before each fftw_plan is created
+#ifdef PURIFY_OPENMP_FFTW
+  fftw_plan_with_nthreads(omp_get_max_threads());
+#endif
+  const std::shared_ptr<fftw_plan_s> m_plan_inverse(
+      fftw_plan_dft_1d(ftsizeu_, reinterpret_cast<fftw_complex *>(src.data()),
+                       reinterpret_cast<fftw_complex *>(dst.data()), FFTW_BACKWARD, plan_flag),
+      del);
+  auto const direct = [m_plan_forward, ftsizeu_](T &output, const T &input) {
+    assert(input.size() == ftsizeu_);
+    output = Vector<typename T::Scalar>::Zero(input.size());
+    fftw_execute_dft(
+        m_plan_forward.get(),
+        const_cast<fftw_complex *>(reinterpret_cast<const fftw_complex *>(input.data())),
+        reinterpret_cast<fftw_complex *>(output.data()));
+    output /= std::sqrt(output.size());
+  };
+  auto const indirect = [m_plan_inverse, ftsizeu_](T &output, const T &input) {
+    assert(input.size() == ftsizeu_);
+    output = Vector<typename T::Scalar>::Zero(input.size());
+    fftw_execute_dft(
+        m_plan_inverse.get(),
+        const_cast<fftw_complex *>(reinterpret_cast<const fftw_complex *>(input.data())),
+        reinterpret_cast<fftw_complex *>(output.data()));
+    output /= std::sqrt(output.size());
+  };
+  return std::make_tuple(direct, indirect);
+}
+
+template <class T>
+std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T>> base_padding_and_FFT_1d(
+    const std::function<t_real(t_real)> ftkernelu, const t_uint imsizex,
+    const t_real oversample_ratio = 2, const fftw_plan &ft_plan = fftw_plan::measure) {
+  sopt::OperatorFunction<T> directZ, indirectZ;
+  sopt::OperatorFunction<T> directFFT, indirectFFT;
+
+  const Image<t_complex> S =
+      purify::rm_details::init_correction1d(oversample_ratio, imsizex, ftkernelu);
+  PURIFY_LOW_LOG("Building Measurement Operator: WGFZDB");
+  PURIFY_LOW_LOG(
+      "Constructing Zero Padding "
+      "and Correction Operator: "
+      "ZDB");
+  PURIFY_MEDIUM_LOG("Faraday Depth size (width): {}", imsizex);
+  PURIFY_MEDIUM_LOG("Oversampling Factor: {}", oversample_ratio);
+  std::tie(directZ, indirectZ) = purify::operators::init_zero_padding_1d<T>(S, oversample_ratio);
+  PURIFY_LOW_LOG("Constructing FFT operator: F");
+  switch (ft_plan) {
+  case fftw_plan::measure:
+    PURIFY_MEDIUM_LOG("Measuring Plans...");
+    break;
+  case fftw_plan::estimate:
+    PURIFY_MEDIUM_LOG("Estimating Plans...");
+    break;
+  }
+  std::tie(directFFT, indirectFFT) =
+      purify::operators::init_FFT_1d<T>(imsizex, oversample_ratio, ft_plan);
+  auto direct = sopt::chained_operators<T>(directFFT, directZ);
+  auto indirect = sopt::chained_operators<T>(indirectZ, indirectFFT);
+  return std::make_tuple(direct, indirect);
+}
+
+template <class T>
+std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T>> base_degrid_operator_1d(
+    const Vector<t_real> &u, const Vector<t_complex> &weights, const t_uint imsizex,
+    const t_real oversample_ratio, const kernels::kernel kernel, const t_uint Ju,
+    const fftw_plan ft_plan) {
+  std::function<t_real(t_real)> kernelu, ftkernelu;
+  std::tie(ftkernelu, kernelu) = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio);
+  sopt::OperatorFunction<T> directFZ, indirectFZ;
+  std::tie(directFZ, indirectFZ) =
+      base_padding_and_FFT_1d<T>(ftkernelu, imsizex, oversample_ratio, ft_plan);
+  sopt::OperatorFunction<T> directG, indirectG;
+  PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG");
+  PURIFY_MEDIUM_LOG("Number of channels: {}", u.size());
+  std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_1d<T>(
+      u, weights, imsizex, oversample_ratio, kernelu, Ju);
+  auto direct = sopt::chained_operators<T>(directG, directFZ);
+  auto indirect = sopt::chained_operators<T>(indirectFZ, indirectG);
+  PURIFY_LOW_LOG("Finished consturction of RM Φ.");
+  return std::make_tuple(direct, indirect);
+}
+
+template <class T>
+std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T>> base_degrid_operator_1d(
+    const Vector<t_real> &u, const Vector<t_real> &widths, const Vector<t_complex> &weights,
+    const t_uint imsizex, const t_real cell, const t_real oversample_ratio,
+    const kernels::kernel kernel, const t_uint Ju, const t_uint J_max, const t_real abs_error,
+    const t_real rel_error, const fftw_plan ft_plan) {
+  std::function<t_real(t_real)> kernelu, ftkernelu;
+  std::tie(ftkernelu, kernelu) = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio);
+  sopt::OperatorFunction<T> directFZ, indirectFZ;
+  std::tie(directFZ, indirectFZ) =
+      base_padding_and_FFT_1d<T>(ftkernelu, imsizex, oversample_ratio, ft_plan);
+  sopt::OperatorFunction<T> directG, indirectG;
+  PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG");
+  PURIFY_MEDIUM_LOG("Number of channels: {}", u.size());
+  std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_1d<T>(
+      u, widths, weights, imsizex, cell, oversample_ratio, ftkernelu, Ju, J_max, abs_error,
+      rel_error);
+  auto direct = sopt::chained_operators<T>(directG, directFZ);
+  auto indirect = sopt::chained_operators<T>(indirectFZ, indirectG);
+  PURIFY_LOW_LOG("Finished consturction of RM Φ.");
+  return std::make_tuple(direct, indirect);
+}
+
+}  // namespace operators
+
+namespace measurementoperator {
+
+//! Returns linear transform that is the standard degridding operator
+template <class T>
+std::shared_ptr<sopt::LinearTransform<T>> init_degrid_operator_1d(
+    const Vector<t_real> &u, const Vector<t_real> &widths, const Vector<t_complex> &weights,
+    const t_uint imsizex, const t_real cell, const t_real oversample_ratio = 2,
+    const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4,
+    const t_uint J_max = 30, const t_real abs_error = 1e-6, const t_real rel_error = 1e-6) {
+  PURIFY_LOW_LOG("Cell size in Faraday Depth is {} rad/m^2", cell);
+  const t_real min_sensitivity = 1.895 / widths.maxCoeff();
+  const t_real max_sensitivity = 1.895 / widths.minCoeff();
+
+  const t_real min_null = constant::pi / widths.maxCoeff();
+  const t_real max_null = constant::pi / widths.minCoeff();
+
+  const t_real imaging_sensitivity = imsizex * cell / 2;
+  PURIFY_LOW_LOG("The max range of sensitivity is +/- {} rad/m^2", max_sensitivity);
+  PURIFY_LOW_LOG("The min range of sensitivity is +/- {} rad/m^2", min_sensitivity);
+  PURIFY_LOW_LOG("The last sensitivity null is +/- {} rad/m^2", max_null);
+  PURIFY_LOW_LOG("The first sensitivity null is +/- {} rad/m^2", min_null);
+  PURIFY_LOW_LOG("The imaging range is +/- {} rad/m^2", imsizex * cell / 2);
+  const t_real du = rm_details::pixel_to_lambda2(cell, imsizex, oversample_ratio);
+  const operators::fftw_plan ft_plan = operators::fftw_plan::measure;
+  std::array<t_int, 3> N = {0, 1, static_cast<t_int>(imsizex)};
+  std::array<t_int, 3> M = {0, 1, static_cast<t_int>(u.size())};
+  sopt::OperatorFunction<T> directDegrid, indirectDegrid;
+  std::tie(directDegrid, indirectDegrid) = purify::operators::base_degrid_operator_1d<T>(
+      -u / du, widths, weights, imsizex, cell, oversample_ratio, kernel, Ju, J_max, abs_error,
+      rel_error, ft_plan);
+  return std::make_shared<sopt::LinearTransform<T>>(directDegrid, M, indirectDegrid, N);
+}
+//! Return degridding operator with no channgel averaging
+template <class T>
+std::shared_ptr<sopt::LinearTransform<T>> init_degrid_operator_1d(
+    const Vector<t_real> &u, const Vector<t_complex> &weights, const t_uint imsizex,
+    const t_real cell, const t_real oversample_ratio = 2,
+    const kernels::kernel kernel = kernels::kernel::kb, const t_uint Ju = 4,
+    const t_uint J_max = 30, const t_real abs_error = 1e-6, const t_real rel_error = 1e-6) {
+  PURIFY_LOW_LOG("Cell size in Faraday Depth is {} rad/m^2", cell);
+  const t_real du = rm_details::pixel_to_lambda2(cell, imsizex, oversample_ratio);
+  const operators::fftw_plan ft_plan = operators::fftw_plan::measure;
+  std::array<t_int, 3> N = {0, 1, static_cast<t_int>(imsizex)};
+  std::array<t_int, 3> M = {0, 1, static_cast<t_int>(u.size())};
+  sopt::OperatorFunction<T> directDegrid, indirectDegrid;
+  std::tie(directDegrid, indirectDegrid) = purify::operators::base_degrid_operator_1d<T>(
+      -u / du, weights, imsizex, oversample_ratio, kernel, Ju, ft_plan);
+  return std::make_shared<sopt::LinearTransform<T>>(directDegrid, M, indirectDegrid, N);
+}
+
+}  // namespace measurementoperator
+};  // namespace purify
+#endif
diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt
index 1aa2c424d..b43b0ed84 100644
--- a/cpp/tests/CMakeLists.txt
+++ b/cpp/tests/CMakeLists.txt
@@ -19,6 +19,7 @@ include_directories("${PROJECT_SOURCE_DIR}/cpp" "${CMAKE_CURRENT_BINARY_DIR}/inc
 file(MAKE_DIRECTORY "${PROJECT_BINARY_DIR}/outputs")
 
 add_catch_test(measurement_operator LIBRARIES libpurify)
+add_catch_test(rm_operator LIBRARIES libpurify)
 add_catch_test(purify_fitsio LIBRARIES libpurify)
 add_catch_test(distribute LIBRARIES libpurify)
 add_catch_test(utils LIBRARIES libpurify)
diff --git a/cpp/tests/rm_operator.cc b/cpp/tests/rm_operator.cc
new file mode 100644
index 000000000..aa02c13bd
--- /dev/null
+++ b/cpp/tests/rm_operator.cc
@@ -0,0 +1,257 @@
+#include <iomanip>
+#include "catch.hpp"
+
+#include "purify/directories.h"
+#include "purify/kernels.h"
+#include "purify/pfitsio.h"
+#include "purify/rm_operators.h"
+#include "purify/test_data.h"
+#include "purify/utilities.h"
+#include <sopt/power_method.h>
+
+using namespace purify::notinstalled;
+using namespace purify;
+
+TEST_CASE("1d zeropad") {
+  const t_int imsize = 8;
+  const t_real oversample_ratio = 2;
+  Image<t_complex> S = Vector<t_complex>::Constant(imsize, 2);
+  auto pad = purify::operators::init_zero_padding_1d<Vector<t_complex>>(S, oversample_ratio);
+  const t_int boundary_size = static_cast<t_int>(imsize * oversample_ratio * 0.5 - imsize * 0.5);
+  Vector<t_complex> output;
+  const Vector<t_complex> input = Vector<t_complex>::Constant(imsize, 5);
+  std::get<0>(pad)(output, input);
+  CHECK(output.size() == static_cast<t_int>(imsize * oversample_ratio));
+  CAPTURE(boundary_size);
+  CAPTURE(input);
+  CAPTURE(output);
+  CHECK(output.head(boundary_size).isApprox(Vector<t_complex>::Zero(boundary_size)));
+  CHECK(output.tail(boundary_size).isApprox(Vector<t_complex>::Zero(boundary_size)));
+  CHECK(output.segment(boundary_size, imsize).isApprox(Vector<t_complex>::Constant(imsize, 10)));
+  // apply adjoint
+  Vector<t_complex> b_output;
+  std::get<1>(pad)(b_output, output);
+  CHECK(b_output.size() == static_cast<t_int>(imsize));
+  CHECK(b_output.isApprox(Vector<t_complex>::Constant(imsize, 20)));
+}
+
+TEST_CASE("1d FFT") {
+  const t_int imsize = 128;
+  const t_real oversample_factor = 2;
+  const operators::fftw_plan fftw_plan_flag_ = operators::fftw_plan::measure;
+  auto fft = operators::init_FFT_1d<Vector<t_complex>>(imsize, oversample_factor, fftw_plan_flag_);
+  SECTION("forward") {
+    Vector<t_complex> input = Vector<t_complex>::Zero(imsize * oversample_factor);
+    input(static_cast<t_int>(0)) = 1.;
+    Vector<t_complex> output;
+    std::get<0>(fft)(output, input);
+    const Vector<t_complex> expected =
+        Vector<t_complex>::Ones(imsize * oversample_factor) / std::sqrt(imsize * oversample_factor);
+    CAPTURE(output.head(10));
+    CAPTURE(expected.head(10));
+    CHECK(output.size() == input.size());
+    CHECK(output.isApprox(expected, 1e-6));
+  }
+  SECTION("backward") {
+    Vector<t_complex> input = Vector<t_complex>::Zero(imsize * oversample_factor);
+    input(static_cast<t_int>(0)) = 1.;
+    Vector<t_complex> output;
+    std::get<1>(fft)(output, input);
+    const Vector<t_complex> expected =
+        Vector<t_complex>::Ones(imsize * oversample_factor) / std::sqrt(imsize * oversample_factor);
+    CAPTURE(output.head(10));
+    CAPTURE(expected.head(10));
+    CHECK(output.size() == input.size());
+    CHECK(output.isApprox(expected, 1e-6));
+  }
+}
+
+TEST_CASE("regression_degrid") {
+  const t_int imsizex = 256;
+  const t_real cell = constant::pi / imsizex;
+  const Vector<t_real> u = Vector<t_real>::Random(10) * imsizex / 2;
+  const Vector<t_real> widths = Vector<t_real>::Zero(u.size());
+  Vector<t_complex> input = Vector<t_complex>::Zero(imsizex);
+  input(static_cast<t_int>(imsizex * 0.5)) = 1.;
+  const t_uint M = u.size();
+  const t_real oversample_ratio = 2;
+  const t_int ftsizeu = std::floor(imsizex * oversample_ratio);
+  const t_uint Ju = 8;
+
+  const Vector<t_complex> y = Vector<t_complex>::Ones(u.size());
+  CAPTURE(u);
+
+  SECTION("kb") {
+    const kernels::kernel kernel = kernels::kernel::kb;
+    const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
+        measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+            u, widths, Vector<t_complex>::Ones(M), imsizex, cell, oversample_ratio, kernel, Ju, 30,
+            1e-6, 1e-6);
+    const Vector<t_complex> y_test = *measure_op * input;
+    CAPTURE(y_test.array() / y.array());
+    const t_real max_test = y_test.cwiseAbs().mean();
+    CAPTURE(y_test / max_test);
+    CAPTURE(y);
+    CHECK((y_test / max_test).isApprox((y), 1e-6));
+  }
+  SECTION("pswf") {
+    const kernels::kernel kernel = kernels::kernel::pswf;
+    const auto measure_op = measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+        u, widths, Vector<t_complex>::Ones(M), imsizex, cell, oversample_ratio, kernel, 6, 30, 1e-6,
+        1e-6);
+    const Vector<t_complex> y_test = *measure_op * input;
+    CAPTURE(y_test.array() / y.array());
+    const t_real max_test = y_test.cwiseAbs().mean();
+    CAPTURE(y_test / max_test);
+    CAPTURE(y);
+    CHECK((y_test / max_test).isApprox((y), 1e-3));
+  }
+}
+
+TEST_CASE("flux units") {
+  const t_real oversample_ratio = 2;
+  Vector<t_real> u = Vector<t_real>::Random(10);
+  const Vector<t_real> widths = Vector<t_real>::Zero(u.size());
+  const t_uint M = u.size();
+  const Vector<t_complex> y = Vector<t_complex>::Ones(u.size());
+  SECTION("kb") {
+    const kernels::kernel kernel = kernels::kernel::kb;
+    for (auto& J : {4, 5, 6, 7, 8}) {
+      for (auto& imsize : {128, 256, 512}) {
+        const t_real cell = constant::pi / imsize;
+        const auto measure_op = measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+            u * imsize / 2, widths, Vector<t_complex>::Ones(M), imsize, cell, oversample_ratio,
+            kernel, J, 30, 1e-6, 1e-6);
+        Vector<t_complex> input = Vector<t_complex>::Zero(imsize);
+        input(static_cast<t_int>(imsize * 0.5)) = 1.;
+        const Vector<t_complex> y_test =
+            (*measure_op * input).eval() *
+            std::sqrt(static_cast<t_int>(input.size()) * oversample_ratio);
+        CAPTURE(y_test.cwiseAbs().mean());
+        CAPTURE(y);
+        CAPTURE(y_test);
+        CAPTURE(J);
+        CAPTURE(imsize)
+        CHECK(y_test.isApprox(y, 1e-3));
+        const Vector<t_complex> psf =
+            (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M;
+        CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5))) == Approx(1.).margin(0.001));
+      }
+    }
+  }
+  SECTION("pswf") {
+    const kernels::kernel kernel = kernels::kernel::pswf;
+    for (auto& J : {4, 5, 6, 7, 8}) {
+      for (auto& imsize : {128, 256, 512}) {
+        const t_real cell = constant::pi / imsize;
+        if (J != 6)
+          CHECK_THROWS(measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+              u * imsize / 2, widths, Vector<t_complex>::Ones(M), imsize, cell, oversample_ratio,
+              kernel, J, 30, 1e-6, 1e-6));
+        else {
+          const auto measure_op = measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+              u * imsize / 2, widths, Vector<t_complex>::Ones(M), imsize, cell, oversample_ratio,
+              kernel, J);
+          Vector<t_complex> input = Vector<t_complex>::Zero(imsize);
+          input(static_cast<t_int>(imsize * 0.5)) = 1.;
+          const Vector<t_complex> y_test =
+              (*measure_op * input).eval() *
+              std::sqrt(static_cast<t_int>(input.size()) * oversample_ratio);
+          CAPTURE(y_test.cwiseAbs().mean());
+          CAPTURE(y);
+          CAPTURE(y_test);
+          CAPTURE(J);
+          CAPTURE(imsize)
+          CHECK(y_test.isApprox(y, 1e-3));
+          const Vector<t_complex> psf =
+              (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M;
+          CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5))) == Approx(1.).margin(0.001));
+        }
+      }
+    }
+  }
+}
+
+TEST_CASE("normed operator") {
+  const t_real oversample_ratio = 2;
+  const t_int power_iters = 1000;
+  const t_real power_tol = 1e-4;
+  Vector<t_real> u = Vector<t_real>::Random(10);
+  const Vector<t_real> widths = Vector<t_real>::Zero(u.size());
+  const t_uint M = u.size();
+  const Vector<t_complex> y = Vector<t_complex>::Ones(u.size());
+  SECTION("kb") {
+    const kernels::kernel kernel = kernels::kernel::kb;
+    for (auto& J : {4, 5, 6, 7, 8}) {
+      for (auto& imsize : {128, 256, 512}) {
+        const t_real cell = constant::pi / imsize;
+        const Vector<t_complex> init = Vector<t_complex>::Ones(imsize);
+        auto power_method_result = sopt::algorithm::normalise_operator<Vector<t_complex>>(
+            measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+                u * imsize / 2, widths, Vector<t_complex>::Ones(M), imsize, cell, oversample_ratio,
+                kernel, J, 30, 1e-6, 1e-6),
+            power_iters, power_tol, init);
+        const t_real norm = std::get<0>(power_method_result);
+        const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
+            std::get<2>(power_method_result);
+
+        Vector<t_complex> input = Vector<t_complex>::Zero(imsize);
+        input(static_cast<t_int>(imsize * 0.5)) = 1.;
+        const Vector<t_complex> y_test =
+            (*measure_op * input).eval() *
+            std::sqrt(static_cast<t_int>(input.size()) * oversample_ratio);
+        CAPTURE(y_test.cwiseAbs().mean());
+        CAPTURE(y);
+        CAPTURE(y_test);
+        CAPTURE(J);
+        CAPTURE(imsize)
+        CHECK((y_test * norm).isApprox(y, 1e-3));
+        const Vector<t_complex> psf =
+            (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M * norm;
+        CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5))) == Approx(1.).margin(0.001));
+      }
+    }
+  }
+}
+
+TEST_CASE("normed operator with no widths") {
+  const t_real oversample_ratio = 2;
+  const t_int power_iters = 1000;
+  const t_real power_tol = 1e-4;
+  Vector<t_real> u = Vector<t_real>::Random(10);
+  const Vector<t_real> widths = Vector<t_real>::Zero(u.size());
+  const t_uint M = u.size();
+  const Vector<t_complex> y = Vector<t_complex>::Ones(u.size());
+  SECTION("kb") {
+    const kernels::kernel kernel = kernels::kernel::kb;
+    for (auto& J : {4, 5, 6, 7, 8}) {
+      for (auto& imsize : {128, 256, 512}) {
+        const t_real cell = constant::pi / imsize;
+        const Vector<t_complex> init = Vector<t_complex>::Ones(imsize);
+        auto power_method_result = sopt::algorithm::normalise_operator<Vector<t_complex>>(
+            measurementoperator::init_degrid_operator_1d<Vector<t_complex>>(
+                u * imsize / 2, Vector<t_complex>::Ones(M), imsize, cell, oversample_ratio, kernel,
+                J),
+            power_iters, power_tol, init);
+        const t_real norm = std::get<0>(power_method_result);
+        const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
+            std::get<2>(power_method_result);
+
+        Vector<t_complex> input = Vector<t_complex>::Zero(imsize);
+        input(static_cast<t_int>(imsize * 0.5)) = 1.;
+        const Vector<t_complex> y_test =
+            (*measure_op * input).eval() *
+            std::sqrt(static_cast<t_int>(input.size()) * oversample_ratio);
+        CAPTURE(y_test.cwiseAbs().mean());
+        CAPTURE(y);
+        CAPTURE(y_test);
+        CAPTURE(J);
+        CAPTURE(imsize)
+        CHECK((y_test * norm).isApprox(y, 1e-3));
+        const Vector<t_complex> psf =
+            (measure_op->adjoint() * y) * std::sqrt(input.size() * oversample_ratio) / M * norm;
+        CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5))) == Approx(1.).margin(0.001));
+      }
+    }
+  }
+}