diff --git a/.github/workflows/ci-checks.yml b/.github/workflows/ci-checks.yml
index fa1544e23..d94a312ff 100644
--- a/.github/workflows/ci-checks.yml
+++ b/.github/workflows/ci-checks.yml
@@ -29,14 +29,19 @@ jobs:
                   cppcheck --force --error-exitcode=1 --language=c++ --std=c++20 --enable=all --inconclusive --suppress=knownConditionTrueFalse --suppress=constParameterCallback --suppress=redundantAssignment --suppress=knownArgument --suppress=missingIncludeSystem --suppress=unusedFunction --check-level=exhaustive .
 
     clang-format:
-        runs-on: ubuntu-24.04
+        runs-on: macos-14
         steps:
             - uses: actions/checkout@v4
 
+            - name: install clang-format
+              run: |
+                  brew update
+                  brew install clang-format || true
+
             - name: Run clang-format style check
               uses: jidicula/clang-format-action@v4.13.0
               with:
-                  clang-format-version: '18'
+                  clang-format-version: '19'
                   check-path: .
 
     clang-tidy:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index a082f1615..994f31fa1 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -22,7 +22,7 @@ repos:
             args: [--fix=lf]
           - id: trailing-whitespace
     - repo: https://github.com/pre-commit/mirrors-clang-format
-      rev: 'v18.1.8'
+      rev: 'v19.1.0'
       hooks:
           - id: clang-format
             files: '\.(c|cc|cpp|h|hpp|cxx|hh|inc)$'
