diff --git a/mngseg.py b/mngseg.py
index 325c317..b19270f 100644
--- a/mngseg.py
+++ b/mngseg.py
@@ -1,4 +1,4 @@
-from pywavan import powspec, nb_scale
+from pywavan import powspec, fan_trans, nb_scale
 from astropy.io import fits
 import numpy as np
 from matplotlib import pyplot as plt
@@ -236,7 +236,7 @@ def load_results(label):
     return wav_k, S1a, wt, S11a, q
 
 
-def analyse_sim(im):
+def analyse_sim(im, load=False):
     """ Do the MnGseg analysis """
     meanim = np.mean(im)
     imzm = im - meanim
@@ -247,16 +247,19 @@ def analyse_sim(im):
     label = "final"  # label to identify parameter setup
 
     # after a lot of trials I found a fixed q=2 is a good value
-    q = [2.0] * nb_scale(imzm.shape)
-    # wt, S11a, wav_k, S1a, q = fan_trans(imzm, reso=1, q=q, qdyn=False)
-    # alternatively you can let pywavan determine it automatically by setting skewl
-    # q=[3.0]*nb_scale(imzm.shape)
-    # wt, S11a, wav_k, S1a, q = fan_trans(imzm, reso=1, q=q, qdyn=True, skewl=0.4)
-    # print(q)
+    if not load:
+        q = [2.0] * nb_scale(imzm.shape)
+        wt, S11a, wav_k, S1a, q = fan_trans(imzm, reso=1, q=q, qdyn=False)
+        # alternatively you can let pywavan determine it automatically by setting skewl
+        q = [3.0] * nb_scale(imzm.shape)
+        wt, S11a, wav_k, S1a, q = fan_trans(imzm, reso=1, q=q, qdyn=True, skewl=0.4)
+        print(q)
 
-    # save results because it can be long to calculate (especially if qdyn=True). Remark that wt and S11a are quite big
-    # save_results(wt, S11a, wav_k, S1a, q, label)
-    wav_k, S1a, wt, S11a, q = load_results(label)
+        # save results because it can be long to calculate (especially if qdyn=True).
+        # Remark that wt and S11a are quite big
+        save_results(wt, S11a, wav_k, S1a, q, label)
+    if load:
+        wav_k, S1a, wt, S11a, q = load_results(label)
 
     # make images of the Gaussian and coherent part
     make_images(im, wt, M, meanim, label)