Merge branch 'adaptive_ray_directions' of github.com:Magritte-code/Ma…

…gritte into adaptive_ray_directions

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.");

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)];
 }

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);
     }
 

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;
 };