diff --git a/src/bindings/pybindings.cpp b/src/bindings/pybindings.cpp
index 560e4994..0ca5accc 100644
--- a/src/bindings/pybindings.cpp
+++ b/src/bindings/pybindings.cpp
@@ -408,16 +408,24 @@ PYBIND11_MODULE(core, module) {
         .def_readwrite("use_adaptive_directions", &Rays::use_adaptive_directions,
             "Whether to use a different set of directions for each ray.")
         .def("get_direction_index", &Rays::get_direction_index<false>,
-            "Get the linearized direction index for a given point and ray, when not using an adaptive ray discretization.")
+            "Get the linearized direction index for a given point and ray, when not using an "
+            "adaptive ray discretization.")
         .def("get_direction_index_adaptive", &Rays::get_direction_index<true>,
-            "Get the linearized direction index for a given point and ray, when using the adaptive ray discretization.")
-        .def("get_direction", &Rays::get_direction<false>, "Get the direction of a ray, when not using an adaptive ray discretization.")
-        .def("get_direction_adaptive", &Rays::get_direction<true>, "Get the direction of a ray, when using the adaptive ray discretization.")
-        .def("get_antipod", &Rays::get_antipod<false>, "Get the antipodal direction of a ray, when not using an adaptive ray discretization.")
-        .def("get_antipod_adaptive", &Rays::get_antipod<true>, "Get the antipodal direction of a ray, when using the adaptive ray discretization.")
+            "Get the linearized direction index for a given point and ray, when using the adaptive "
+            "ray discretization.")
+        .def("get_direction", &Rays::get_direction<false>,
+            "Get the direction of a ray, when not using an adaptive ray discretization.")
+        .def("get_direction_adaptive", &Rays::get_direction<true>,
+            "Get the direction of a ray, when using the adaptive ray discretization.")
+        .def("get_antipod", &Rays::get_antipod<false>,
+            "Get the antipodal direction of a ray, when not using an adaptive ray discretization.")
+        .def("get_antipod_adaptive", &Rays::get_antipod<true>,
+            "Get the antipodal direction of a ray, when using the adaptive ray discretization.")
         .def("get_antipod_index", &Rays::get_antipod_index, "Get the index of the antipodal ray.")
-        .def("get_weight", &Rays::get_weight<false>, "Get the weight of a ray, when not using an adaptive ray discretization.")
-        .def("get_weight_adaptive", &Rays::get_weight<true>, "Get the weight of a ray, when using the adaptive ray discretization.")
+        .def("get_weight", &Rays::get_weight<false>,
+            "Get the weight of a ray, when not using an adaptive ray discretization.")
+        .def("get_weight_adaptive", &Rays::get_weight<true>,
+            "Get the weight of a ray, when using the adaptive ray discretization.")
         // io
         .def("read", &Rays::read, "Read object from file.")
         .def("write", &Rays::write, "Write object to file.");
diff --git a/src/model/geometry/rays/rays.cpp b/src/model/geometry/rays/rays.cpp
index 3d4a3ce1..708f97be 100644
--- a/src/model/geometry/rays/rays.cpp
+++ b/src/model/geometry/rays/rays.cpp
@@ -15,7 +15,8 @@ void Rays ::read(const Io& io) {
         parameters->set_hnrays(parameters->nrays() / 2);
         use_adaptive_directions = true;
     } else {
-        parameters->set_nrays(len_dir);//will error, and give some more info on the sizes of the arrays
+        parameters->set_nrays(
+            len_dir); // will error, and give some more info on the sizes of the arrays
     }
 
     antipod.resize(parameters->nrays());
@@ -33,7 +34,7 @@ void Rays ::read(const Io& io) {
             direction[r] =
                 Vector3D(direction_buffer[r][0], direction_buffer[r][1], direction_buffer[r][2]);
         }
