diff --git a/.github/workflows/test_nightly.yml b/.github/workflows/test_nightly.yml
index f32773617..87bf7d711 100644
--- a/.github/workflows/test_nightly.yml
+++ b/.github/workflows/test_nightly.yml
@@ -93,7 +93,7 @@ jobs:
       run: |
         source venv/bin/activate
         cd main_class/python
-        TEST_LEVEL=3 CLASS_VERBOSE=1 python "$(which nosetests)" -a dump_ini_files test_class.py
+        TEST_LEVEL=3 VALGRIND_MODE=1 CLASS_VERBOSE=1 python "$(which nosetests)" -a dump_ini_files test_class.py
         deactivate
     - name: Valgrind 🤖
       run: |
diff --git a/.github/workflows/test_on_pull_request.yml b/.github/workflows/test_on_pull_request.yml
index ad0e6cb5d..f6cc3145b 100644
--- a/.github/workflows/test_on_pull_request.yml
+++ b/.github/workflows/test_on_pull_request.yml
@@ -38,7 +38,7 @@ jobs:
       run: |
         source venv/bin/activate
         cd main_class/python
-        OMP_NUM_THREADS=16 COMPARE_OUTPUT_REF=1 TEST_LEVEL=2 nosetests -v -a test_scenario test_class.py --nologcapture --nocapture
+        OMP_NUM_THREADS=16 COMPARE_OUTPUT_REF=1 TEST_LEVEL=3 nosetests -v -a test_scenario test_class.py --nologcapture --nocapture
         deactivate
     - name: Upload plots 📤
       if: success() || failure()
diff --git a/explanatory.ini b/explanatory.ini
index b483cd901..c74f297bd 100644
--- a/explanatory.ini
+++ b/explanatory.ini
@@ -151,6 +151,9 @@ gauge = synchronous
 #      (default: set to 'no')
 #nbody_gauge_transfer_functions = yes
 
+# 4.c) Do you want the source functions for total non-relativistic matter, delta_m and theta_m, and baryon+cdm, delta_cb and theta_cb, to be outputed in the current gauge (the one selected in 4.a), instead of being automatically expressed as a gauge-invariant (GI) quantity, likeL: delta_m^GI=delta_m + 3 a H theta_m/k2, theta_m^GI=theta_m+alpha*k2 (default: no, that is, convert to GI)
+#matter_source_in_current_gauge= no
+
 # 5) Hubble parameter : either 'H0' in km/s/Mpc or 'h' (or 'theta_s_100'), where
 #    the latter is the peak scale parameter defined exactly as 100(ds_dec/da_dec)
 #    with a decoupling time given by maximum of visibility function (quite different
diff --git a/external/RealSpaceInterface/Calc2D/CalculationClass.py b/external/RealSpaceInterface/Calc2D/CalculationClass.py
index 0a6b5bb40..c63af4204 100644
--- a/external/RealSpaceInterface/Calc2D/CalculationClass.py
+++ b/external/RealSpaceInterface/Calc2D/CalculationClass.py
@@ -71,7 +71,7 @@ def getData(self, redshiftindex):
         Valuenew = dict()
         FValue_abs = np.abs(self.FValue)
         _min, _max = FValue_abs.min(), FValue_abs.max()