diff --git a/decay.cc b/decay.cc
index 09064e9c1..93efd7f19 100644
--- a/decay.cc
+++ b/decay.cc
@@ -501,8 +501,8 @@ auto sample_decaytime(const int decaypathindex, const double tdecaymin, const do
 }
 
 constexpr auto calculate_decaychain(const double firstinitabund, const std::vector<double> &lambdas,
-                                    const int num_nuclides, const double timediff,
-                                    const bool useexpansionfactor) -> double {
+                                    const int num_nuclides, const double timediff, const bool useexpansionfactor)
+    -> double {
   // calculate final number abundance from multiple decays, e.g., Ni56 -> Co56 -> Fe56 (nuc[0] -> nuc[1] -> nuc[2])
   // the top nuclide initial abundance is set and the chain-end abundance is returned (all intermediates nuclides
   // are assumed to start with zero abundance)
@@ -615,8 +615,8 @@ auto get_nuc_massfrac(const int modelgridindex, const int z, const int a, const
   return nuctotal;
 }
 
-auto get_endecay_to_tinf_per_ejectamass_at_time(const int modelgridindex, const int decaypathindex,
-                                                const double time) -> double
+auto get_endecay_to_tinf_per_ejectamass_at_time(const int modelgridindex, const int decaypathindex, const double time)
+    -> double
 // returns decay energy [erg/g] that would be released from time tstart [s] to time infinity by a given decaypath
 {
   // e.g. Ni56 -> Co56, represents the decay of Co56 nuclei
@@ -653,8 +653,8 @@ auto get_endecay_to_tinf_per_ejectamass_at_time(const int modelgridindex, const
 }
 
 auto get_endecay_per_ejectamass_t0_to_time_withexpansion_chain_numerical(const int modelgridindex,
-                                                                         const int decaypathindex,
-                                                                         const double tstart) -> double
+                                                                         const int decaypathindex, const double tstart)
+    -> double
 // just here as as check on the analytic result from get_endecay_per_ejectamass_t0_to_time_withexpansion()
 // this version does an Euler integration
 {
@@ -1112,8 +1112,8 @@ void free_decaypath_energy_per_mass() {
 }
 
 // energy release rate in form of kinetic energy of positrons, electrons, and alpha particles in [erg/s/g]
-[[nodiscard]] auto get_particle_injection_rate(const int modelgridindex, const double t,
-                                               const int decaytype) -> double {
+[[nodiscard]] auto get_particle_injection_rate(const int modelgridindex, const double t, const int decaytype)
+    -> double {
   double dep_sum = 0.;
   const auto num_nuclides = get_num_nuclides();
   for (int nucindex = 0; nucindex < num_nuclides; nucindex++) {
diff --git a/grid.cc b/grid.cc
index 84b23a08a..3c5f88c78 100644
--- a/grid.cc
+++ b/grid.cc
@@ -264,7 +264,9 @@ void set_elem_untrackedstable_abund_from_total(const int mgi, const int element,
   }
 
   // if (globals::rank_in_node == 0)
-  { modelgrid[mgi].initmassfracuntrackedstable[element] = massfrac_untrackedstable; }
+  {
+    modelgrid[mgi].initmassfracuntrackedstable[element] = massfrac_untrackedstable;
+  }
 
   // (isofracsum + massfracstable) might not exactly match elemabundance if we had to boost it to reach isofracsum
   set_elem_abundance(mgi, element, isofracsum + massfrac_untrackedstable);
@@ -1598,11 +1600,12 @@ template <size_t S1>
   return -1.;
 }
 
-auto get_coordboundary_distances_cylindrical2d(
-    const std::array<double, 3> &pkt_pos, const std::array<double, 3> &pkt_dir,
-    const std::array<double, 3> &pktposgridcoord, const std::array<double, 3> &pktvelgridcoord, const int cellindex,
-    const double tstart,
-    const std::array<double, 3> &cellcoordmax) -> std::tuple<std::array<double, 3>, std::array<double, 3>> {
+auto get_coordboundary_distances_cylindrical2d(const std::array<double, 3> &pkt_pos,
+                                               const std::array<double, 3> &pkt_dir,
+                                               const std::array<double, 3> &pktposgridcoord,
+                                               const std::array<double, 3> &pktvelgridcoord, const int cellindex,
+                                               const double tstart, const std::array<double, 3> &cellcoordmax)
+    -> std::tuple<std::array<double, 3>, std::array<double, 3>> {
   // to get the cylindrical intersection, get the spherical intersection with Z components set to zero, and the
   // propagation speed set to the xy component of the 3-velocity
 
@@ -2421,9 +2424,10 @@ auto get_totmassradionuclide(const int z, const int a) -> double {
 }
 
 // compute distance to a cell boundary.
-[[nodiscard]] __host__ __device__ auto boundary_distance(
-    const std::array<double, 3> &dir, const std::array<double, 3> &pos, const double tstart, const int cellindex,
-    enum cell_boundary *pkt_last_cross) -> std::tuple<double, int> {
+[[nodiscard]] __host__ __device__ auto boundary_distance(const std::array<double, 3> &dir,
+                                                         const std::array<double, 3> &pos, const double tstart,
+                                                         const int cellindex, enum cell_boundary *pkt_last_cross)
+    -> std::tuple<double, int> {
   if constexpr (FORCE_SPHERICAL_ESCAPE_SURFACE) {
     if (get_cell_r_inner(cellindex) > globals::vmax * globals::tmin) {
       return {0., -99};
diff --git a/ltepop.cc b/ltepop.cc
index 8cac5b535..022cda193 100644
--- a/ltepop.cc
+++ b/ltepop.cc
@@ -514,8 +514,8 @@ auto calculate_levelpop(const int modelgridindex, const int element, const int i
 }
 
 // Calculate the population of a level from either LTE or NLTE information
-__host__ __device__ auto get_levelpop(const int modelgridindex, const int element, const int ion,
-                                      const int level) -> double {
+__host__ __device__ auto get_levelpop(const int modelgridindex, const int element, const int ion, const int level)
+    -> double {
   double nn = 0.;
   if (use_cellcache) {
     assert_testmodeonly(modelgridindex == globals::cellcache[cellcacheslotid].cellnumber);
diff --git a/ltepop.h b/ltepop.h
index 14a930846..1cbdd039c 100644
--- a/ltepop.h
+++ b/ltepop.h
@@ -7,13 +7,13 @@
 [[nodiscard]] auto calculate_levelpop(int modelgridindex, int element, int ion, int level) -> double;
 [[nodiscard]] auto calculate_levelpop_lte(int modelgridindex, int element, int ion, int level) -> double;
 [[nodiscard]] auto get_levelpop(int modelgridindex, int element, int ion, int level) -> double;
-[[nodiscard]] auto calculate_sahafact(int element, int ion, int level, int upperionlevel, double T,
-                                      double E_threshold) -> double;
+[[nodiscard]] auto calculate_sahafact(int element, int ion, int level, int upperionlevel, double T, double E_threshold)
+    -> double;
 [[nodiscard]] auto get_nnion(int modelgridindex, int element, int ion) -> double;
 void calculate_ion_balance_nne(int modelgridindex);
 void calculate_cellpartfuncts(int modelgridindex, int element);
-[[nodiscard]] auto calculate_ionfractions(int element, int modelgridindex, double nne,
-                                          bool use_phi_lte) -> std::vector<double>;
+[[nodiscard]] auto calculate_ionfractions(int element, int modelgridindex, double nne, bool use_phi_lte)
+    -> std::vector<double>;
 void set_groundlevelpops(int modelgridindex, int element, float nne, bool force_lte);
 
 #endif  // LTEPOP_H
diff --git a/macroatom.cc b/macroatom.cc
index 4693ce74c..8e9140a30 100644
--- a/macroatom.cc
+++ b/macroatom.cc
@@ -151,8 +151,8 @@ auto calculate_macroatom_transitionrates(const int modelgridindex, const int ele
   return processrates;
 }
 
-auto do_macroatom_internal_down_same(const int element, const int ion, const int level,
-                                     const CellCacheLevels &chlevel) -> int {
+auto do_macroatom_internal_down_same(const int element, const int ion, const int level, const CellCacheLevels &chlevel)
+    -> int {
   const int ndowntrans = get_ndowntrans(element, ion, level);
 
   // printout("[debug] do_ma:   internal downward jump within current ionstage\n");
@@ -904,8 +904,8 @@ auto col_deexcitation_ratecoeff(const float T_e, const float nne, const double e
 
 // multiply by lower level population to get a rate per second
 auto col_excitation_ratecoeff(const float T_e, const float nne, const int element, const int ion, const int lower,
-                              const int uptransindex, const double epsilon_trans,
-                              const double lowerstatweight) -> double {
+                              const int uptransindex, const double epsilon_trans, const double lowerstatweight)
+    -> double {
   const auto &uptr = get_uptranslist(element, ion, lower)[uptransindex];
   const double coll_strength = uptr.coll_str;
   const double eoverkt = epsilon_trans / (KB * T_e);
diff --git a/macroatom.h b/macroatom.h
index 3265012c3..3b1428a67 100644
--- a/macroatom.h
+++ b/macroatom.h
@@ -10,11 +10,11 @@ void macroatom_close_file();
 void do_macroatom(Packet &pkt, const MacroAtomState &pktmastate);
 
 [[nodiscard]] auto rad_deexcitation_ratecoeff(int modelgridindex, int element, int ion, int lower, double epsilon_trans,
-                                              float A_ul, double upperstatweight, double nnlevelupper,
-                                              double t_current) -> double;
+                                              float A_ul, double upperstatweight, double nnlevelupper, double t_current)
+    -> double;
 [[nodiscard]] auto rad_excitation_ratecoeff(int modelgridindex, int element, int ion, int lower, int uptransindex,
-                                            double epsilon_trans, double nnlevel_lower, int lineindex,
-                                            double t_current) -> double;
+                                            double epsilon_trans, double nnlevel_lower, int lineindex, double t_current)
+    -> double;
 [[nodiscard]] auto rad_recombination_ratecoeff(float T_e, float nne, int element, int upperion, int upperionlevel,
                                                int lowerionlevel, int modelgridindex) -> double;
 [[nodiscard]] auto stim_recombination_ratecoeff(float nne, int element, int upperion, int upper, int lower,
diff --git a/nltepop.cc b/nltepop.cc
index 8bf3f9048..381cd0edb 100644
--- a/nltepop.cc
+++ b/nltepop.cc
@@ -120,8 +120,8 @@ void filter_nlte_matrix(const int element, gsl_matrix *rate_matrix, gsl_vector *
 }
 
 [[nodiscard]] auto get_total_rate(const int index_selected, const gsl_matrix *rate_matrix, const gsl_vector *popvec,
-                                  const bool into_level, const bool only_levels_below,
-                                  const bool only_levels_above) -> double {
+                                  const bool into_level, const bool only_levels_below, const bool only_levels_above)
+    -> double {
   double total_rate = 0.;
   assert_always(!only_levels_below || !only_levels_above);
 
diff --git a/nonthermal.cc b/nonthermal.cc
index 4a2ccbc1f..f95db1d01 100644
--- a/nonthermal.cc
+++ b/nonthermal.cc
@@ -694,8 +694,8 @@ auto get_xs_ionization_vector(std::array<double, SFPTS> &xs_vec, const collionro
 
 // distribution of secondary electron energies for primary electron with energy e_p
 // Opal, Peterson, & Beaty (1971)
-[[nodiscard]] constexpr auto Psecondary(const double e_p, const double epsilon, const double I,
-                                        const double J) -> double {
+[[nodiscard]] constexpr auto Psecondary(const double e_p, const double epsilon, const double I, const double J)
+    -> double {
   const double e_s = epsilon - I;
 
   if (e_p <= I || e_s < 0.) {
@@ -1081,8 +1081,8 @@ auto get_oneoverw(const int element, const int ion, const int modelgridindex) ->
 
 // the fraction of deposited energy that goes into ionising electrons in a particular shell
 auto calculate_nt_frac_ionization_shell(const int modelgridindex, const int element, const int ion,
-                                        const collionrow &collionrow,
-                                        const std::array<double, SFPTS> &yfunc) -> double {
+                                        const collionrow &collionrow, const std::array<double, SFPTS> &yfunc)
+    -> double {
   const double nnion = get_nnion(modelgridindex, element, ion);
   const double ionpot_ev = collionrow.ionpot_ev;
 
@@ -1105,8 +1105,8 @@ auto nt_ionization_ratecoeff_wfapprox(const int modelgridindex, const int elemen
 }
 
 auto calculate_nt_ionization_ratecoeff(const int modelgridindex, const int element, const int ion,
-                                       const bool assumeshellpotentialisvalence,
-                                       const std::array<double, SFPTS> &yfunc) -> double
+                                       const bool assumeshellpotentialisvalence, const std::array<double, SFPTS> &yfunc)
+    -> double
 // Integrate the ionization cross section over the electron degradation function to get the ionization rate coefficient
 // i.e. multiply this by ion population to get a rate of ionizations per second
 // Do not call during packet propagation, as the y vector may not be in memory!
@@ -1286,8 +1286,8 @@ auto nt_ionization_ratecoeff_sf(const int modelgridindex, const int element, con
 // returns the index of the first valid cross section point (en >= epsilon_trans)
 // all elements below this index are invalid and should not be used
 auto get_xs_excitation_vector(const int element, const int ion, const int lower, const int uptransindex,
-                              const double statweight_lower,
-                              const double epsilon_trans) -> std::tuple<std::array<double, SFPTS>, int> {
+                              const double statweight_lower, const double epsilon_trans)
+    -> std::tuple<std::array<double, SFPTS>, int> {
   std::array<double, SFPTS> xs_excitation_vec{};
   const auto &uptr = get_uptranslist(element, ion, lower)[uptransindex];
   if (uptr.coll_str >= 0) {
@@ -1344,8 +1344,8 @@ auto get_xs_excitation_vector(const int element, const int ion, const int lower,
 // returns the rate coefficient in s^-1 divided by deposition rate density in erg/cm^3/s
 auto calculate_nt_excitation_ratecoeff_perdeposition(const std::array<double, SFPTS> &yvec, const int element,
                                                      const int ion, const int lower, const int uptransindex,
-                                                     const double statweight_lower,
-                                                     const double epsilon_trans) -> double {
+                                                     const double statweight_lower, const double epsilon_trans)
+    -> double {
   const auto [xs_excitation_vec, xsstartindex] =
       get_xs_excitation_vector(element, ion, lower, uptransindex, statweight_lower, epsilon_trans);
 
@@ -2264,8 +2264,8 @@ __host__ __device__ auto nt_ionization_ratecoeff(const int modelgridindex, const
 }
 
 __host__ __device__ auto nt_excitation_ratecoeff(const int modelgridindex, const int element, const int ion,
-                                                 const int lowerlevel, const int uptransindex,
-                                                 const int lineindex) -> double {
+                                                 const int lowerlevel, const int uptransindex, const int lineindex)
+    -> double {
   if constexpr (!NT_EXCITATION_ON) {
     return 0.;
   }
diff --git a/radfield.cc b/radfield.cc
index a7a3a2b41..c50324a0d 100644
--- a/radfield.cc
+++ b/radfield.cc
@@ -905,8 +905,8 @@ __host__ __device__ auto radfield(const double nu, const int modelgridindex) ->
 
 // return the integral of nu^3 / (exp(h nu / k T) - 1) from nu_lower to nu_upper
 // or if times_nu is true, the integral of nu^4 / (exp(h nu / k T) - 1) from nu_lower to nu_upper
-auto planck_integral_analytic(const double T_R, const double nu_lower, const double nu_upper,
-                              const bool times_nu) -> double {
+auto planck_integral_analytic(const double T_R, const double nu_lower, const double nu_upper, const bool times_nu)
+    -> double {
   double integral = 0.;
 
   if (times_nu) {
diff --git a/radfield.h b/radfield.h
index 8a335c1ae..406c785ca 100644
--- a/radfield.h
+++ b/radfield.h
@@ -34,8 +34,8 @@ void do_MPI_Bcast(int modelgridindex, int root, int root_node_id);
 void write_restart_data(FILE *gridsave_file);
 void read_restart_data(FILE *gridsave_file);
 void normalise_bf_estimators(int nts, int nts_prev, int titer, double deltat);
-[[nodiscard]] auto get_bfrate_estimator(int element, int lowerion, int lower, int phixstargetindex,
-                                        int modelgridindex) -> double;
+[[nodiscard]] auto get_bfrate_estimator(int element, int lowerion, int lower, int phixstargetindex, int modelgridindex)
+    -> double;
 void print_bfrate_contributions(int element, int lowerion, int lower, int phixstargetindex, int modelgridindex,
                                 double nnlowerlevel, double nnlowerion);
 void reset_bfrate_contributions(int modelgridindex);
diff --git a/ratecoeff.cc b/ratecoeff.cc
index 942c43942..f6715a1cd 100644
--- a/ratecoeff.cc
+++ b/ratecoeff.cc
@@ -996,8 +996,8 @@ __host__ __device__ auto get_spontrecombcoeff(int element, const int ion, const
 // multiply by upper ion population (or ground population if per_groundmultipletpop is true) and nne to get a rate
 auto calculate_ionrecombcoeff(const int modelgridindex, const float T_e, const int element, const int upperion,
                               const bool assume_lte, const bool collisional_not_radiative, const bool printdebug,
-                              const bool lower_superlevel_only, const bool per_groundmultipletpop,
-                              const bool stimonly) -> double {
+                              const bool lower_superlevel_only, const bool per_groundmultipletpop, const bool stimonly)
+    -> double {
   const int lowerion = upperion - 1;
   if (lowerion < 0) {
     return 0.;
diff --git a/ratecoeff.h b/ratecoeff.h
index 05aa066bb..c351b6979 100644
--- a/ratecoeff.h
+++ b/ratecoeff.h
@@ -17,19 +17,19 @@ void setup_photoion_luts();
 
 [[nodiscard]] auto select_continuum_nu(int element, int lowerion, int lower, int upperionlevel, float T_e) -> double;
 
-[[nodiscard]] auto interpolate_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex,
-                                                 double T) -> double;
+[[nodiscard]] auto interpolate_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex, double T)
+    -> double;
 
 [[nodiscard]] auto get_spontrecombcoeff(int element, int ion, int level, int phixstargetindex, float T_e) -> double;
-[[nodiscard]] auto get_stimrecombcoeff(int element, int lowerion, int level, int phixstargetindex,
-                                       int modelgridindex) -> double;
+[[nodiscard]] auto get_stimrecombcoeff(int element, int lowerion, int level, int phixstargetindex, int modelgridindex)
+    -> double;
 
 [[nodiscard]] auto get_bfcoolingcoeff(int element, int ion, int level, int phixstargetindex, float T_e) -> double;
 
-[[nodiscard]] auto get_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex,
-                                         int modelgridindex) -> double;
-[[nodiscard]] auto get_corrphotoioncoeff_ana(int element, int ion, int level, int phixstargetindex,
-                                             int modelgridindex) -> double;
+[[nodiscard]] auto get_corrphotoioncoeff(int element, int ion, int level, int phixstargetindex, int modelgridindex)
+    -> double;
+[[nodiscard]] auto get_corrphotoioncoeff_ana(int element, int ion, int level, int phixstargetindex, int modelgridindex)
+    -> double;
 
 [[nodiscard]] auto iongamma_is_zero(int nonemptymgi, int element, int ion) -> bool;
 
diff --git a/rpkt.h b/rpkt.h
index 077bdecf4..871bf59ff 100644
--- a/rpkt.h
+++ b/rpkt.h
@@ -102,8 +102,8 @@ constexpr auto closest_transition(const double nu_cmf, const int next_trans, con
   return matchindex;
 }
 
-[[nodiscard]] inline auto get_ionestimindex_nonemptymgi(const int nonemptymgi, const int element,
-                                                        const int ion) -> int {
+[[nodiscard]] inline auto get_ionestimindex_nonemptymgi(const int nonemptymgi, const int element, const int ion)
+    -> int {
   assert_testmodeonly(ion >= 0);
   assert_testmodeonly(ion < get_nions(element) - 1);
   const int groundcontindex = globals::elements[element].ions[ion].groundcontindex;
@@ -111,8 +111,8 @@ constexpr auto closest_transition(const double nu_cmf, const int next_trans, con
   return (nonemptymgi * globals::nbfcontinua_ground) + groundcontindex;
 }
 
-inline auto keep_this_cont(int element, const int ion, const int level, const int modelgridindex,
-                           const float nnetot) -> bool {
+inline auto keep_this_cont(int element, const int ion, const int level, const int modelgridindex, const float nnetot)
+    -> bool {
   if constexpr (DETAILED_BF_ESTIMATORS_ON) {
     return grid::get_elem_abundance(modelgridindex, element) > 0;
   }
diff --git a/sn3d.h b/sn3d.h
index 458e501c9..3fbe4a42a 100644
--- a/sn3d.h
+++ b/sn3d.h
@@ -313,8 +313,8 @@ inline void check_already_running() {
 #endif
 }
 
-constexpr auto get_range_chunk(const ptrdiff_t size, const ptrdiff_t nchunks,
-                               const ptrdiff_t nchunk) -> std::tuple<ptrdiff_t, ptrdiff_t> {
+constexpr auto get_range_chunk(const ptrdiff_t size, const ptrdiff_t nchunks, const ptrdiff_t nchunk)
+    -> std::tuple<ptrdiff_t, ptrdiff_t> {
   assert_always(size >= 0);
   assert_always(nchunks >= 0);
   assert_always(nchunk >= 0);
diff --git a/stats.cc b/stats.cc
index 198db4e7b..4784986aa 100644
--- a/stats.cc
+++ b/stats.cc
@@ -111,8 +111,8 @@ void increment_ion_stats_contabsorption(const Packet &pkt, const int modelgridin
   }
 }
 
-auto get_ion_stats(const int modelgridindex, const int element, const int ion,
-                   enum ionstattypes ionstattype) -> double {
+auto get_ion_stats(const int modelgridindex, const int element, const int ion, enum ionstattypes ionstattype)
+    -> double {
   assert_always(ion < get_nions(element));
   assert_always(ionstattype < ION_STAT_COUNT);
   const int uniqueionindex = get_uniqueionindex(element, ion);
diff --git a/thermalbalance.h b/thermalbalance.h
index e99b3b806..91a44c73d 100644
--- a/thermalbalance.h
+++ b/thermalbalance.h
@@ -26,8 +26,8 @@ struct HeatingCoolingRates {
 
 void call_T_e_finder(int modelgridindex, int timestep, double t_current, double T_min, double T_max,
                      HeatingCoolingRates *heatingcoolingrates, const std::vector<double> &bfheatingcoeffs);
-[[nodiscard]] auto get_bfheatingcoeff_ana(int element, int ion, int level, int phixstargetindex, double T_R,
-                                          double W) -> double;
+[[nodiscard]] auto get_bfheatingcoeff_ana(int element, int ion, int level, int phixstargetindex, double T_R, double W)
+    -> double;
 void calculate_bfheatingcoeffs(int modelgridindex, std::vector<double> &bfheatingcoeffs);
 
 #endif  // THERMALBALANCE_H
diff --git a/update_grid.cc b/update_grid.cc
index 3d375f2ce..502eca330 100644
--- a/update_grid.cc
+++ b/update_grid.cc
@@ -1155,7 +1155,9 @@ void update_grid(FILE *estimators_file, const int nts, const int nts_prev, const
 #ifdef _OPENMP
 #pragma omp critical(estimators_file)
 #endif
-      { write_to_estimators_file(estimators_file, mgi, nts, titer, &heatingcoolingrates); }
+      {
+        write_to_estimators_file(estimators_file, mgi, nts, titer, &heatingcoolingrates);
+      }
     }  // end parallel for loop over all modelgrid cells
 
   }  // end OpenMP parallel section
diff --git a/vectors.h b/vectors.h
index 8bee8c776..68616eca1 100644
--- a/vectors.h
+++ b/vectors.h
@@ -55,8 +55,8 @@ template <size_t S1, size_t S2>
 //   dir1: direction unit vector in frame1
 //   vel: velocity of frame2 relative to frame1
 //   dir2: direction vector in frame2
-[[nodiscard]] constexpr auto angle_ab(const std::array<double, 3> &dir1,
-                                      const std::array<double, 3> &vel) -> std::array<double, 3> {
+[[nodiscard]] constexpr auto angle_ab(const std::array<double, 3> &dir1, const std::array<double, 3> &vel)
+    -> std::array<double, 3> {
   const double vsqr = dot(vel, vel) / CLIGHTSQUARED;
   const double gamma_rel = 1. / std::sqrt(1 - vsqr);
 
@@ -77,8 +77,8 @@ template <size_t S1, size_t S2>
 //   prop_time: the propagation time of the packet
 // returns: the ratio f = (nu_cmf / nu_rf) ^ 2
 [[nodiscard]] constexpr auto doppler_squared_nucmf_on_nurf(const std::array<double, 3> &pos_rf,
-                                                           const std::array<double, 3> &dir_rf,
-                                                           const double prop_time) -> double {
+                                                           const std::array<double, 3> &dir_rf, const double prop_time)
+    -> double {
   // velocity of the comoving frame relative to the rest frame
   const auto vel_rf = get_velocity(pos_rf, prop_time);
 
@@ -212,8 +212,8 @@ constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double {
 
 // Rotation angle from the scattering plane
 [[nodiscard]] constexpr auto get_rot_angle(const std::array<double, 3> &n1, const std::array<double, 3> &n2,
-                                           const std::array<double, 3> &ref1,
-                                           const std::array<double, 3> &ref2) -> double {
+                                           const std::array<double, 3> &ref1, const std::array<double, 3> &ref2)
+    -> double {
   // We need to rotate Stokes Parameters to (or from) the scattering plane from (or to)
   // the meridian frame such that Q=1 is in the scattering plane and along ref1
 
@@ -292,8 +292,8 @@ constexpr auto move_pkt_withtime(Packet &pkt, const double distance) -> double {
 }
 
 // Routine to transform the Stokes Parameters from RF to CMF
-constexpr auto frame_transform(const std::array<double, 3> &n_rf, double *Q, double *U,
-                               const std::array<double, 3> &v) -> std::array<double, 3> {
+constexpr auto frame_transform(const std::array<double, 3> &n_rf, double *Q, double *U, const std::array<double, 3> &v)
+    -> std::array<double, 3> {
   // Meridian frame in the RF
   const auto [ref1_rf, ref2_rf] = meridian(n_rf);