-    } else {//use adaptive directions == true
+    } else { // use adaptive directions == true
         weight.resize(parameters->nrays() * parameters->npoints());
         io.read_list(prefix + "weight", weight);
 
@@ -86,48 +87,43 @@ void Rays ::write(const Io& io) const {
     io.write_list(prefix + "weight", weight);
 }
 
-// Note: I can't seem to get constexpr to work correctly (linker error), thus manual template specialization instead
-// Note2: It might be that pybind11 expects the manual specialization, so all other template functions in this file are manually specialized
+// Note: I can't seem to get constexpr to work correctly (linker error), thus manual template
+// specialization instead Note2: It might be that pybind11 expects the manual specialization, so all
+// other template functions in this file are manually specialized
 
-template <>
-Size Rays ::get_direction_index<true>(const Size pointidx, const Size rayidx) const {
+template <> Size Rays ::get_direction_index<true>(const Size pointidx, const Size rayidx) const {
     return pointidx * parameters->nrays() + rayidx;
 }
 
-template <>
-Size Rays ::get_direction_index<false>(const Size pointidx, const Size rayidx) const {
+template <> Size Rays ::get_direction_index<false>(const Size pointidx, const Size rayidx) const {
     return rayidx;
 }
 
 /// Get the direction of a ray
 /// @param[in] pointidx: Index of the point
 /// @param[in] rayidx: Index of the ray
-template <>
-Vector3D Rays ::get_direction<true>(const Size pointidx, const Size rayidx) const {
+template <> Vector3D Rays ::get_direction<true>(const Size pointidx, const Size rayidx) const {
     return direction[get_direction_index<true>(pointidx, rayidx)];
 }
 
 /// Get the direction of a ray
 /// @param[in] pointidx: Index of the point
 /// @param[in] rayidx: Index of the ray
-template <>
-Vector3D Rays ::get_direction<false>(const Size pointidx, const Size rayidx) const {
+template <> Vector3D Rays ::get_direction<false>(const Size pointidx, const Size rayidx) const {
     return direction[get_direction_index<false>(pointidx, rayidx)];
 }
 
 /// Get the antipodal direction of a ray
 /// @param[in] pointidx: Index of the point
 /// @param[in] rayidx: Index of the ray
-template <>
-Vector3D Rays ::get_antipod<true>(const Size pointidx, const Size rayidx) const {
+template <> Vector3D Rays ::get_antipod<true>(const Size pointidx, const Size rayidx) const {
     return direction[get_direction_index<true>(pointidx, get_antipod_index(rayidx))];
 }
 
 /// Get the antipodal direction of a ray
 /// @param[in] pointidx: Index of the point
 /// @param[in] rayidx: Index of the ray
-template <>
-Vector3D Rays ::get_antipod<false>(const Size pointidx, const Size rayidx) const {
+template <> Vector3D Rays ::get_antipod<false>(const Size pointidx, const Size rayidx) const {
     return direction[get_direction_index<false>(pointidx, get_antipod_index(rayidx))];
 }
 
@@ -138,15 +134,13 @@ Size Rays ::get_antipod_index(const Size rayidx) const { return antipod[rayidx];
 /// Get the weight of a ray
 /// @param[in] pointidx: Index of the point
 /// @param[in] rayidx: Index of the ray
-template <>
-Real Rays ::get_weight<true>(const Size pointidx, const Size rayidx) const {
+template <> Real Rays ::get_weight<true>(const Size pointidx, const Size rayidx) const {
     return weight[get_direction_index<true>(pointidx, rayidx)];
 }
 
 /// Get the weight of a ray
 /// @param[in] pointidx: Index of the point
 /// @param[in] rayidx: Index of the ray
-template <>
-Real Rays ::get_weight<false>(const Size pointidx, const Size rayidx) const {
+template <> Real Rays ::get_weight<false>(const Size pointidx, const Size rayidx) const {
     return weight[get_direction_index<false>(pointidx, rayidx)];
 }
diff --git a/src/model/model.cpp b/src/model/model.cpp
index f103a6df..e99ada18 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -315,12 +315,11 @@ int Model ::compute_radiation_field_shortchar_order_0() {
 
     Solver solver;
     const bool use_adaptive_directions = geometry.rays.use_adaptive_directions;
-    //TODO: create macro for this substitution instead of doing the if-else manually
+    // TODO: create macro for this substitution instead of doing the if-else manually
     if (use_adaptive_directions) {
         cout << "Using adaptive directions" << endl;
         solver.setup<CoMoving, true>(*this);
-    }
-    else {
+    } else {
         solver.setup<CoMoving, false>(*this);
     }
 
@@ -356,12 +355,11 @@ int Model ::compute_radiation_field_feautrier_order_2() {
 
     Solver solver;
     const bool use_adaptive_directions = geometry.rays.use_adaptive_directions;
-    //TODO: create macro for this substitution instead of doing the if-else manually
+    // TODO: create macro for this substitution instead of doing the if-else manually
     if (use_adaptive_directions) {
         cout << "Using adaptive directions" << endl;
         solver.setup<CoMoving, true>(*this);
-    }
-    else {
+    } else {
         solver.setup<CoMoving, false>(*this);
     }
 
@@ -402,8 +400,7 @@ int Model ::compute_radiation_field_feautrier_order_2_uv() {
     if (use_adaptive_directions) {
         cout << "Using adaptive directions" << endl;
         solver.setup<CoMoving, true>(*this);
-    }
-    else {
+    } else {
         solver.setup<CoMoving, false>(*this);
     }
 
@@ -443,8 +440,7 @@ int Model ::compute_radiation_field_feautrier_order_2_anis() {
     if (use_adaptive_directions) {
         cout << "Using adaptive directions" << endl;
         solver.setup<CoMoving, true>(*this);
-    }
-    else {
+    } else {
         solver.setup<CoMoving, false>(*this);
     }
 
@@ -483,8 +479,7 @@ int Model ::compute_radiation_field_feautrier_order_2_sparse() {
     if (use_adaptive_directions) {
         cout << "Using adaptive directions" << endl;
         solver.setup<CoMoving, true>(*this);
-    }
-    else {
+    } else {
         solver.setup<CoMoving, false>(*this);
     }
 
@@ -1090,7 +1085,8 @@ int Model ::compute_image_optical_depth_new(const Size ray_nr) {
 ///  Wrapper for the new imager
 ///////////////////////////////
 int Model ::compute_image_optical_depth_new(const Size ray_nr, const Size Nxpix, const Size Nypix) {
-    return compute_image_optical_depth_new(geometry.rays.get_direction<false>(0, ray_nr), Nxpix, Nypix);
+    return compute_image_optical_depth_new(
+        geometry.rays.get_direction<false>(0, ray_nr), Nxpix, Nypix);
 }
 
 ///  Wrapper for the new imager
@@ -1150,8 +1146,7 @@ int Model ::set_column() {
     if (use_adaptive_directions) {
         cout << "Using adaptive directions" << endl;
         solver.set_column<true>(*this);
-    }
-    else {
+    } else {
         solver.set_column<false>(*this);
     }
 
diff --git a/src/solver/solver.hpp b/src/solver/solver.hpp
index 1c373736..9510986a 100644
--- a/src/solver/solver.hpp
+++ b/src/solver/solver.hpp
@@ -80,7 +80,8 @@ struct Solver {
 
     template <Frame frame, bool use_adaptive_directions> inline void get_ray_lengths(Model& model);
 
-    template <Frame frame, bool use_adaptive_directions> inline Size get_ray_lengths_max(Model& model);
+    template <Frame frame, bool use_adaptive_directions>
+    inline Size get_ray_lengths_max(Model& model);
 
     // template <Frame frame>
     inline Size get_ray_lengths_max_new_imager(Model& model, Image& image, const Vector3D& ray_dir);
@@ -168,7 +169,7 @@ struct Solver {
 
     // Solvers only computing u
     ///////////////////////////
-    template <ApproximationType approx, bool use_adaptive_directions> 
+    template <ApproximationType approx, bool use_adaptive_directions>
     accel inline void solve_feautrier_order_2(Model& model);
 
     template <ApproximationType approx, bool use_adaptive_directions>
@@ -202,8 +203,7 @@ struct Solver {
 
     // Solvers for column densities
     ///////////////////////////////
-    template <bool use_adaptive_directions>
-    accel inline void set_column(Model& model) const;
+    template <bool use_adaptive_directions> accel inline void set_column(Model& model) const;
     template <bool use_adaptive_directions>
     accel inline Real get_column(const Model& model, const Size o, const Size r) const;
 };
diff --git a/src/solver/solver.tpp b/src/solver/solver.tpp
index 242e607d..990a383f 100644
--- a/src/solver/solver.tpp
+++ b/src/solver/solver.tpp
@@ -91,7 +91,7 @@ accel inline Real Solver ::get_dshift_max(const Model& model, const Size o) {
     return dshift_max;
 }
 
-template <Frame frame, bool use_adaptive_directions> 
+template <Frame frame, bool use_adaptive_directions>
 inline void Solver ::get_ray_lengths(Model& model) {
     for (Size rr = 0; rr < model.parameters->hnrays(); rr++) {
 
@@ -100,8 +100,9 @@ inline void Solver ::get_ray_lengths(Model& model) {
 
             const Real dshift_max = get_dshift_max(model, o);
 
-            model.geometry.lengths(rr, o) = model.geometry.get_ray_length<frame, use_adaptive_directions>(o, rr, dshift_max)
-                                          + model.geometry.get_ray_length<frame, use_adaptive_directions>(o, ar, dshift_max);
+            model.geometry.lengths(rr, o) =
+                model.geometry.get_ray_length<frame, use_adaptive_directions>(o, rr, dshift_max)
+                + model.geometry.get_ray_length<frame, use_adaptive_directions>(o, ar, dshift_max);
         })
 
         pc::accelerator::synchronize();
@@ -110,7 +111,8 @@ inline void Solver ::get_ray_lengths(Model& model) {
     model.geometry.lengths.copy_ptr_to_vec();
 }
 
-template <Frame frame, bool use_adaptive_directions> inline Size Solver ::get_ray_lengths_max(Model& model) {
+template <Frame frame, bool use_adaptive_directions>
+inline Size Solver ::get_ray_lengths_max(Model& model) {
     get_ray_lengths<frame, use_adaptive_directions>(model);
 
     Geometry& geo = model.geometry;
@@ -260,11 +262,12 @@ inline void Solver ::solve_feautrier_order_2_uv(Model& model) {
             nr_()[centre]    = o;
             shift_()[centre] = 1.0;
 
-            first_() =
-                trace_ray<CoMoving, use_adaptive_directions>(model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1)
-                + 1;
-            last_() =
-                trace_ray<CoMoving, use_adaptive_directions>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre) - 1;
+            first_() = trace_ray<CoMoving, use_adaptive_directions>(
+                           model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1)
+                     + 1;
+            last_() = trace_ray<CoMoving, use_adaptive_directions>(
+                          model.geometry, o, ar, dshift_max, +1, centre + 1, centre)
+                    - 1;
             n_tot_() = (last_() + 1) - first_();
 
             if (n_tot_() > 1) {
@@ -312,9 +315,11 @@ inline void Solver ::solve_feautrier_order_2_sparse(Model& model) {
 
         for (LineProducingSpecies& lspec : model.lines.lineProducingSpecies) {
             threaded_for(o, model.parameters->npoints(), {
-                const Vector3D nn = model.geometry.rays.get_direction<use_adaptive_directions>(o, rr);
-                const Size ar     = model.geometry.rays.get_antipod_index(rr);
-                const Real wt     = model.geometry.rays.get_weight<use_adaptive_directions>(o, rr) * two;
+                const Vector3D nn =
+                    model.geometry.rays.get_direction<use_adaptive_directions>(o, rr);
+                const Size ar = model.geometry.rays.get_antipod_index(rr);
+                const Real wt =
+                    model.geometry.rays.get_weight<use_adaptive_directions>(o, rr) * two;
 
                 const Real dshift_max = get_dshift_max(model, o);
 
@@ -324,9 +329,9 @@ inline void Solver ::solve_feautrier_order_2_sparse(Model& model) {
                 first_() = trace_ray<CoMoving, use_adaptive_directions>(
                                model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1)
                          + 1;
-                last_() =
-                    trace_ray<CoMoving, use_adaptive_directions>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre)
-                    - 1;
+                last_() = trace_ray<CoMoving, use_adaptive_directions>(
+                              model.geometry, o, ar, dshift_max, +1, centre + 1, centre)
+                        - 1;
                 n_tot_() = (last_() + 1) - first_();
 
                 if (n_tot_() > 1) {
@@ -337,7 +342,8 @@ inline void Solver ::solve_feautrier_order_2_sparse(Model& model) {
 
                             lspec.J(o, k) += lspec.quadrature.weights[z] * wt * Su_()[centre];
 
-                            update_Lambda<approx, use_adaptive_directions>(model, rr, lspec.nr_line[o][k][z]);
+                            update_Lambda<approx, use_adaptive_directions>(
+                                model, rr, lspec.nr_line[o][k][z]);
                         }
                     }
                 } else {
@@ -386,9 +392,10 @@ inline void Solver ::solve_feautrier_order_2_anis(Model& model) {
 
         for (LineProducingSpecies& lspec : model.lines.lineProducingSpecies) {
             threaded_for(o, model.parameters->npoints(), {
-                const Size ar     = model.geometry.rays.get_antipod_index(rr);
-                const Real wt     = model.geometry.rays.get_weight<use_adaptive_directions>(o, rr);
-                const Vector3D nn = model.geometry.rays.get_direction<use_adaptive_directions>(o, rr);
+                const Size ar = model.geometry.rays.get_antipod_index(rr);
+                const Real wt = model.geometry.rays.get_weight<use_adaptive_directions>(o, rr);
+                const Vector3D nn =
+                    model.geometry.rays.get_direction<use_adaptive_directions>(o, rr);
 
                 const Real wt_0    = inv_sqrt2 * (three * nn.z() * nn.z() - one);
                 const Real wt_1_Re = -sqrt3 * nn.x() * nn.z();
@@ -404,9 +411,9 @@ inline void Solver ::solve_feautrier_order_2_anis(Model& model) {
                 first_() = trace_ray<CoMoving, use_adaptive_directions>(
                                model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1)
                          + 1;
-                last_() =
-                    trace_ray<CoMoving, use_adaptive_directions>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre)
-                    - 1;
+                last_() = trace_ray<CoMoving, use_adaptive_directions>(
+                              model.geometry, o, ar, dshift_max, +1, centre + 1, centre)
+                        - 1;
                 n_tot_() = (last_() + 1) - first_();
 
                 if (n_tot_() > 1) {
@@ -481,9 +488,9 @@ inline void Solver ::solve_feautrier_order_2(Model& model) {
                 first_() = trace_ray<CoMoving, use_adaptive_directions>(
                                model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1)
                          + 1;
-                last_() =
-                    trace_ray<CoMoving, use_adaptive_directions>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre)
-                    - 1;
+                last_() = trace_ray<CoMoving, use_adaptive_directions>(
+                              model.geometry, o, ar, dshift_max, +1, centre + 1, centre)
+                        - 1;
                 n_tot_() = (last_() + 1) - first_();
 
                 if (n_tot_() > 1) {
@@ -492,7 +499,8 @@ inline void Solver ::solve_feautrier_order_2(Model& model) {
 
                         model.radiation.u(rr_loc, o, f) = Su_()[centre];
                         model.radiation.J(o, f) +=
-                            Su_()[centre] * two * model.geometry.rays.get_weight<use_adaptive_directions>(o, rr);
+                            Su_()[centre] * two
+                            * model.geometry.rays.get_weight<use_adaptive_directions>(o, rr);
 
                         update_Lambda<approx, use_adaptive_directions>(model, rr, f);
                     }
@@ -500,8 +508,9 @@ inline void Solver ::solve_feautrier_order_2(Model& model) {
                     for (Size f = 0; f < model.parameters->nfreqs(); f++) {
                         model.radiation.u(rr_loc, o, f) =
                             boundary_intensity(model, o, model.radiation.frequencies.nu(o, f));
-                        model.radiation.J(o, f) += two * model.geometry.rays.get_weight<use_adaptive_directions>(o, rr)
-                                                 * model.radiation.u(rr_loc, o, f);
+                        model.radiation.J(o, f) +=
+                            two * model.geometry.rays.get_weight<use_adaptive_directions>(o, rr)
+                            * model.radiation.u(rr_loc, o, f);
                     }
                 }
             })
@@ -538,8 +547,10 @@ inline void Solver ::image_feautrier_order_2(Model& model, const Size rr) {
         ;
 
         first_() =
-            trace_ray<Rest, false>(model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1) + 1;
-        last_()  = trace_ray<Rest, false>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre) - 1;
+            trace_ray<Rest, false>(model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1)
+            + 1;
+        last_() =
+            trace_ray<Rest, false>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre) - 1;
         n_tot_() = (last_() + 1) - first_();
 
         if (n_tot_() > 1) {
@@ -629,8 +640,9 @@ inline void Solver ::image_feautrier_order_2_for_point(Model& model, const Size
     nr_()[centre]    = o;
     shift_()[centre] = model.geometry.get_shift<Rest, false>(o, rr, o, 0.0);
 
-    first_() = trace_ray<Rest, false>(model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1) + 1;
-    last_()  = trace_ray<Rest, false>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre) - 1;
+    first_() =
+        trace_ray<Rest, false>(model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1) + 1;
+    last_() = trace_ray<Rest, false>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre) - 1;
     n_tot_() = (last_() + 1) - first_();
 
     model.S_ray.resize(n_tot_(), model.parameters->nfreqs());
@@ -658,8 +670,10 @@ inline void Solver ::image_optical_depth(Model& model, const Size rr) {
         ;
 
         first_() =
-            trace_ray<Rest, false>(model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1) + 1;
-        last_()  = trace_ray<Rest, false>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre) - 1;
+            trace_ray<Rest, false>(model.geometry, o, rr, dshift_max, -1, centre - 1, centre - 1)
+            + 1;
+        last_() =
+            trace_ray<Rest, false>(model.geometry, o, ar, dshift_max, +1, centre + 1, centre) - 1;
         n_tot_() = (last_() + 1) - first_();
 
         if (n_tot_() > 1) {
@@ -1465,7 +1479,8 @@ accel inline void Solver ::update_Lambda(Model& model, const Size rr, const Size
         Matrix<Real>& L_lower = L_lower_();
         // Vector<Real  >& inverse_chi = inverse_chi_();
 
-        const Real w_ang = two * model.geometry.rays.get_weight<use_adaptive_directions>(nr[centre], rr);
+        const Real w_ang =
+            two * model.geometry.rays.get_weight<use_adaptive_directions>(nr[centre], rr);
 
         const Size l = freqs.corresponding_l_for_spec[f]; // index of species
         const Size k = freqs.corresponding_k_for_tran[f]; // index of transition
@@ -2520,8 +2535,7 @@ accel inline void Solver ::set_boundary_condition(Model& model) const {
     }
 }
 
-template <bool use_adaptive_directions>
-inline void Solver ::set_column(Model& model) const {
+template <bool use_adaptive_directions> inline void Solver ::set_column(Model& model) const {
     model.column.resize(model.parameters->nrays(), model.parameters->npoints());
 
     for (Size rr = 0; rr < model.parameters->hnrays(); rr++) {