diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py
index 28e1456..fe4616a 100644
--- a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py
+++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py
@@ -171,9 +171,9 @@ show()
 #! Test Gaussianity of a stochastic process
 #!------------------------------------------
 #! TESTGAUSSIAN simulates  e(g(u)-u) = int (g(u)-u)^2 du  for Gaussian processes 
-#!  given the spectral density, S. The result is plotted if test0 is given.
-#!  This is useful for testing if the process X(t) is Gaussian.
-#!  If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.
+#! given the spectral density, S. The result is plotted if test0 is given.
+#! This is useful for testing if the process X(t) is Gaussian.
+#! If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.
 #! 
 #! As we see from the figure below: none of the simulated values of test1 is 
 #! above 1.00. Thus the data significantly departs from a Gaussian distribution. 
@@ -182,7 +182,7 @@ test0 = glc.dist2gauss()
 #! the following test takes time
 N = len(xx)
 test1 = S1.testgaussian(ns=N, cases=50, test0=test0)
-is_gaussian = sum(test1 > test0) < 5 
+is_gaussian = sum(test1 > test0) > 5 
 print(is_gaussian)
 show()
 
@@ -295,49 +295,31 @@ show()
 ##!wafostamp('','(CR)')
 #disp('Block = 24'),pause(pstate)
 # 
-##!#! Estimated spectrum compared to Torsethaugen spectrum
-#clf
-#Tp = 1.1;
-#H0 = 4*sqrt(spec2mom(S1,1))
-#St = torsethaugen([0:0.01:5],[H0  2*pi/Tp]);
-#plotspec(S1)
-#hold on
-#plotspec(St,'-.')
-#axis([0 6 0 0.4])
-##!wafostamp('','(ER)')
-#disp('Block = 25'),pause(pstate)
-#
-##!#!
-#clf
-#Snorm = St;
-#Snorm.S = Snorm.S/sa^2;
-#dt = spec2dt(Snorm)
-#disp('Block = 26'),pause(pstate)
-#
-##!#!
-#clf
-#[Sk Su] = spec2skew(St);
-#sa = sqrt(spec2mom(St,1));
-#gh = hermitetr([],[sa sk ku me]);
-#Snorm.tr = gh;
-#disp('Block = 27'),pause(pstate)
-#
-##!#! Transformed Gaussian model compared to Gaussian model
-#clf
-#dt = 0.5;
-#ysim_t = spec2sdat(Snorm,240,dt);
-#xsim_t = dat2gaus(ysim_t,Snorm.tr);
-#disp('Block = 28'),pause(pstate)
-#
-##!#! Compare
-##! In order to compare the Gaussian and non-Gaussian models we need to scale  
-##! \verb+xsim_t+ #!{\tt xsim$_t$} 
-##! to have the same first spectral moment as 
-##! \verb+ysim_t+,  #!{\tt ysim$_t$},  Since the process xsim_t has variance one
-##! which will be done by the following commands. 
-#clf
-#xsim_t(:,2) = sa*xsim_t(:,2);
-#waveplot(xsim_t,ysim_t,5,1,sa,4.5,'r.','b')
-##!wafostamp('','(CR)')
-#disp('Block = 29, Last block'),pause(pstate)
+#! Estimated spectrum compared to Torsethaugen spectrum
+#!-------------------------------------------------------
+ 
+clf()
+fp = 1.1;dw = 0.01
+H0 = S1.characteristic('Hm0')[0]
+St = wsm.Torsethaugen(Hm0=H0,Tp=2*pi/fp).tospecdata(np.arange(0,5+dw/2,dw))  
+S1.plot()
+St.plot('-.')
+axis([0, 6, 0, 0.4])
+show()
+ 
 
