diff --git a/fsl_mrs_mce/mc_gen.py b/fsl_mrs_mce/mc_gen.py
index c86f1ef..a739a8e 100644
--- a/fsl_mrs_mce/mc_gen.py
+++ b/fsl_mrs_mce/mc_gen.py
@@ -1,10 +1,13 @@
-#  Synthetise new data and analyse it directly to return fit parameters
+# mc_gen.py - Generating data for Monte-Carlo style simulations
+#
+# Author: Konstantin E Bosbach <konstantin.bosbach@mars.uni-freiburg.de>
 
 import numpy as np
 import pandas as pd
+import time
 
 from fsl_mrs.core import MRS
-from fsl_mrs.utils import plotting, mrs_io
+from fsl_mrs.utils import mrs_io
 from fsl_mrs.utils.synthetic import synthetic_from_basis as synth
 from fsl_mrs.utils.misc import parse_metab_groups
 from fsl_mrs.utils.fitting import fit_FSLModel
@@ -31,18 +34,20 @@ def synth_and_ana(noise_cov,
     coilphase = [synth_parameter["Phi0"]]
 
     # Generate synthetic data
-    fidS, headerS, concentrationsS = synth.syntheticFromBasisFile(basis_path,
-                                                                  concentrations=synth_parameter,
-                                                                  ignore=['Gly'], ind_scaling=['mm'],
-                                                                  metab_groups='mm', broadening=broadening,  shifting=shifting,
-                                                                  # correct for complex noise
-                                                                  noisecovariance=np.divide(
-                                                                      noise_cov, 2),
-                                                                  # CAVE: baseline chosen manually
-                                                                  baseline=baseline, baseline_ppm=(
-                                                                      .2, 4.2),
-                                                                  coilphase=coilphase,
-                                                                  bandwidth=6000)
+    fidS, headerS, concentrationsS = synth.syntheticFromBasisFile(
+        basis_path,
+        concentrations=synth_parameter,
+        ignore=['Gly'], ind_scaling=['mm'],
+        metab_groups='mm', broadening=broadening,
+        shifting=shifting,
+        # correct for complex noise
+        noisecovariance=np.divide(noise_cov, 2),
+        # CAVE: baseline chosen manually
+        baseline=baseline, baseline_ppm=(.2, 4.2),
+        coilphase=coilphase,
+        bandwidth=6000
+    )
+
     # Create mrs object for further use
     mrsA = MRS(FID=fidS, header=headerS, basis=basis,
                basis_hdr=Bheader[0], names=names)
@@ -88,3 +93,42 @@ def synth_and_ana(noise_cov,
     df_params['noise_var'] = fit_varnoise
 
     return df_params
+
+
+def mc(
+    n, noise_sd, df_parameter_synth, basis_path, output_path,
+    fit_parameters=[], fit_snr=[], fit_sdnoise=[], fit_varnoise=[]
+):
+    """Function for calling synth_and_ana repeatedly,
+     as in Monte-Carlo approach"""
+    runtime = time.time()
+
+    # Define workspace output path
+    file_out_path = str(
+        output_path + "noise_sd_" +
+        str(round(noise_sd, 3)) + "_runs_"+str(n) + ".csv"
+    )
+
+    print(
+        "Starting  noise_sd", round(noise_sd, 2), " with ",
+        round(n, 2), "repetitions"
+    )
+
+    # Call function generation the desired amount of times
+    for k in range(0, n):
+        noise_fit = synth_and_ana(
+            [[np.square(noise_sd)]],
+            fit_parameters=fit_parameters,
+            fit_snr=fit_snr, fit_sdnoise=fit_sdnoise,
+            fit_varnoise=fit_varnoise,
+            df_parameter_synth=df_parameter_synth, basis_path=basis_path
+        )
+
+    # Print results to file
+    noise_fit.to_csv(file_out_path)
+
+    print(str("Finishing noise_sd", round(noise_sd, 2), " with ",
+              str(round(n, 2)), "repetitions, Runtime took ",
+              round(time.time()-runtime, 2), "[s]"))
+
+    return noise_fit, file_out_path