diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 7180778a2d005acad6983a1c62b89ea8f9c71a56..13cf380f2973b37704986c81847d6ba3ef1ba5d5 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -22,7 +22,7 @@ Tasks:
 import math
 import numpy as np
 from typing import Optional, Union
-from binarycpython.utils.useful_funcs import calc_period_from_sep
+from binarycpython.utils.useful_funcs import calc_period_from_sep, calc_sep_from_period
 
 ###
 # File containing probability distributions
@@ -1425,12 +1425,6 @@ def Moe_de_Stefano_2017_pdf(options):
 
     prob = [] # Value that we will return
 
-    # Separation of the inner binary
-    options['sep'] = calc_sep_from_period(options['M1'], options['M2'], options['P'])
-
-    # Total mass inner binary:
-    options['M1+M2'] = options['M1'] + options['M2']
-
     multiplicity = options['multiplicity']
     if not options.get('multiplicity', None):
         multiplicity = 1
@@ -1446,7 +1440,7 @@ def Moe_de_Stefano_2017_pdf(options):
 
     ############################################################
     # multiplicity fraction    
-    prob.append(Moe_de_Stefano_2017_multiplicity_fractions(opts)[multiplicity-1])
+    prob.append(Moe_de_Stefano_2017_multiplicity_fractions(options)[multiplicity-1])
 
     ############################################################
     # always require an IMF for the primary star
@@ -1458,6 +1452,13 @@ def Moe_de_Stefano_2017_pdf(options):
     prob.append(Kroupa2001(options['M1']) * options['M1'])
 
     if(multiplicity >= 2):