+#! Transformed Gaussian model compared to Gaussian model
+#!--------------------------------------------------------
+dt = St.sampling_period()
+va, sk, ku = St.stats_nl(moments='vsk' )
+#sa = sqrt(va)
+gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku, ysigma=sa)
+ 
+ysim_t = St.sim(ns=240, dt=0.5)
+xsim_t = ysim_t.copy()
+xsim_t[:,1] = gh.gauss2dat(ysim_t[:,1])
+
+ts_y = wo.mat2timeseries(ysim_t)
+ts_x = wo.mat2timeseries(xsim_t)
+ts_y.plot_wave(sym1='r.', ts=ts_x, sym2='b', sigma=sa, nsub=5, nfig=1)
+show()
diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py
index 930ec3f..ec94b60 100644
--- a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py
+++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py
@@ -45,49 +45,51 @@ show()
 #! Extreme waves - model check: the highest and steepest wave
 #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 clf()
-method = 0;
-rate = 8;
-[S, H, Ac, At, Tcf, Tcb, z_ind, yn] = ...
-   dat2steep(xx,rate,method);
-disp('Block = 4'), pause(pstate)
-
-clf
-[Smax indS]=max(S)
-[Amax indA]=max(Ac)
-spwaveplot(yn,[indA indS],'k.')
-wafostamp([],'(ER)')
-disp('Block = 5'), pause(pstate)
-
-#!#! Does the highest wave contradict a transformed Gaussian model?
-clf
-inds1 = (5965:5974)'; #! points to remove
-Nsim = 10;
-[y1, grec1, g2, test, tobs, mu1o, mu1oStd] = ...
-   reconstruct(xx,inds1,Nsim);
-spwaveplot(y1,indA-10)
-hold on
-plot(xx(inds1,1),xx(inds1,2),'+')
-lamb = 2.;
-muLstd = tranproc(mu1o-lamb*mu1oStd,fliplr(grec1));
-muUstd = tranproc(mu1o+lamb*mu1oStd,fliplr(grec1));
-plot (y1(inds1,1), [muLstd muUstd],'b-')
-axis([1482 1498 -1 3]), 
-wafostamp([],'(ER)')
-disp('Block = 6'), 
-pause(pstate)
-
-#!#! Expected value (solid) compared to data removed
-clf
-plot(xx(inds1,1),xx(inds1,2),'+'), hold on
-mu = tranproc(mu1o,fliplr(grec1));
-plot(y1(inds1,1), mu), hold off
-disp('Block = 7'), pause(pstate)
-
-#!#! Crest height PDF
+S, H = ts.wave_height_steepness(method=0)
+indS = S.argmax()
+indH = H.argmax()
+ts.plot_sp_wave([indH, indS],'k.')
+show()
+
+#! Does the highest wave contradict a transformed Gaussian model?
+#!----------------------------------------------------------------
+
+# TODO: Fix this
+
+#clf
+#inds1 = (5965:5974)'; #! points to remove
+#Nsim = 10;
+#[y1, grec1, g2, test, tobs, mu1o, mu1oStd] = ...
+#   reconstruct(xx,inds1,Nsim);
+#spwaveplot(y1,indA-10)
+#hold on
+#plot(xx(inds1,1),xx(inds1,2),'+')
+#lamb = 2.;
+#muLstd = tranproc(mu1o-lamb*mu1oStd,fliplr(grec1));
+#muUstd = tranproc(mu1o+lamb*mu1oStd,fliplr(grec1));
+#plot (y1(inds1,1), [muLstd muUstd],'b-')
+#axis([1482 1498 -1 3]), 
+#wafostamp([],'(ER)')
+#disp('Block = 6'), 
+#pause(pstate)
+#
+##!#! Expected value (solid) compared to data removed
+#clf
+#plot(xx(inds1,1),xx(inds1,2),'+'), hold on
+#mu = tranproc(mu1o,fliplr(grec1));
+#plot(y1(inds1,1), mu), hold off
+#disp('Block = 7'), pause(pstate)
+
+#! Crest height PDF
+#!------------------
 #! Transform data so that kde works better
-clf
+clf()
+wave_data = ts.wave_parameters()
+Ac = wave_data['Ac']
 L2 = 0.6;
-plotnorm(Ac.^L2)
+import pylab
+ws.probplot(Ac**L2, dist='norm', plot=pylab)
+show()
 
 #!#!
 #!
diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py
index 1209b03..a4017f9 100644
--- a/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py
+++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py
@@ -65,10 +65,10 @@ mM_rfc = tp.cycle_pairs(h=0.3)
 clf()
 subplot(122), 
 mM.plot() 
-title('min-max cycle count')
+title('min-max cycle pairs')
 subplot(121), 
 mM_rfc.plot()
-title('Rainflow cycle count')
+title('Rainflow filtered cycles')
 show()
 
 #! Min-max and rainflow cycle distributions
diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py
index 9f63c66..5813629 100644
--- a/pywafo/src/wafo/objects.py
+++ b/pywafo/src/wafo/objects.py
@@ -63,7 +63,7 @@ class LevelCrossings(WafoData):
     Member variables
     ----------------
     data : array-like
-        number of upcrossings
+        number of upcrossings or upcrossingintensity
     args : array-like
         crossing levels
         
@@ -1858,7 +1858,7 @@ class TimeSeries(WafoData):
         plot = plotbackend.plot
         subplot = plotbackend.subplot
         figs = []
-        for iz in xrange(nfig):
+        for unused_iz in xrange(nfig):
             figs.append(plotbackend.figure())
             plotbackend.title('Surface elevation from mean water level (MWL).')
             for ix in xrange(nsub):
@@ -1883,10 +1883,12 @@ class TimeSeries(WafoData):
         return figs
         
 
-    def plot_sp_wave(self, wave_idx_, tz_idx=None, *args, **kwds):
+    def plot_sp_wave(self, wave_idx_, *args, **kwds):
         """
         Plot specified wave(s) from timeseries
-
+        
+        Parameters
+        ----------
         wave_idx : integer vector
             of indices to waves we want to plot, i.e., wave numbers.
         tz_idx : integer vector
@@ -1907,8 +1909,9 @@ class TimeSeries(WafoData):
         plot_wave, findtc
         """
         wave_idx = atleast_1d(wave_idx_).flatten()
+        tz_idx = kwds.pop('tz_idx', None)
         if tz_idx is None:
-            tc_ind, tz_idx = findtc(self.data, 0, 'tw') # finding trough to trough waves
+            unused_tc_ind, tz_idx = findtc(self.data, 0, 'tw') # finding trough to trough waves
 
         dw = nonzero(abs(diff(wave_idx)) > 1)[0]
         Nsub = dw.size + 1
@@ -1932,7 +1935,7 @@ class TimeSeries(WafoData):
         Nfig = int(ceil(Nsub / 6))
         Nsub = min(6, int(ceil(Nsub / Nfig)))
         figs = []
-        for iy in range(Nfig):
+        for unused_iy in range(Nfig):
             figs.append(plotbackend.figure())
             for ix in range(Nsub):
                 plotbackend.subplot(Nsub, 1, mod(ix, Nsub) + 1)
diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py
index d05a98b..abebfa6 100644
--- a/pywafo/src/wafo/transform/core.py
+++ b/pywafo/src/wafo/transform/core.py
@@ -166,19 +166,17 @@ class TrData(WafoData, TrCommon):
     >>> g.dist2gauss() < 1e-16
     True
     """
-    def __init__(self, *args, **kwds):
-        self.ymean = kwds.pop('ymean', 0e0)
-        self.ysigma = kwds.pop('ysigma', 1e0)
-        self.mean = kwds.pop('mean', None)
-        self.sigma = kwds.pop('sigma', None)
-        
+    def __init__(self, *args, **kwds): 
         options = dict(title='Transform',
                             xlab='x', ylab='g(x)',
                             plot_args=['r'],
                             plot_args_children=['g--'],)
         options.update(**kwds)
         super(TrData, self).__init__(*args, **options)
-        
+        self.ymean = kwds.get('ymean', 0e0)
+        self.ysigma = kwds.get('ysigma', 1e0)
+        self.mean = kwds.get('mean', None)
+        self.sigma = kwds.get('sigma', None)
         
         if self.mean is None: 
             #self.mean = np.mean(self.args) #