-        dimensions = (self.endshape[0] / 2, self.endshape[1])
+        dimensions = (self.endshape[0] // 2, self.endshape[1])
         for quantity, FT in FValuenew.items():
             FT_abs = np.abs(FT)
             FT_normalized = cv2.resize(FT_abs, dimensions).ravel()
@@ -98,7 +98,7 @@ def getInitialData(self):
         return Value.ravel(), cv2.resize(
             (np.abs(self.FValue) - np.abs(self.FValue).min()) /
             (np.abs(self.FValue).max() - np.abs(self.FValue).min()),
-            (self.endshape[0] / 2, self.endshape[1])).ravel(), (minimum,
+            (self.endshape[0] // 2, self.endshape[1])).ravel(), (minimum,
                                                                 maximum)
 
     def getTransferData(self, redshiftindex):
diff --git a/external/RealSpaceInterface/Calc2D/Database.py b/external/RealSpaceInterface/Calc2D/Database.py
index bc8f006c2..af5e4176c 100644
--- a/external/RealSpaceInterface/Calc2D/Database.py
+++ b/external/RealSpaceInterface/Calc2D/Database.py
@@ -14,22 +14,22 @@ def __init__(self, directory, db_file="database.dat"):
         self.db_path = os.path.join(directory, db_file)
         if not os.path.exists(self.db_path):
             logging.info("No database found; Creating one at {}.".format(self.db_path))
-            with open(self.db_path, "w") as f:
+            with open(self.db_path, "wb") as f:
                 pickle.dump(dict(), f)
 
         self.db = self.__read_database()
 
     def __read_database(self):
-        with open(self.db_path) as f:
-            return pickle.load(f)
+        with open(self.db_path, "rb") as f:
+            return pickle.load(f, encoding='latin1')
 
     def __write_database(self):
-        with open(self.db_path, "w") as f:
-            pickle.dump(self.db, f)
+        with open(self.db_path, "wb") as f:
+            pickle.dump(self.db, f, protocol=pickle.HIGHEST_PROTOCOL)
 
     def __create_file(self, data):
         filename = str(uuid.uuid4())
-        with open(os.path.join(self.directory, filename), "w") as f:
+        with open(os.path.join(self.directory, filename), "wb") as f:
             pickle.dump(data, f)
         return filename
 
@@ -40,8 +40,8 @@ def __getitem__(self, key):
         frozen_key = self.__get_frozen_key(key)
         if frozen_key in self.db:
             filename = self.db[frozen_key]
-            with open(os.path.join(self.directory, filename)) as f:
-                return pickle.load(f)
+            with open(os.path.join(self.directory, filename), "rb") as f:
+                return pickle.load(f, encoding='latin1')
         else:
             raise KeyError("No data for key: {}".format(key))
 
@@ -55,4 +55,4 @@ def __contains__(self, key):
         Return whether `self` contains a record
         for the given `key`.
         """
-        return self.__get_frozen_key(key) in self.db
\ No newline at end of file
+        return self.__get_frozen_key(key) in self.db
diff --git a/external/RealSpaceInterface/Calc2D/TransferFunction.py b/external/RealSpaceInterface/Calc2D/TransferFunction.py
index 8747f007f..135bc195f 100644
--- a/external/RealSpaceInterface/Calc2D/TransferFunction.py
+++ b/external/RealSpaceInterface/Calc2D/TransferFunction.py
@@ -42,6 +42,7 @@ def ComputeTransferFunctionList(cosmologicalParameters, redshift, kperdecade=200
     class_settings.update({
         "output": "mTk",
         "gauge": "newtonian",
+        "matter_source_in_current_gauge": "yes",
         "evolver": "1",
         "P_k_max_h/Mpc": P_k_max,
         "k_per_decade_for_pk": kperdecade,
diff --git a/external/RealSpaceInterface/README b/external/RealSpaceInterface/README
index 4fc84d3e3..83b3ca585 100644
--- a/external/RealSpaceInterface/README
+++ b/external/RealSpaceInterface/README
@@ -1,8 +1,14 @@
+CLASS real space interface
+07.2015: started by Max Beutelspacher
+09.2018: improved by Georgios Samaras and released in class v2.7.0
+10.2025: for v3.2.5: made compatible with python 3, credits Timothy Morton and Mowen Zhao
+
+
 For installation of python packages, run
 
     pip install -r requirements.txt
 
-Launch the application with
+To see the real space interface, go to this directory (external/RealSpaceInterface/) and launch the application with
 
     python tornadoserver.py
 
@@ -15,4 +21,3 @@ Then in any browser open the URL
 Cache files are located in cache/, so to clear the cache, run
 
     rm cache/*
-
diff --git a/external/RealSpaceInterface/tornadoserver.py b/external/RealSpaceInterface/tornadoserver.py
index 808119e25..299369ba4 100644
--- a/external/RealSpaceInterface/tornadoserver.py
+++ b/external/RealSpaceInterface/tornadoserver.py
@@ -143,8 +143,8 @@ def send_frame(self, redindex):
         self.write_message(json.dumps({'type': 'extrema', 'extrema': extrema}))
         progress = float(redindex) / len(self.calc.redshift)
 
-        real = {quantity: base64.b64encode(data.astype(np.float32)) for quantity, data in Valuenew.iteritems()}
-        transfer = {quantity: base64.b64encode(data.astype(np.float32)) for quantity, data in TransferData.iteritems()}
+        real = {quantity: base64.b64encode(data.astype(np.float32)).decode() for quantity, data in Valuenew.items()}
+        transfer = {quantity: base64.b64encode(data.astype(np.float32)).decode() for quantity, data in TransferData.items()}
         self.write_message(
             json.dumps({
                 'type': 'data',
@@ -166,9 +166,9 @@ def send_initial_state(self):
         extremastring = json.dumps({'type': 'extrema', 'extrema': extrema})
         datastring = json.dumps({
             'type': 'data',
-            'real': base64.b64encode(Value.astype(np.float32)),
+            'real': base64.b64encode(Value.astype(np.float32)).decode(),
             'fourier': [],
-            'transfer': base64.b64encode(TransferData.astype(np.float32)),
+            'transfer': base64.b64encode(TransferData.astype(np.float32)).decode(),
             'k': krange.tolist()
             })
         self.write_message(extremastring)
@@ -212,6 +212,7 @@ def set_cosmological_parameters(self, cosmologicalParameters):
 
             z_of_decoupling = self.calc.z_dec
             frame_of_decoupling = np.argmin(np.abs(z_of_decoupling - self.calc.redshift))
+            frame_of_decoupling = int(frame_of_decoupling) # this is np.int64 object, which json doesn't like
             if self.calc.redshift[frame_of_decoupling] > z_of_decoupling:
                 frame_of_decoupling -= 1
             messages.append({
diff --git a/include/common.h b/include/common.h
index 2100f1cb6..699fcd18c 100644
--- a/include/common.h
+++ b/include/common.h
@@ -15,7 +15,7 @@
 #ifndef __COMMON__
 #define __COMMON__
 
-#define _VERSION_ "v3.2.4"
+#define _VERSION_ "v3.2.5"
 
 /* @cond INCLUDE_WITH_DOXYGEN */
 
diff --git a/include/fourier.h b/include/fourier.h
index af8356d4d..e63136cfc 100644
--- a/include/fourier.h
+++ b/include/fourier.h
@@ -356,6 +356,7 @@ extern "C" {
 
   int fourier_get_k_list(
                          struct precision *ppr,
+                         struct primordial *ppm,
                          struct perturbations * ppt,
                          struct fourier * pfo
                          );
diff --git a/include/perturbations.h b/include/perturbations.h
index bf4c4e976..be28af61e 100644
--- a/include/perturbations.h
+++ b/include/perturbations.h
@@ -198,6 +198,8 @@ struct perturbations
 
   enum possible_gauges gauge; /**< gauge in which to perform this calculation */
 
+  short has_matter_source_in_current_gauge; /**< whether to keep matter and baryon+CDM sources in current gauge, instead of automatic conversion to gauge-invariant variables */
+
   //@}
 
   /** @name - indices running on modes (scalar, vector, tensor) */
diff --git a/python/test_class.py b/python/test_class.py
index 4e0b6baf6..a28b0f066 100644
--- a/python/test_class.py
+++ b/python/test_class.py
@@ -93,6 +93,8 @@
 COMPARE_OUTPUT_GAUGE = bool(int(os.getenv('COMPARE_OUTPUT_GAUGE', '0'))) # Compare synchronous and Newtonian gauge outputs?
 COMPARE_OUTPUT_REF = bool(int(os.getenv('COMPARE_OUTPUT_REF', '0'))) # Compare classy with classyref?
 POWER_ALL = bool(int(os.getenv('POWER_ALL', '0'))) # Combine every extension with each other? (Very slow!)
+VALGRIND_MODE = bool(int(os.getenv('VALGRIND_MODE', '0')))
+
 TEST_LEVEL = int(os.getenv('TEST_LEVEL', '0')) # 0 <= TEST_LEVEL <= 3
 
 if COMPARE_OUTPUT_REF:
@@ -108,6 +110,7 @@
 COMPARE_PK_RELATIVE_ERROR = 1e-2
 COMPARE_PK_RELATIVE_ERROR_GAUGE = 5*1e-2
 COMPARE_PK_ABSOLUTE_ERROR = 1e-20
+COMPARE_SHOOTING_RELATIVE_ERROR = 1e-5
 
 # Dictionary of models to test the wrapper against. Each of these scenario will
 # be run against all the possible output choices (nothing, tCl, mPk, etc...),
@@ -116,6 +119,10 @@
 # against. Indeed, when not specifying a field, CLASS takes the default input.
 CLASS_INPUT = {}
 
+# power = compare ALL settings with each of the entries
+# normal = compare only the power settings with this entry
+# onlyfull = compare only for 'full' output settings (see below)
+
 CLASS_INPUT['Output_spectra'] = (
     [{'output': 'mPk', 'P_k_max_1/Mpc': 2},
      {'output': 'tCl'},
@@ -126,7 +133,8 @@
     'power')
 
 CLASS_INPUT['Nonlinear'] = (
-    [{'non linear': 'halofit'}],
+    [{'non_linear': 'halofit'},
+    {'non_linear': 'hmcode'}],
     'power')
 
 CLASS_INPUT['Lensing'] = (
@@ -148,44 +156,65 @@
     CLASS_INPUT['modes'] = (
         [{'modes': 't'},
         {'modes': 's, t'}],
-        'power')
+        'normal')
 
     CLASS_INPUT['Tensor_method'] = (
         [{'tensor method': 'exact'},
         {'tensor method': 'photons'}],
-        'power')
+        'onlyfull')
 
 if TEST_LEVEL > 2:
     CLASS_INPUT['Isocurvature_modes'] = (
         [{'ic': 'ad,nid,cdi', 'c_ad_cdi': -0.5}],
-        'normal')
+        'onlyfull')
 
     CLASS_INPUT['Scalar_field'] = (
         [{'Omega_scf': 0.1, 'attractor_ic_scf': 'yes',
         'scf_parameters': '10, 0, 0, 0'}],
-        'normal')
+        'onlyfull')
 
     CLASS_INPUT['Inflation'] = (
-        [{'P_k_ini type': 'inflation_V'},
-        {'P_k_ini type': 'inflation_H'}],
-        'normal')
-    # CLASS_INPUT['Inflation'] = (
-    #     [{'P_k_ini type': 'inflation_V'},
-    #     {'P_k_ini type': 'inflation_H'},
-    #     {'P_k_ini type': 'inflation_V_end'}],
-    #     'normal')
+        [{'Pk_ini_type': 'inflation_V'},
+        {'Pk_ini_type': 'inflation_H'}],
+        'onlyfull')
+
+    CLASS_INPUT['Shooting'] = (
+        [{'sigma8': 0.8},
+        {'S8': 0.8},
+        {'theta_s_100': 1.04},
+        {'Neff': 4.0},
+        {'Neff': 2.0},
+        {'Neff': 2.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 3.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06, 'theta_s_100':1.04,'sigma8':0.8},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06, 'theta_s_100':1.04,'S8':0.8},
+        {'tau_reio':0.05}
+        ],
+        'onlyfull')
+# TODO :: possibly another one with nonlinear=hmcode in the future, but for now this is more stable
+FULL_SETTINGS = {'output':'mPk mTk vTk tCl pCl lCl sCl nCl','P_k_max_1/Mpc': 2,'lensing':'yes','modes':'s,t'}
 
 if POWER_ALL:
     for k, v in iteritems(CLASS_INPUT):
         models, state = v
         CLASS_INPUT[k] = (models, 'power')
+if VALGRIND_MODE:
+    for k, v in iteritems(CLASS_INPUT):
+        if k == 'Output_spectra':
+            continue
+        models, state = v
+        CLASS_INPUT[k] = (models, 'normal')
 
 INPUTPOWER = []
 INPUTNORMAL = [{}]
+INPUTONLYFULL = [{}]
 for key, value in list(CLASS_INPUT.items()):
     models, state = value
     if state == 'power':
         INPUTPOWER.append([{}]+models)
+    elif state == 'onlyfull':
+        INPUTONLYFULL.extend(models)
     else:
         INPUTNORMAL.extend(models)
 
@@ -199,6 +228,10 @@
                 temp_dict.update(elem)
             DICTARRAY.append(temp_dict)
 
+    for onlyfull in INPUTONLYFULL:
+        temp_dict = onlyfull.copy()
+        temp_dict.update(FULL_SETTINGS)
+        DICTARRAY.append(temp_dict)
 
 TUPLE_ARRAY = []
 for e in DICTARRAY:
@@ -212,7 +245,7 @@ def powerset(iterable):
         itertools.combinations(xs, n) for n in range(1, len(xs)+1))
 
 def custom_name_func(testcase_func, param_num, param):
-    special_keys = ['N_ncdm']
+    special_keys = ['N_ncdm','Pk_ini_type','ic','Neff','sigma8','S8','theta_s_100','tau_reio']
     somekeys = []
     for key in param.args[0]:
         if key in special_keys:
@@ -435,9 +468,9 @@ def has_incompatible_input(self):
             if not has_tensor(self.scenario):
                 should_fail = True
 
-        # If we have specified non linear, we must have some form of
+        # If we have specified non_linear, we must have some form of
         # perturbations output.
-        if 'non linear' in self.scenario:
+        if 'non_linear' in self.scenario:
             if 'output' not in self.scenario:
                 should_fail = True
 
@@ -460,7 +493,7 @@ def has_incompatible_input(self):
 
         # If we use inflation module, we must have scalar modes,
         # tensor modes, no vector modes and we should only have adiabatic IC:
-        if 'P_k_ini type' in self.scenario and self.scenario['P_k_ini type'].find('inflation') != -1:
+        if 'Pk_ini_type' in self.scenario and self.scenario['Pk_ini_type'].find('inflation') != -1:
             if 'modes' not in self.scenario:
                 should_fail = True
             else:
@@ -472,6 +505,8 @@ def has_incompatible_input(self):
                     should_fail = True
             if 'ic' in self.scenario and self.scenario['ic'].find('i') != -1:
                 should_fail = True
+            if 'non_linear' in self.scenario and self.scenario['non_linear'].find('hmcode') != -1:
+                should_fail = True
 
 
         return should_fail
@@ -545,6 +580,36 @@ def compare_output(self, reference, reference_name, candidate, candidate_name, r
                     self.pk_faulty_plot(k, reference_pk, reference_name, candidate_pk, candidate_name, rtol_pk)
                     status_pass = False
 
+        if 'output' in self.scenario and self.scenario['output'].find('mPk') != -1:
+            ref_sigma8 = reference.sigma8()
+            can_sigma8 = candidate.sigma8()
+            ref_sigma = reference.sigma(8, 0, h_units=True)
+            can_sigma = candidate.sigma(8, 0, h_units=True)
+            try:
+                np.testing.assert_allclose([ref_sigma8, ref_sigma, can_sigma8],
+                                           [can_sigma8, can_sigma, can_sigma],
+                                           rtol=rtol_pk,
+                                           atol=COMPARE_PK_ABSOLUTE_ERROR)
+            except AssertionError:
+                status_pass = False
+
+        if 'sigma8' in self.scenario:
+            if not np.isclose(self.scenario['sigma8'], candidate.sigma8(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.sigma8(), reference.sigma8(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+
+        if 'Neff' in self.scenario:
+            if not np.isclose(self.scenario['Neff'], candidate.Neff(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.Neff(), reference.Neff(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+
+        if 'theta_s_100' in self.scenario:
+            if not np.isclose(self.scenario['theta_s_100'], candidate.theta_s_100(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.theta_s_100(), reference.theta_s_100(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
         return status_pass
 
     def store_ini_file(self, path):
diff --git a/source/fourier.c b/source/fourier.c
index 392b4983c..1d2186ee1 100644
--- a/source/fourier.c
+++ b/source/fourier.c
@@ -1776,7 +1776,7 @@ int fourier_indices(
 
   /** - get list of k values */
 
-  class_call(fourier_get_k_list(ppr,ppt,pfo),
+  class_call(fourier_get_k_list(ppr,ppm,ppt,pfo),
              pfo->error_message,
              pfo->error_message);
 
@@ -1857,6 +1857,7 @@ int fourier_indices(
 
 int fourier_get_k_list(
                        struct precision *ppr,
+                       struct primordial * ppm,
                        struct perturbations * ppt,
                        struct fourier * pfo
                        ) {
@@ -1881,6 +1882,7 @@ int fourier_get_k_list(
                "could not reach extrapolated value k = %.10e starting from k = %.10e with k_per_decade of %.10e in _MAX_NUM_INTERPOLATION_=%i steps",
                ppr->hmcode_max_k_extra,k_max,ppr->k_per_decade_for_pk,_MAX_NUM_EXTRAPOLATION_
                );
+
     pfo->k_size_extra = pfo->k_size+index_k;
   }
   /** - otherwise, same number of values as in perturbation module */
@@ -1906,6 +1908,13 @@ int fourier_get_k_list(
     pfo->ln_k[index_k] = log(k) + exponent*log(10.);
   }
 
+  class_test(pfo->k[pfo->k_size_extra-1]>exp(ppm->lnk[ppm->lnk_size-1]) && ppm->primordial_spec_type != analytic_Pk,
+             pfo->error_message,
+             "Setting the output to HMcode with a large 'hmcode_max_k_extra' and using the primordial spectrum to not analytic is incompatible. Either use the analytic power spectrum or set a smaller 'hmcode_max_k_extra' (k_max_hmcode=%.5e , k_max_primordial=%.5e)",
+             pfo->k[pfo->k_size_extra-1],
+             exp(ppm->lnk[ppm->lnk_size-1])
+             )
+
   return _SUCCESS_;
 }
 
diff --git a/source/input.c b/source/input.c
index 581946223..f35aa82cc 100644
--- a/source/input.c
+++ b/source/input.c
@@ -412,6 +412,9 @@ int input_read_from_file(struct file_content * pfc,
   class_read_int("input_verbose",input_verbose);
   if (input_verbose >0) printf("Reading input parameters\n");
 
+  /** -- Special setting of parameter, before anything else: did shooting fail? */
+  pba->shooting_failed = _FALSE_;
+
   /** Find out if shooting necessary and, eventually, shoot and initialize
       read parameters */
   class_call(input_shooting(pfc,ppr,pba,pth,ppt,ptr,ppm,phr,pfo,ple,psd,pop,
@@ -641,10 +644,9 @@ int input_shooting(struct file_content * pfc,
 
       /* We can do 1 dimensional root finding */
       if (input_verbose > 0) {
-        fprintf(stdout,
-                "Computing unknown input parameter '%s' using input parameter '%s'\n",
-                fzw.fc.name[fzw.unknown_parameters_index[0]],
-                target_namestrings[fzw.target_name[0]]);
+        printf("Computing unknown input parameter '%s' using input parameter '%s'\n",
+               fzw.fc.name[fzw.unknown_parameters_index[0]],
+               target_namestrings[fzw.target_name[0]]);
       }
 
       /* If shooting fails, postpone error to background module to play nice with MontePython. */
@@ -662,9 +664,14 @@ int input_shooting(struct file_content * pfc,
       // precision of around 1e-16, so 1e-20 should be good enough for the shooting
       class_sprintf(fzw.fc.value[fzw.unknown_parameters_index[0]],"%.20e",xzero);
       if (input_verbose > 0) {
-        fprintf(stdout," -> found '%s = %s'\n",
-                fzw.fc.name[fzw.unknown_parameters_index[0]],
-                fzw.fc.value[fzw.unknown_parameters_index[0]]);
+        if (shooting_failed == _FALSE_){
+          printf(" -> found '%s = %s'\n",
+                  fzw.fc.name[fzw.unknown_parameters_index[0]],
+                  fzw.fc.value[fzw.unknown_parameters_index[0]]);
+        }
+        else{
+          printf("Shooting failed! Aborting...\n");
+        }
       }
 
     }
@@ -673,7 +680,7 @@ int input_shooting(struct file_content * pfc,
 
       /* We need to do multidimensional root finding */
       if (input_verbose > 0) {
-        fprintf(stdout,"Computing unknown input parameters\n");
+        printf("Computing unknown input parameters\n");
       }
 
       /* Allocate local variables */
@@ -710,9 +717,14 @@ int input_shooting(struct file_content * pfc,
         class_sprintf(fzw.fc.value[fzw.unknown_parameters_index[counter]],
                 "%.20e",x_inout[counter]);
         if (input_verbose > 0) {
-          fprintf(stdout," -> found '%s = %s'\n",
-                  fzw.fc.name[fzw.unknown_parameters_index[counter]],
-                  fzw.fc.value[fzw.unknown_parameters_index[counter]]);
+          if (shooting_failed == _FALSE_){
+            printf(" -> found '%s = %s'\n",
+                    fzw.fc.name[fzw.unknown_parameters_index[counter]],
+                    fzw.fc.value[fzw.unknown_parameters_index[counter]]);
+          }
+          else{
+            printf("Shooting failed! Aborting...\n");
+          }
         }
       }
 
@@ -721,8 +733,8 @@ int input_shooting(struct file_content * pfc,
       free(dxdF);
     }
 
-    if (input_verbose > 1) {
-      fprintf(stdout,"Shooting completed using %d function evaluations\n",fevals);
+    if (input_verbose > 1 && shooting_failed == _FALSE_) {
+      printf("Shooting completed using %d function evaluations\n",fevals);
     }
 
     /** Set status of shooting */
@@ -757,7 +769,7 @@ int input_shooting(struct file_content * pfc,
     /* Create file content structure with additional entries */
     class_call(parser_extend(pfc, 1, errmsg),
                errmsg,errmsg);
-    
+
     class_call(parser_init_from_pfc(pfc, &(fzw.fc), errmsg),
                errmsg,errmsg);
 
@@ -790,10 +802,9 @@ int input_shooting(struct file_content * pfc,
 
     /* Print to the user */
     if (input_verbose > 0) {
-      fprintf(stdout,
-              "Computing unknown input parameter '%s' using input parameter '%s'\n",
-              (flag1 ==_TRUE_?"sigma8":"S8"),
-              "A_s");
+      printf("Computing unknown input parameter '%s' using input parameter '%s'\n",
+             (flag1 ==_TRUE_?"sigma8":"S8"),
+             "A_s");
     }
 
     /* Set a guess for A_s from LCDM (doesn't need to be super accurate) */
@@ -820,11 +831,11 @@ int input_shooting(struct file_content * pfc,
     /* Store the derived value with high enough accuracy */
     class_sprintf(fzw.fc.value[pfc->size - 1],"%.20e",A_s);
     if (input_verbose > 0) {
-      fprintf(stdout," -> found '%s = %s'\n",
-              fzw.fc.name[pfc->size - 1],
-              fzw.fc.value[pfc->size - 1]);
+      printf(" -> found '%s = %s'\n",
+             fzw.fc.name[pfc->size - 1],
+             fzw.fc.value[pfc->size - 1]);
     }
-  
+
     parser_copy(&(fzw.fc), pfc, pfc->size - 1, pfc->size);
 
     /** Free arrays allocated */
@@ -868,7 +879,6 @@ int input_needs_shooting_for_target(struct file_content * pfc,
   case Omega_scf:
   case Omega_ini_dcdm:
   case omega_ini_dcdm:
-  case Neff:
     /* Check that Omega's or omega's are nonzero: */
     if (target_value == 0.)
       *needs_shooting = _FALSE_;
@@ -926,8 +936,6 @@ int input_find_root(double *xzero,
   /* Try fifteen times to go above and below the root (i.e. where shooting succeeds) */
   for (iter=1; iter<=15; iter++){
     x2 = x1 - dx;
-    // Early stopping if the two points are infinitessimaly close together, i.e. the answer is already reached
-    if(abs(x2/x1-1) < tol_x_rel){*xzero = x1; return _SUCCESS_;}
     /* Try three times to get a 'reasonable' value, i.e. no CLASS error */
     for (iter2=1; iter2 <= 3; iter2++) {
       return_function = input_fzerofun_1d(x2, pfzw, &f2, errmsg);
@@ -1154,6 +1162,9 @@ int input_get_guess(double *xguess,
   /* Cheat to read only known parameters: */
   pfzw->fc.size -= pfzw->target_size;
 
+  /* Assume for now shooting did not fail */
+  ba.shooting_failed = _FALSE_;
+
   class_call(input_read_precisions(&(pfzw->fc),&pr,&ba,&th,&pt,&tr,&pm,&hr,&fo,&le,&sd,&op,
                                    errmsg),
              errmsg,
@@ -1319,6 +1330,9 @@ int input_try_unknown_parameters(double * unknown_parameter,
   int param;
   short compute_sigma8 = _FALSE_;
 
+  /* Assume for now shooting did not fail */
+  ba.shooting_failed = _FALSE_;
+
   pfzw = (struct fzerofun_workspace *) voidpfzw;
   /** Read input parameters */
   // This needs to be done with enough accuracy. A standard double has a relative
@@ -2077,6 +2091,9 @@ int input_read_parameters_general(struct file_content * pfc,
 
   }
 
+  /** 4.c) Do we want matter and baryon+CDM sources in current gauge, instead of automatic conversion to gauge-invariant variables? */
+  class_read_flag("matter_source_in_current_gauge",ppt->has_matter_source_in_current_gauge);
+
   /** 5) h in [-] and H_0/c in [1/Mpc = h/2997.9 = h*10^5/c] */
   /* Read */
   class_call(parser_read_double(pfc,"H0",&param1,&flag1,errmsg),
@@ -2430,7 +2447,7 @@ int input_read_parameters_species(struct file_content * pfc,
       pba->Omega0_ur = param3/pba->h/pba->h;
     }
   }
-  class_test(pba->Omega0_ur<0,errmsg,"You cannot set the density of ultra-relativistic relics (dark radiation/neutrinos) to negative values.");
+  class_test(pba->Omega0_ur<0,errmsg,"You cannot set the density of ultra-relativistic relics (dark radiation/neutrinos) to negative values. You might have input a total Neff smaller than what your massive neutrinos require minimally (around 1.02 * N_ncdm * deg_ncdm).");
 
   /** 3.a) Case of non-standard properties */
   /* Read */
@@ -2543,18 +2560,28 @@ int input_read_parameters_species(struct file_content * pfc,
                errmsg,
                errmsg);
 
-    /** 5.d) Mass or Omega of each ncdm species */
+    /** 5.d) Mass and/or Omega of each ncdm species */
     /* Read */
     class_read_list_of_doubles_or_default("m_ncdm",pba->m_ncdm_in_eV,0.0,N_ncdm);
+    for (n=0; n<N_ncdm; n++){
+      class_test(pba->m_ncdm_in_eV[n]<0,
+                 errmsg,
+                 "You entered a negative non-CDM mass m_ncdm[%d], which makes no sense. This error was not caught in previous CLASS versions because the mass is always squared in the code, so CLASS returned the exact same results form +m_ncdm and -m_ncdm. If you want to define an 'effective negative neutrino mass' in the sense of e.g. 2405.00836 or 2407.10965, you can implement it in a python script following e.g. eq.(3) of 2407.10965",n);
+    }
+
     class_read_list_of_doubles_or_default("Omega_ncdm",pba->Omega0_ncdm,0.0,N_ncdm);
+    // the name M_ncdm is borrowed temporarily to store omega_ncdm
     class_read_list_of_doubles_or_default("omega_ncdm",pba->M_ncdm,0.0,N_ncdm);
     for (n=0; n<N_ncdm; n++){
       if (pba->M_ncdm[n]!=0.0){
         /* Test */
         class_test(pba->Omega0_ncdm[n]!=0,errmsg,
                    "You can only enter one of 'Omega_ncdm' or 'omega_ncdm' for ncdm species %d.",n);
-        /* Complete set of parameters */
+        /* Complete set of parameters: if the user passed either
+           Omega_ncdm or omega_ncdm, now it's stored anyway as
+           Omega_0_ncdm */
         pba->Omega0_ncdm[n] = pba->M_ncdm[n]/pba->h/pba->h;
+        // the name M_ncdm is now available for its real destination
       }
       /* Set default value
          this is the right place for passing the default value of the mass
@@ -5658,6 +5685,8 @@ int input_default_params(struct background *pba,
   ppt->gauge=synchronous;
   /** 4.b) N-body gauge */
   ppt->has_Nbody_gauge_transfers = _FALSE_;
+  /** 4.c) keep delta_m, theta_m, delta_cb, theta_cb in current gauge */
+  ppt->has_matter_source_in_current_gauge = _FALSE_;
 
   /** 5) Hubble parameter */
   pba->h = 0.67810;
@@ -5826,8 +5855,6 @@ int input_default_params(struct background *pba,
   pba->phi_prime_ini_scf = 1;          //     factors of the radiation attractor values
   /** 9.b.3) Tuning parameter */
   pba->scf_tuning_index = 0;
-  /** 9.b.4) Shooting parameter */
-  pba->shooting_failed = _FALSE_;
 
   /**
    * Deafult to input_read_parameters_heating
diff --git a/source/perturbations.c b/source/perturbations.c
index d88662b21..9eaf301ce 100644
--- a/source/perturbations.c
+++ b/source/perturbations.c
@@ -866,9 +866,9 @@ int perturbations_init(
                ppt->z_max_pk,
                pth->z_rec);
 
-    class_test(ppt->has_source_delta_m == _TRUE_,
+    class_test((ppt->has_source_delta_m == _TRUE_) && (ppt->has_matter_source_in_current_gauge == _FALSE_),
                ppt->error_message,
-               "You requested a very high z_pk=%e, higher than z_rec=%e. This works very well when you ask only transfer functions, e.g. with 'output=mTk' or 'output=mTk,vTk'. But if you need the total matter (e.g. with 'mPk', 'dCl', etc.) there is an issue with the calculation of delta_m at very early times. By default, delta_m is a gauge-invariant variable (the density fluctuation in comoving gauge) and this quantity is hard to get accurately at very early times. The solution is to define delta_m as the density fluctuation in the current gauge, synchronous or newtonian. For the moment this must be done manually by commenting the line 'ppw->delta_m += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2;' in perturbations_sources(). In the future there will be an option for doing it in an easier way.",
+               "You requested a very high z_pk=%e, higher than z_rec=%e. This works very well when you ask only transfer functions, e.g. with 'output=mTk' or 'output=mTk,vTk'. But if you need the total matter (e.g. with 'mPk', 'dCl', etc.) there is an issue with the calculation of delta_m (or delta_cb) at very early times. By default, delta_m and delta_cb are expressed by the code as gauge-invariant variables (the density fluctuation in comoving gauge). This quantity is hard to get accurately at very early times. The solution is to define delta_m and delta_cb as density fluctuations in the current gauge, e.g. synchronous or newtonian. This is done by setting the input flag 'matter_source_in_current_gauge' to 'yes'.",
                ppt->z_max_pk,
                pth->z_rec);
 
@@ -6610,31 +6610,33 @@ int perturbations_einstein(
        gauge-independent variables (you could comment this out if you
        really want gauge-dependent results) */
 
-    if (ppt->has_source_delta_m == _TRUE_) {
-      ppw->delta_m += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2;
-      // note: until 2.4.3 there was a typo, the factor was (-2 H'/H) instead
-      // of (3 aH). There is the same typo in the CLASSgal paper
-      // 1307.1459v1,v2,v3. It came from a confusion between (1+w_total)
-      // and (1+w_matter)=1 [the latter is the relevant one here].
-      //
-      // note2: at this point this gauge-invariant variable is only
-      // valid if all matter components are pressureless and
-      // stable. This relation will be generalized soon to the case
-      // of decaying dark matter.
-    }
+    if (ppt->has_matter_source_in_current_gauge == _FALSE_) {
+      if (ppt->has_source_delta_m == _TRUE_) {
+        ppw->delta_m += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2;
+        // note: until 2.4.3 there was a typo, the factor was (-2 H'/H) instead
+        // of (3 aH). There is the same typo in the CLASSgal paper
+        // 1307.1459v1,v2,v3. It came from a confusion between (1+w_total)
+        // and (1+w_matter)=1 [the latter is the relevant one here].
+        //
+        // note2: at this point this gauge-invariant variable is only
+        // valid if all matter components are pressureless and
+        // stable. This relation will be generalized soon to the case
+        // of decaying dark matter.
+      }
 
-    if (ppt->has_source_delta_cb == _TRUE_) {
-      ppw->delta_cb += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_cb/k2;//check gauge transformation
-    }
+      if (ppt->has_source_delta_cb == _TRUE_) {
+        ppw->delta_cb += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_cb/k2;//check gauge transformation
+      }
 
-    if (ppt->has_source_theta_m == _TRUE_) {
-      if  (ppt->gauge == synchronous) {
-        ppw->theta_m += ppw->pvecmetric[ppw->index_mt_alpha]*k2;
+      if (ppt->has_source_theta_m == _TRUE_) {
+        if  (ppt->gauge == synchronous) {
+          ppw->theta_m += ppw->pvecmetric[ppw->index_mt_alpha]*k2;
+        }
       }
-    }
-    if (ppt->has_source_theta_cb == _TRUE_){
-      if  (ppt->gauge == synchronous) {
-        ppw->theta_cb += ppw->pvecmetric[ppw->index_mt_alpha]*k2; //check gauge transformation
+      if (ppt->has_source_theta_cb == _TRUE_){
+        if  (ppt->gauge == synchronous) {
+          ppw->theta_cb += ppw->pvecmetric[ppw->index_mt_alpha]*k2; //check gauge transformation
+        }
       }
     }
   }