+
+        # Separation of the inner binary
+        options['sep'] = calc_sep_from_period(options['M1'], options['M2'], options['P'])
+
+        # Total mass inner binary:
+        options['M1+M2'] = options['M1'] + options['M2']
+
         # binary, triple or quadruple system
         if not Moecache.get("rinterpolator_log10P", None):
             Moecache['rinterpolator_log10P'] = py_rinterpolate.Rinterpolate(
diff --git a/binarycpython/utils/functions.py b/binarycpython/utils/functions.py
index e96fc6b36ff7a2eb5fa2cc8e67d52b4c929a4311..4fcb1b7fa435063276b6e1a9369ed50b7032c951 100644
--- a/binarycpython/utils/functions.py
+++ b/binarycpython/utils/functions.py
@@ -1257,6 +1257,96 @@ def merge_dicts(dict_1: dict, dict_2: dict) -> dict:
     #
     return new_dict
 
+def update_dicts(dict_1: dict, dict_2: dict) -> dict:
+    """
+    Function to update dict_1 with values of dict_2 in a recursive way.
+
+    Behaviour:
+
+    When dict keys are only present in one of either:
+        - we just add the content to the new dict
+
+    When dict keys are present in both, we decide based on the value types how to combine them:
+        - value of dict2 will be taken
+
+    Args:
+        dict_1: first dictionary
+        dict_2: second dictionary
+
+    Returns:
+        New dictionary with Updated values
+
+    """
+
+    # Set up new dict
+    new_dict = {}
+
+    #
+    keys_1 = dict_1.keys()
+    keys_2 = dict_2.keys()
+
+    # Find overlapping keys of both dicts
+    overlapping_keys = set(keys_1).intersection(set(keys_2))
+
+    # Find the keys that are unique
+    unique_to_dict_1 = set(keys_1).difference(set(keys_2))
+    unique_to_dict_2 = set(keys_2).difference(set(keys_1))
+
+    # Add the unique keys to the new dict
+    for key in unique_to_dict_1:
+        # If these items are ints or floats, then just put them in
+        if isinstance(dict_1[key], (float, int)):
+            new_dict[key] = dict_1[key]
+        # Else, to be safe we should deepcopy them
+        else:
+            copy_dict = copy.deepcopy(dict_1[key])
+            new_dict[key] = copy_dict
+
+    for key in unique_to_dict_2:
+        # If these items are ints or floats, then just put them in
+        if isinstance(dict_2[key], (float, int)):
+            new_dict[key] = dict_2[key]
+        # Else, to be safe we should deepcopy them
+        else:
+            copy_dict = copy.deepcopy(dict_2[key])
+            new_dict[key] = copy_dict
+
+    # Go over the common keys:
+    for key in overlapping_keys:
+
+        # See whether the types are actually the same
+        if not type(dict_1[key]) is type(dict_2[key]):
+            # Exceptions:
+            if (type(dict_1[key]) in [int, float]) and (
+                type(dict_2[key]) in [int, float]
+            ):
+                new_dict[key] = dict_2[key]
+
+            else:
+                print(
+                    "Error key: {} value: {} type: {} and key: {} value: {} type: {} are not of the same type and cannot be merged".format(
+                        key,
+                        dict_1[key],
+                        type(dict_1[key]),
+                        key,
+                        dict_2[key],
+                        type(dict_2[key]),
+                    )
+                )
+                raise ValueError
+
+        # Here we check for the cases that we want to explicitly catch. Ints will be added,
+        # floats will be added, lists will be appended (though that might change) and dicts will be
+        # dealt with by calling this function again.
+        else:
+            # dicts
+            if isinstance(dict_1[key], dict) and isinstance(dict_2[key], dict):
+                new_dict[key] = update_dicts(dict_1[key], dict_2[key])
+            else:
+                new_dict[key] = dict_2[key]
+
+    #
+    return new_dict
 
 def extract_ensemble_json_from_string(binary_c_output: str) -> dict:
     """
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 9e8836f02c0b972a6d6671dc22f33c37dc249495..9f44a0931f67b7e03bec5b00a4839f05b8b28c6b 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -31,10 +31,12 @@ import json
 import datetime
 import time
 import uuid
+import copy
 import logging
 import argparse
 import subprocess
 import importlib.util
+
 from typing import Union, Any
 from pathos.helpers import mp as pathos_multiprocess
 
@@ -57,6 +59,7 @@ from binarycpython.utils.functions import (
     verbose_print,
     binarycDecoder,
     merge_dicts,
+    update_dicts,
     BinaryCEncoder,
     extract_ensemble_json_from_string,
 )
@@ -69,8 +72,19 @@ from binarycpython.utils.functions import (
 #     get_python_details,
 # )
 
+from binarycpython.utils.distribution_functions import (
+    Moecache,
+    LOG_LN_CONVERTER
+)
+
 from binarycpython import _binary_c_bindings
-import copy
+
+import types
+def imports():
+    for name, val in globals().items():
+        if isinstance(val, types.ModuleType):
+            yield val.__name__
+
 
 class Population:
     """
@@ -89,6 +103,7 @@ class Population:
         # Setting stuff will check against the defaults to see if the input is correct.
         self.defaults = get_defaults()
         self.cleaned_up_defaults = self._cleanup_defaults()
+        self.available_keys = list(self.defaults.keys())
 
         # make the input dictionary
         self.bse_options = {}  # bse_options is just empty.
@@ -317,7 +332,8 @@ class Population:
 
         number = len(self.grid_options["_grid_variables"])
         for grid_variable in self.grid_options["_grid_variables"]:
-            print(self.grid_options["_grid_variables"][grid_variable]['grid_variable_number'])
+            if self.grid_options["_grid_variables"][grid_variable]['grid_variable_number'] == number-1:
+                return grid_variable
 
     def add_grid_variable(
         self,
@@ -327,8 +343,9 @@ class Population:
         resolution: str,
         spacingfunc: str,
         probdist: str,
-        dphasevol: str,
+        dphasevol: Union[str, int],
         parameter_name: str,
+        gridtype: str,
         precode: Union[str, None] = None,
         condition: Union[str, None] = None,
     ) -> None:
@@ -366,11 +383,13 @@ class Population:
                 FUnction determining the probability that gets asigned to the sampled parameter
                 example: probdist = 'Kroupa2001(M_1)*M_1'
             dphasevol:
-                part of the parameter space that the total probability is calculated with
+                part of the parameter space that the total probability is calculated with. Put to -1 if you want to ignore any dphasevol calculations and set the value to 1
                 example: dphasevol = 'dlnm1'
             condition:
                 condition that has to be met in order for the grid generation to continue
                 example: condition = 'self.grid_options['binary']==1'
+            gridtype:
+                Method on how the value range is sampled. Can be either edge (steps starting at the lower edge of the value range) or centered (steps starting at lower edge + 0.5 * stepsize). 
         """
 
         # Add grid_variable
@@ -385,6 +404,7 @@ class Population:
             "dphasevol": dphasevol,
             "parameter_name": parameter_name,
             "condition": condition,
+            "gridtype": gridtype,
             "grid_variable_number": len(self.grid_options["_grid_variables"]),
         }
 
@@ -919,6 +939,9 @@ class Population:
                 "\tUsing persistent_data memaddr: {}".format(persistent_data_memaddr), self.grid_options["verbosity"], 0
             )
 
+        #
+        self._generate_grid_code(dry_run=False)
+
         # apparently we have to re-load this for every process, otherwise NameErrors arise (seems like a bug but I'm not sure)
         self._load_grid_function()
 
@@ -972,6 +995,15 @@ class Population:
 
                     start_runtime_binary_c = time.time()
 
+                    # In the first system, explicitly check all the keys that are passed to see if they match the keys known to binary_c.
+                    # Won't do that every system cause that is a bit of a waste of computing time. 
+                    if localcounter == 0:
+                        for key in full_system_dict.keys():
+                            if not key in self.available_keys:
+                                print("Error: Found a parameter unknown to binary_c: {}. Abort".format(key))
+                                raise ValueError
+
+                    # TODO: build flag to actually evolve the system
                     # Evolve the system
                     self._evolve_system_mp(full_system_dict)
 
@@ -1223,12 +1255,15 @@ class Population:
                 "_start_time_evolution"
             ] = time.time()  # Setting start time of grid
 
-            #
-            self._generate_grid_code(dry_run=False)
+            # # Making sure the loaded grid code isn't lingering in the main PID
+            # self._generate_grid_code(dry_run=False)
 
-            #
-            self._load_grid_function()
+            # #
+            # self._load_grid_function()
 
+            # 
+            self.grid_options["_system_generator"] = None
+            
         # Source file
         elif self.grid_options["evolution_type"] == "source_file":
             #######################
@@ -1322,6 +1357,7 @@ class Population:
         # TODO: add sensible description to this function.
         # TODO: Check whether all the probability and phasevol values are correct.
         # TODO: import only the necessary packages/functions
+        # TODO: Put all the masses, eccentricities and periods in there already
 
         Results in a generated file that contains a system_generator function.
         """
@@ -1381,6 +1417,22 @@ class Population:
         code_string += indent * depth + "phasevol = 1\n"
         code_string += indent * depth + "\n"
 
+        # Setting up the system parameters
+        code_string += indent * depth + "m1 = None\n"
+        code_string += indent * depth + "m2 = None\n"
+        code_string += indent * depth + "m3 = None\n"
+        code_string += indent * depth + "m4 = None\n"
+
+        code_string += indent * depth + "per = None\n"
+        code_string += indent * depth + "per2 = None\n"
+        code_string += indent * depth + "per3 = None\n"
+
+        code_string += indent * depth + "eccentricity = None\n"
+        code_string += indent * depth + "eccentricity2 = None\n"
+        code_string += indent * depth + "eccentricity3 = None\n"
+        code_string += indent * depth + "\n"
+
+
         # Prepare the probability
         code_string += indent * depth + "# setting probability lists\n"
         for grid_variable_el in sorted(
@@ -1408,37 +1460,6 @@ class Population:
             print("Constructing/adding: {}".format(grid_variable_el[0]))
             grid_variable = grid_variable_el[1]
 
-            #################################################################################
-            # Check condition and generate forloop
-
-            # If the grid variable has a condition, write the check and the action
-            if grid_variable["condition"]:
-                # Add comment
-                code_string += (
-                    indent * depth
-                    + "# Condition for {}".format(grid_variable["parameter_name"])
-                    + "\n"
-                )
-
-                # Add condition check
-                code_string += (
-                    indent * depth
-                    + "if not {}:".format(grid_variable["condition"])
-                    + "\n"
-                )
-
-                # Add condition failed action:
-                code_string += (
-                    indent * (depth + 1)
-                    + 'print("Condition for {} not met!")'.format(
-                        grid_variable["parameter_name"]
-                    )
-                    + "\n"
-                )
-                code_string += indent * (depth + 1) + "raise ValueError" + "\n"
-
-                # Add some whiteline
-                code_string += indent * (depth + 1) + "\n"
 
             #########################
             # Setting up the forloop
@@ -1460,13 +1481,19 @@ class Population:
             # TODO: Make clear that the phasevol only works good
             # TODO: add option to ignore this phasevol calculation and set it to 1
             #   if you sample linearly in that thing.
-            code_string += (
-                indent * depth
-                + "phasevol_{} = sampled_values_{}[1]-sampled_values_{}[0]".format(
-                    grid_variable["name"], grid_variable["name"], grid_variable["name"]
+            if not grid_variable['dphasevol']==-1: # A method to turn off this calculation and allow a phasevol = 1
+                code_string += (
+                    indent * depth
+                    + "phasevol_{} = sampled_values_{}[1]-sampled_values_{}[0]".format(
+                        grid_variable["name"], grid_variable["name"], grid_variable["name"]
+                    )
+                    + "\n"
+                )
+            else:
+                code_string += (
+                    indent * depth
+                    + "phasevol_{} = 1\n".format(grid_variable["name"])
                 )
-                + "\n"
-            )
 
             # # Some print statement
             # code_string += (
@@ -1485,6 +1512,38 @@ class Population:
                 + "\n"
             )
 
+            #################################################################################
+            # Check condition and generate forloop
+
+            # If the grid variable has a condition, write the check and the action
+            if grid_variable["condition"]:
+                # Add comment
+                code_string += (
+                    indent * (depth + 1)
+                    + "# Condition for {}".format(grid_variable["parameter_name"])
+                    + "\n"
+                )
+
+                # Add condition check
+                code_string += (
+                    indent * (depth + 1)
+                    + "if not {}:".format(grid_variable["condition"])
+                    + "\n"
+                )
+
+                # Add condition failed action:
+                code_string += (
+                    indent * (depth + 2)
+                    + 'print("Condition for {} not met!")'.format(
+                        grid_variable["parameter_name"]
+                    )
+                    + "\n"
+                )
+                code_string += indent * (depth + 2) + "continue" + "\n"
+
+                # Add some whiteline
+                code_string += indent * (depth + 1) + "\n"
+
             #########################
             # Setting up pre-code and value in some cases
             # Add pre-code
@@ -2685,6 +2744,8 @@ class Population:
             if not options['file'].endswith('.json'):
                 print("Provided filename does not end with .json")
 
+        # TODO: put this in defaults
+
         # Set default values
         default_options = {
             'resolutions':
@@ -2803,7 +2864,7 @@ class Population:
         }
 
         # Take the option dictionary that was given and override.
-        options = {**default_options, **options}
+        options = update_dicts(default_options, options)
 
         # Write to a file
         with open('/tmp/moeopts.dat', 'w') as f:
@@ -3067,15 +3128,16 @@ class Population:
         # Multiplicity
         self.add_grid_variable(
                 name='multiplicity',
+                parameter_name='multiplicity',
                 longname='multiplicity',
-                range=[1, max_multiplicity],
+                valuerange=[1, max_multiplicity],
                 resolution=1,
-                spacingfunc='number(1.0)',
-                precode='self.grid_options["multiplicity"] = multiplicity',
-                condition='if({}[{}-1] > 0'.format(str(options['multiplicity_modulator'])),
+                spacingfunc='range(1, 5)',
+                precode='self.grid_options["multiplicity"] = multiplicity; options={}'.format(options),
+                condition='({}[multiplicity-1] > 0)'.format(str(options['multiplicity_modulator'])),
                 gridtype='edge',
-                noprobdist=1,
-                nophasevol=1
+                dphasevol=-1,
+                probdist=1,
         )
 
         ############################################################
@@ -3087,17 +3149,18 @@ class Population:
         # log-spaced m1 with given resolution
         self.add_grid_variable(
             name='lnm1',
+            parameter_name='M_1',
             longname='Primary mass',
             resolution="options['resolutions']['M'][0]",
-            spacingsfunc="const(np.log({}), np.log({}),0 {})".format(options['ranges']['M'][0], options['ranges']['M'][1], options['resolutions']['M'][0]),
-            range=[
+            spacingfunc="const(np.log({}), np.log({}), {})".format(options['ranges']['M'][0], options['ranges']['M'][1], options['resolutions']['M'][0]),
+            valuerange=[
                 'np.log({})'.format(options['ranges']['M'][0]),
                 'np.log({})'.format(options['ranges']['M'][1])
             ],
             gridtype='centred',
             dphasevol='dlnm1',
-            precode='m1 = np.exp(lnm1); options["M1"]=m1',
-            probdist='self.Moe_de_Stefano_2017_pdf(options) if multiplicity == 1 else 1'
+            precode='M_1 = np.exp(lnm1); options["M1"]=M_1',
+            probdist='Moe_de_Stefano_2017_pdf({{{}, {}, {}}}) if multiplicity == 1 else 1'.format(str(options)[1:-1], "'multiplicity': multiplicity", "'M1': M_1")
         )
 
         # Go to higher multiplicities
@@ -3105,10 +3168,11 @@ class Population:
             # binaries: period
             self.add_grid_variable(
                 name='log10per',
+                parameter_name='orbital_period',
                 longname='log10(Orbital_Period)',
                 resolution=options['resolutions']['logP'][0],
                 probdist=1.0,
-                condition='self.grid_options["grid_options"]["multiplicity"] >= 2',
+                condition='(self.grid_options["grid_options"]["multiplicity"] >= 2)',
                 branchpoint=1,
                 gridtype='centred',
                 dphasevol='({} * dlog10per)'.format(LOG_LN_CONVERTER),
@@ -3117,7 +3181,7 @@ class Population:
                     options['ranges']['logP'][1]
                 ],
                 spacingfunc="const({}, {}, {})".format(options['ranges']['logP'][0], options['ranges']['logP'][1], options['resolutions']['logP'][0]),
-                precode="""per = 10.0**log10per
+                precode="""orbital_period = 10.0**log10per
     qmin={}/m1
     qmax=self.maximum_mass_ratio_for_RLOF(m1, per)
     """.format(options.get('Mmin', 0.07))
@@ -3126,6 +3190,7 @@ class Population:
             # binaries: mass ratio
             self.add_grid_variable(
                 name='q',
+                parameter_name='M_2',
                 longname='Mass ratio',
                 valuerange=[
                     options['ranges']['q'][0] if options.get('ranges', {}).get('q', None) else "options.get('Mmin', 0.07)/m1",
@@ -3136,8 +3201,8 @@ class Population:
                 gridtype='centred',
                 dphasevol='dq',
                 precode="""
-    m2 = q * m1
-    sep = calc_sep_from_period(m1, m2, per)
+    M_2 = q * m1
+    sep = calc_sep_from_period(M_1, M_2, orbital_period)
     """,
                 spacingfunc="const({}, {}, {})".format(
                     options['ranges']['q'][0] if options.get('ranges', {}).get('q', None) else "options.get('Mmin', 0.07)/m1",
@@ -3150,6 +3215,7 @@ class Population:
             if options['resolutions']['ecc'][0] > 0:
                 self.add_grid_variable(
                     name='ecc',
+                    parameter_name='eccentricity',
                     longname='Eccentricity',
                     resolution=options['resolutions']['ecc'][0],
                     probdist=1,
@@ -3172,10 +3238,11 @@ class Population:
                 # Triple: period
                 self.add_grid_variable(
                     name='log10per2',
+                    parameter_name='orbital_period_triple',
                     longname='log10(Orbital_Period2)',
                     resolution=options['resolutions']['logP'][1],
                     probdist=1.0,
-                    condition='self.grid_options["grid_options"]["multiplicity"] >= 3',
+                    condition='(self.grid_options["grid_options"]["multiplicity"] >= 3)',
                     branchpoint=1,
                     gridtype='centred',
                     dphasevol='({} * dlog10per2)'.format(LOG_LN_CONVERTER),
@@ -3184,9 +3251,9 @@ class Population:
                         options['ranges']['logP'][1]
                     ],
                     spacingfunc="const({}, {}, {})".format(options['ranges']['logP'][0], options['ranges']['logP'][1], options['resolutions']['logP'][1]),
-                    precode="""per2 = 10.0**log10per2
-    q2min={}/(m1+m2)
-    q2max=self.maximum_mass_ratio_for_RLOF(m1+m2, per2)
+                    precode="""orbital_period_triple = 10.0**log10per2
+    q2min={}/(M_1+M_2)
+    q2max=self.maximum_mass_ratio_for_RLOF(M_1+M_2, orbital_period_triple)
     """.format(options.get('Mmin', 0.07))
                 )
 
@@ -3195,6 +3262,7 @@ class Population:
                 print("M2 resolution {}".format(options['resolutions']['M'][2]))
                 self.add_grid_variable(
                     name='q2',
+                    parameter_name='M_3',
                     longname='Mass ratio outer/inner',
                     valuerange=[
                         options['ranges']['q'][0] if options.get('ranges', {}).get('q', None) else "options.get('Mmin', 0.07)/(m1+m2)",
@@ -3205,8 +3273,8 @@ class Population:
                     gridtype='centred',
                     dphasevol='dq2',
                     precode="""
-    m3 = q2 * (m1 + m2)
-    sep2 = calc_sep_from_period((m1+m2), m3, per2)
+    M_3 = q2 * (M_1 + M_2)
+    sep2 = calc_sep_from_period((M_1+M_2), M_3, orbital_period_triple)
     eccentricity2=0
     """,
                     spacingfunc="const({}, {}, {})".format(
@@ -3220,6 +3288,7 @@ class Population:
                 if options['resolutions']['ecc'][1] > 0:
                     self.add_grid_variable(
                         name='ecc2',
+                        parameter_name='eccentricity2',
                         longname='Eccentricity of the triple',
                         resolution=options['resolutions']['ecc'][1],
                         probdist=1,
@@ -3241,10 +3310,11 @@ class Population:
                     # Quadruple: period
                     self.add_grid_variable(
                         name='log10per3',
+                        parameter_name='orbital_period_quadruple',
                         longname='log10(Orbital_Period3)',
                         resolution=options['resolutions']['logP'][2],
                         probdist=1.0,
-                        condition='self.grid_options["grid_options"]["multiplicity"] >= 4',
+                        condition='(self.grid_options["grid_options"]["multiplicity"] >= 4)',
                         branchpoint=1,
                         gridtype='centred',
                         dphasevol='({} * dlog10per3)'.format(LOG_LN_CONVERTER),
@@ -3257,9 +3327,9 @@ class Population:
                             options['ranges']['logP'][1], 
                             options['resolutions']['logP'][2]
                         ),
-                        precode="""per3 = 10.0**log10per3
-    q3min={}/(m3)
-    q3max=self.maximum_mass_ratio_for_RLOF(m3, per3)
+                        precode="""orbital_period_quadruple = 10.0**log10per3
+    q3min={}/(M_3)
+    q3max=self.maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
     """.format(options.get('Mmin', 0.07))
                     )
 
@@ -3267,6 +3337,7 @@ class Population:
                     print("M3 resolution {}".format(options['resolutions']['M'][3]))
                     self.add_grid_variable(
                         name='q3',
+                        parameter_name='M_4',
                         longname='Mass ratio outer low/outer high',
                         valuerange=[
                             options['ranges']['q'][0] if options.get('ranges', {}).get('q', None) else "options.get('Mmin', 0.07)/(m3)",
@@ -3277,8 +3348,8 @@ class Population:
                         gridtype='centred',
                         dphasevol='dq3',
                         precode="""
-    m4 = q3 * m3
-    sep3 = calc_sep_from_period((m3), m4, per3)
+    M_4 = q3 * M_3
+    sep3 = calc_sep_from_period((M_3), M_4, orbital_period_quadruple)
     eccentricity3=0
         """,
                         spacingfunc="const({}, {}, {})".format(
@@ -3293,6 +3364,7 @@ class Population:
                     if options['resolutions']['ecc'][2] > 0:
                         self.add_grid_variable(
                             name='ecc3',
+                            parameter_name='eccentricity3',
                             longname='Eccentricity of the triple+quadruple/outer binary',
                             resolution=options['resolutions']['ecc'][2],
                             probdist=1,
@@ -3315,28 +3387,22 @@ class Population:
         # This will then be passed to the Moe_de_Stefano_2017_pdf to calculate the real probability
         # The trick we use is to strip the options_dict as a string and add some keys to it: 
 
-        # and finally the probability calculator
-        self.last_grid_variable()['probdist']="""
-    self.Moe_de_Stefano_2017_pdf(\{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}\})
-
-        """.format(str(options)[1:-1], 
+        probdist_addition = "Moe_de_Stefano_2017_pdf({{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}})".format(str(options)[1:-1], 
             '"multiplicity": multiplicity',
 
-            '"M1": m1',
-            '"M2": m2',
-            '"M3": m3',
-            '"M4": m4',
+            '"M1": M_1',
+            '"M2": M_2',
+            '"M3": M_3',
+            '"M4": M_4',
 
-            '"P": per',
-            '"P2": per2',
-            '"P3": per3',
+            '"P": orbital_period',
+            '"P2": orbital_period_triple',
+            '"P3": orbital_period_quadruple',
 
             '"ecc": eccentricity',
             '"ecc2": eccentricity2',
-            '"ecc3": eccentricity3',
+            '"ecc3": eccentricity3'
         )
 
-        print(json.dumps(json.load(self.grid_options), indent=4))
-
-# Moe_de_Stefano_2017(options={'file': '/home/david/projects/binary_c_root/binary_c/data/MdS_data.json.v2', 'testing': '1', 'resolutions': {'M': [1]}})
-
+        # and finally the probability calculator
+        # self.grid_options['_grid_variables'][self.last_grid_variable()]['probdist'] = probdist_addition