diff --git a/pywafo/src/wafo/covariance/core.py b/pywafo/src/wafo/covariance/core.py
index a19a7a2..9337937 100644
--- a/pywafo/src/wafo/covariance/core.py
+++ b/pywafo/src/wafo/covariance/core.py
@@ -20,7 +20,7 @@ import warnings
 from numpy import (zeros, sqrt, dot, newaxis, inf, where, pi, nan, #@UnresolvedImport
                    atleast_1d, hstack, vstack, r_, linspace, flatnonzero, size, #@UnresolvedImport
                    isnan, finfo, diag, ceil, floor, random, pi) #@UnresolvedImport
-from numpy.fft import fft
+from numpy.fft import fft #as fft
 from numpy.random import randn
 import scipy.interpolate as interpolate
 from scipy.linalg import toeplitz, sqrtm, svd, cholesky, diagsvd, pinv
@@ -182,7 +182,7 @@ class CovData1D(WafoData):
 ##        wdata = CovData1D(**kwds)
 ##        return wdata
 
-    def tospecdata(self, rate=None, method='linear', nugget=0.0, trunc=1e-5, fast=True):
+    def tospecdata(self, rate=None, method='fft', nugget=0.0, trunc=1e-5, fast=True):
         '''
         Computes spectral density from the auto covariance function
 
@@ -192,7 +192,7 @@ class CovData1D(WafoData):
             1,2,4,8...2^r, interpolation rate for f (default 1)
 
         method: string
-            interpolation method 'stineman', 'linear', 'cubic'
+            interpolation method 'stineman', 'linear', 'cubic', 'fft'
 
         nugget = scalar, real
             nugget effect to ensure that round off errors do not result in
@@ -217,7 +217,8 @@ class CovData1D(WafoData):
         Example:
         >>> import wafo.spectrum.models as sm
         >>> import numpy as np
-        >>> import scipy.signal.signaltools as st
+        >>> import scipy.signal as st
+        >>> import pylab
         >>> L = 129
         >>> t = np.linspace(0,75,L)
         >>> R = np.zeros(L)
@@ -230,7 +231,13 @@ class CovData1D(WafoData):
         >>> S = Sj.tospecdata()
         >>> R2 = S.tocovdata()
         >>> S1 = R2.tospecdata()
-        >>> assert(all(abs(S1.data-S.data)<1e-4) ,'COV2SPEC')
+        >>> abs(S1.data-S.data).max()
+        
+        >>> S1.plot('r-')
+        >>> S.plot('b:')
+        >>> pylab.show()
+        
+        >>> all(abs(S1.data-S.data)<1e-4)
 
         See also
         --------
@@ -238,10 +245,10 @@ class CovData1D(WafoData):
         datastructures
         '''
 
-        dT = self.sampling_period()
-        # dT = time-step between data points.
+        dt = self.sampling_period()
+        # dt = time-step between data points.
 
-        ACF, unused_ti = atleast_1d(self.data, self.args)
+        acf, unused_ti = atleast_1d(self.data, self.args)
 
         if self.lagtype in 't':
             spectype = 'freq'
@@ -258,30 +265,34 @@ class CovData1D(WafoData):
 
         ## add a nugget effect to ensure that round off errors
         ## do not result in negative spectral estimates
-        ACF[0] = ACF[0] + nugget
-        n = ACF.size
+        acf[0] = acf[0] + nugget
+        n = acf.size
         # embedding a circulant vector and Fourier transform
-        if fast:
-            nfft = 2 ** nextpow2(2 * n - 2)
-        else:
-            nfft = 2 * n - 2
-
+        
+        nfft = 2 ** nextpow2(2 * n - 2) if fast else 2 * n - 2
+        
+        if method=='fft':
+            nfft *= rate
+         
         nf = nfft / 2 ## number of frequencies
-        ACF = r_[ACF, zeros(nfft - 2 * n + 2), ACF[n - 1:0:-1]]
-
-        Rper = (fft(ACF, nfft).real).clip(0) ## periodogram
+        acf = r_[acf, zeros(nfft - 2 * n + 2), acf[n - 2:0:-1]]
+        
+        Rper = (fft(acf, nfft).real).clip(0) ## periodogram
+#        import pylab
+#        pylab.semilogy(Rper)
+#        pylab.show()
         RperMax = Rper.max()
         Rper = where(Rper < trunc * RperMax, 0, Rper)
     
-        S = abs(Rper[0:(nf + 1)]) * dT / pi
-        w = linspace(0, pi / dT, nf + 1)
+        S = abs(Rper[0:(nf + 1)]) * dt / pi
+        w = linspace(0, pi / dt, nf + 1)
         So = _wafospec.SpecData1D(S, w, type=spectype, freqtype=ftype)
         So.tr = self.tr
         So.h = self.h
         So.norm = self.norm
 
-        if rate > 1:
-            So.args = linspace(0, pi / dT, nf * rate)
+        if method != 'fft' and rate > 1:
+            So.args = linspace(0, pi / dt, nf * rate)
             if method == 'stineman':
                 So.data = stineman_interp(So.args, w, S)
             else:
@@ -374,15 +385,15 @@ class CovData1D(WafoData):
 
         _set_seed(iseed)
 
-        ACF = self.data.ravel()
-        n = ACF.size
+        acf = self.data.ravel()
+        n = acf.size
 
 
-        I = ACF.argmax()
+        I = acf.argmax()
         if I != 0:
             raise ValueError('ACF does not have a maximum at zero lag')
 
-        ACF.shape = (n, 1)
+        acf.shape = (n, 1)
 
         dT = self.sampling_period()
 
@@ -393,25 +404,24 @@ class CovData1D(WafoData):
 
         ## add a nugget effect to ensure that round off errors
         ## do not result in negative spectral estimates
-        ACF[0] = ACF[0] + nugget
+        acf[0] = acf[0] + nugget
 
         ## Fast and exact simulation of simulation of stationary
         ## Gaussian process throug circulant embedding of the
         ## Covariance matrix
         floatinfo = finfo(float)
-        if (abs(ACF[-1]) > floatinfo.eps): ## assuming ACF(n+1)==0
+        if (abs(acf[-1]) > floatinfo.eps): ## assuming acf(n+1)==0
             m2 = 2 * n - 1
             nfft = 2 ** nextpow2(max(m2, 2 * ns))
-            ACF = r_[ACF, zeros((nfft - m2, 1)), ACF[-1:0:-1, :]]
-            #disp('Warning: I am now assuming that ACF(k)=0 ')
-            #disp('for k>MAXLAG.')
+            acf = r_[acf, zeros((nfft - m2, 1)), acf[-1:0:-1, :]]
+            #warnings,warn('I am now assuming that ACF(k)=0 for k>MAXLAG.')
         else: # # ACF(n)==0
             m2 = 2 * n - 2
             nfft = 2 ** nextpow2(max(m2, 2 * ns))
-            ACF = r_[ACF, zeros((nfft - m2, 1)), ACF[n - 1:1:-1, :]]
+            acf = r_[acf, zeros((nfft - m2, 1)), acf[n - 1:1:-1, :]]
 
         ##m2=2*n-2
-        S = fft(ACF, nfft, axis=0).real ## periodogram
+        S = fft(acf, nfft, axis=0).real ## periodogram
 
         I = S.argmax()
         k = flatnonzero(S < 0)
diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py
index afd2d19..5d09b9a 100644
--- a/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py
+++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py
@@ -1,6 +1,28 @@
 from scipy import *
 from pylab import *
 
+#import wafo.spectrum.models as sm
+#Sj = sm.Jonswap()
+#S = Sj.tospecdata()
+#S.data[0:40] = 0.0
+#S.data[100:-1] = 0.0
+#Nt = len(S.data)-1
+#acf = S.tocovdata(nr=0, nt=Nt)
+#S2 = acf.tospecdata()
+#S.plot('r')
+#S2.plot('b:')
+#
+#show()
+#        
+#import wafo
+#import wafo.objects as wo
+#xn = wafo.data.sea()
+#ts = wo.mat2timeseries(xn)
+#Sest = ts.tospecdata(method='cov')
+#Sest.setplotter('semilogy')
+#Sest.plot()
+#show()
+
 # pyreport -o chapter1.html chapter1.py
 
 #! CHAPTER1 demonstrates some applications of WAFO
@@ -26,6 +48,7 @@ S1 = S.tospecdata()
 S1.plot()
 show()
 
+
 ##
 import wafo.objects as wo
 xs = S1.sim(ns=2000, dt=0.1)
@@ -33,6 +56,7 @@ ts = wo.mat2timeseries(xs)
 ts.plot_wave('-')
 show()
 
+
 #! Estimation of spectrum 
 #!~~~~~~~~~~~~~~~~~~~~~~~
 #! A common situation is that one wants to estimate the spectrum for wave
@@ -42,19 +66,20 @@ clf()
 Fs = 4;  
 xs = S1.sim(ns=fix(20 * 60 * Fs), dt=1. / Fs) 
 ts = wo.mat2timeseries(xs) 
-Sest = ts.tospecdata(NFFT=400)
+Sest = ts.tospecdata(L=400)
 S1.plot()
 Sest.plot('--')
 axis([0, 3, 0, 5]) # This may depend on the simulation
 show()
 
-## Section 1.4.2 Probability distributions of wave characteristics.
-## Probability distribution of wave trough period
-# WAFO gives the possibility of computing the exact probability
-# distributions for a number of characteristics given a spectral density.
-# In the following example we study the trough period extracted from the
-# time series and compared with the theoretical density computed with exact
-# spectrum, S1, and the estimated spectrum, Sest.
+#! Section 1.4.2 Probability distributions of wave characteristics.
+#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#! Probability distribution of wave trough period:
+#! WAFO gives the possibility of computing the exact probability
+#! distributions for a number of characteristics given a spectral density.
+#! In the following example we study the trough period extracted from the
+#! time series and compared with the theoretical density computed with exact
+#! spectrum, S1, and the estimated spectrum, Sest.
 clf()
 import wafo.misc as wm
 dtyex = S1.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)
@@ -81,11 +106,11 @@ Sp = 15;   # spreading parameter
 D1 = sm.Spreading(type='cos', theta0=th0, method=None) # frequency independent
 D12 = sm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent
 
-#SD1 = mkdspec(S1, D1)
-#SD12 = mkdspec(S1, D12);
-#plotspec(SD1, plotflag), hold on, plotspec(SD12, plotflag, '-.'); hold off
-#wafostamp('', '(ER)')
-#disp('Block = 5'), pause(pstate)
+SD1 = D1.tospecdata2d(S1)
+SD12 = D12.tospecdata2d(S1)
+SD1.plot()
+SD12.plot()#linestyle='dashdot')
+show()
 
 #! 3D Simulation of the sea surface 
 #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py
index ef4291b..e9d4892 100644
--- a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py
+++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py
@@ -180,6 +180,7 @@ glc, gemp = lc.trdata(mean=me, sigma=sa)
 g.plot('r')
 glc.plot('b-') #! Transf. estimated from level-crossings
 gh.plot('b-.')  #! Hermite Transf. estimated from moments
+grid('on')
 show()
  
 #!  Test Gaussianity of a stochastic process.
@@ -205,7 +206,7 @@ show()
 #! "light" lower tail. 
 clf()
 import pylab
-ws.probplot(ts.data, dist='norm', plot=pylab)
+ws.probplot(ts.data.ravel(), dist='norm', plot=pylab)
 show()
 #! Section 2.2.3 Spectral densities of sea data
 #!-----------------------------------------------
diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py
index c467ebe..48be7a4 100644
--- a/pywafo/src/wafo/misc.py
+++ b/pywafo/src/wafo/misc.py
@@ -2157,4 +2157,4 @@ if __name__ == "__main__":
         import doctest
         doctest.testmod()
     else:
-        test_find_cross()
+        _test_find_cross()
diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py
index e9e9cf1..2ed9622 100644
--- a/pywafo/src/wafo/objects.py
+++ b/pywafo/src/wafo/objects.py
@@ -197,7 +197,7 @@ class LevelCrossings(WafoData):
         integral = zeros(u.shape, dtype=float)
         for i in range(len(integral)):
             y = factor1 * exp(-u[i] * u[i] * factor2)
-            integral[i] = trapz(x, y)
+            integral[i] = trapz(y, x)
         #end
         G = G - integral / (2 * pi)
         G = G / max(G)
@@ -215,7 +215,7 @@ class LevelCrossings(WafoData):
         
         g = lc2.trdata()
 
-        f = [u, u]
+        #f = [u, u]
         f = g.dat2gauss(Z)
         G = TrData(f, u)
         
@@ -320,11 +320,11 @@ class LevelCrossings(WafoData):
         >>> g1, g1emp = lc.trdata(gvar=0.5 ) # Equal weight on all points
         >>> g2, g2emp = lc.trdata(gvar=[3.5, 0.5, 3.5])  # Less weight on the ends
         >>> int(S.tr.dist2gauss()*100)
-        593
+        141
         >>> int(g0emp.dist2gauss()*100)
-        492
+        1544
         >>> int(g0.dist2gauss()*100)
-        361
+        340
         >>> int(g1.dist2gauss()*100)
         352
         >>> int(g2.dist2gauss()*100)
@@ -855,12 +855,13 @@ class TimeSeries(WafoData):
         c0 = R[0]
         if norm:
             R = R / c0
+        r0 = R[0]
         if dt is None:
             dt = self.sampling_period()
         t = linspace(0, lag * dt, lag + 1)
         #cumsum = np.cumsum
         acf = _wafocov.CovData1D(R[lags], t)
-        acf.stdev = sqrt(r_[ 0, 1 , 1 + 2 * cumsum(R[1:] ** 2)] / Ncens)
+        acf.stdev = sqrt(r_[ 0, r0**2 , r0**2 + 2 * cumsum(R[1:] ** 2)] / Ncens)
         acf.children = [WafoData(-2. * acf.stdev[lags], t), WafoData(2. * acf.stdev[lags], t)]
         acf.plot_args_children = ['r:']
         acf.norm = norm
@@ -910,7 +911,7 @@ class TimeSeries(WafoData):
         fact = 2.0 * pi
         w = fact * f
         return _wafospec.SpecData1D(S / fact, w)
-    def tospecdata(self, L=None, tr=None, method='cov', detrend=detrend_mean, window=parzen, noverlap=0, pad_to=None, ftype='w', alpha=None):
+    def tospecdata(self, L=None, tr=None, method='cov', detrend=detrend_mean, window=parzen, noverlap=0, ftype='w', alpha=None):
         '''
         Estimate one-sided spectral density from data.
      
@@ -930,9 +931,14 @@ class TimeSeries(WafoData):
             'cov' :  Frequency smoothing using a parzen window function
                     on the estimated autocovariance function.  (default)
             'psd' : Welch's averaged periodogram method with no overlapping batches
-        dflag : string 
+        detrend : function
             defining detrending performed on the signal before estimation.
-            'mean','linear' or 'ma' (= moving average)  (default 'mean')   
+            (default detrend_mean)   
+        window : vector of length NFFT or function
+            To create window vectors see numpy.blackman, numpy.hamming,
+            numpy.bartlett, scipy.signal, scipy.signal.get_window etc.
+        noverlap : scalar int
+             gives the length of the overlap between segments.
         ftype : character
             defining frequency type: 'w' or 'f'  (default 'w')    
      
@@ -968,9 +974,9 @@ class TimeSeries(WafoData):
         
         #% Initialize constants 
         #%~~~~~~~~~~~~~~~~~~~~~
-        nugget = 0; #%10^-12;
+        nugget = 1e-12
         rate = 2; #% interpolationrate for frequency
-        tapery = 0; #% taper the data before the analysis
+        
         wdef = 1; #% 1=parzen window 2=hanning window, 3= bartlett window
           
         dt = self.sampling_period()
@@ -981,18 +987,17 @@ class TimeSeries(WafoData):
         L = min(L, n);
          
         max_L = min(300, n); #% maximum lag if L is undetermined
-        change_L = L is None
-        if change_L:
+        estimate_L = L is None
+        if estimate_L:
             L = min(n - 2, int(4. / 3 * max_L + 0.5))
-
             
-        if method == 'cov' or change_L:
+        if method == 'cov' or estimate_L:
             tsy = TimeSeries(yy, self.args)
             R = tsy.tocovdata() 
-            if  change_L:
+            if  estimate_L:
                 #finding where ACF is less than 2 st. deviations.
-                L = max_L - (np.abs(R.data[max_L::-1]) > 2 * R.stdev[max_L::-1]).argmax() # a better L value  
-                if wdef == 1:  # % modify L so that hanning and Parzen give appr. the same result
+                L = max_L + 2 - (np.abs(R.data[max_L::-1]) > 2 * R.stdev[max_L::-1]).argmax() # a better L value  
+                if wdef == 1:  # modify L so that hanning and Parzen give appr. the same result
                     L = min(int(4 * L / 3), n - 2)
                 print('The default L is set to %d' % L)       
         try:
@@ -1014,9 +1019,9 @@ class TimeSeries(WafoData):
             Be = None
               
         if method == 'psd':
-            nf = rate * 2 ** nextpow2(2 * L - 2) #  Interpolate the spectrum with rate 
-            nfft = 2 * nf 
-            S, f = psd(yy, Fs=1. / dt, NFFT=nfft, detrend=detrend, window=window,
+            nfft =  2 ** nextpow2(L) 
+            pad_to = rate*nfft #  Interpolate the spectrum with rate 
+            S, f = psd(yy, Fs=1. / dt, NFFT=nfft, detrend=detrend, window=window(nfft),
                    noverlap=noverlap, pad_to=pad_to, scale_by_freq=True)
             fact = 2.0 * pi
             w = fact * f
@@ -1025,10 +1030,13 @@ class TimeSeries(WafoData):
             # add a nugget effect to ensure that round off errors
             # do not result in negative spectral estimates
              
-            R.data = R.data[:L] * win[L - 1::] 
-            R.args = R.args[:L]
-            
-            spec = R.tospecdata(rate=2, nugget=nugget)
+            R.data[:L] = R.data[:L] * win[L - 1::]
+            R.data[L] = 0.0
+            R.data = R.data[:L+1]
+            R.args = R.args[:L+1]
+            #R.plot()
+            #R.show()
+            spec = R.tospecdata(rate=rate, nugget=nugget)
              
         spec.Bw = Be
         if ftype == 'f':
@@ -1356,7 +1364,9 @@ class TimeSeries(WafoData):
             t = self.args[ind]
         except:
             t = ind
-        return TurningPoints(self.data[ind], t)
+        mean = self.data.mean()
+        stdev = self.data.std()
+        return TurningPoints(self.data[ind], t, mean=mean, stdev=stdev)
     
     def wave_periods(self, vh=None, pdef='d2d', wdef=None, index=None, rate=1):
         """ 
diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py
index 3483d3f..523fc63 100644
--- a/pywafo/src/wafo/spectrum/core.py
+++ b/pywafo/src/wafo/spectrum/core.py
@@ -649,20 +649,17 @@ class SpecData1D(WafoData):
     """
 
     def __init__(self, *args, **kwds):
-        super(SpecData1D, self).__init__(*args, **kwds)
-        self.name = 'WAFO Spectrum Object'
-        self.type = 'freq'
-        self.freqtype = 'w'
+        self.name_ = kwds.pop('name', 'WAFO Spectrum Object') 
+        self.type = kwds.pop('type','freq')
+        self.freqtype = kwds.pop('freqtype','w')
         self.angletype = ''
-        self.h = inf
-        self.tr = None #TrLinear()
-        self.phi = 0.0
-        self.v = 0.0
-        self.norm = False
-        somekeys = ['phi', 'name', 'h', 'tr', 'freqtype', 'v','type', 'norm']
-
-        self.__dict__.update(sub_dict_select(kwds, somekeys))
-
+        self.h = kwds.pop('h',inf)
+        self.tr = kwds.pop('tr',None) #TrLinear()
+        self.phi = kwds.pop('phi',0.0)
+        self.v = kwds.pop('v',0.0)
+        self.norm = kwds.pop('norm',False)
+        super(SpecData1D, self).__init__(*args, **kwds)
+        
         self.setlabels()
 
     def tocov_matrix(self, nr=0, nt=None, dt=None):
@@ -869,7 +866,7 @@ class SpecData1D(WafoData):
         nfft = rate * 2 ** nextpow2(2 * n_f - 2)
 
         # periodogram
-        rper = r_[specn, zeros(nfft - (2 * n_f) + 2), conj(specn[n_f - 1:0:-1])] 
+        rper = r_[specn, zeros(nfft - (2 * n_f) + 2), conj(specn[n_f - 2:0:-1])] 
         time = r_[0:nt + 1] * d_t * (2 * n_f - 2) / nfft
 
         r = fft(rper, nfft).real / (2 * n_f - 2)
@@ -1194,9 +1191,6 @@ class SpecData1D(WafoData):
        
         if self.tr is None:
             g = TrLinear(var=m[0])
-            #y = linspace(-5, 5, 513)
-            #g = _wafotransform.
-            #g = TrData(y, sqrt(m[0]) * y)
         else:
             g = self.tr
         
@@ -2086,7 +2080,7 @@ class SpecData1D(WafoData):
 #        
 #        opt = troptset(opt,'multip',1)
         
-        plotflag = 0 if test0 is None else 1
+        plotflag = False if test0 is None else True
         if cases > 50:
             print('  ... be patient this may take a while')
         
@@ -2104,9 +2098,6 @@ class SpecData1D(WafoData):
                 ts = TimeSeries(xs[:, iy], xs[:, 0].ravel())
                 g, tmp = ts.trdata(method, **opt)
                 test1.append(g.dist2gauss())
-            #xs = cov2sdat(R,[ns Nstep]);
-            #[g, tmp] = dat2tr(xs,method, **opt);
-            #test1 = [test1; tmp(:)]
             if verbose:
                 print('finished %d of %d ' % (ix + 1, rep))
         
@@ -2790,6 +2781,8 @@ class SpecData2D(WafoData):
     """
 
     def __init__(self, *args, **kwds):
+    
+        
         super(SpecData2D, self).__init__(*args, **kwds)
 
         self.name = 'WAFO Spectrum Object'
diff --git a/pywafo/src/wafo/spectrum/models.py b/pywafo/src/wafo/spectrum/models.py
index bd88885..2ae535c 100644
--- a/pywafo/src/wafo/spectrum/models.py
+++ b/pywafo/src/wafo/spectrum/models.py
@@ -1952,6 +1952,8 @@ class Spreading(object):
         w = specdata.args
         S = specdata.data
         D, phi0 = self(theta, w=w, wc=wc)
+        if D.ndim != 2: # frequency dependent spreading
+            D = D[:, None]
             
         SD = D * S[None,:]
         
diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py
index d08b9ef..392856e 100644
--- a/pywafo/src/wafo/transform/core.py
+++ b/pywafo/src/wafo/transform/core.py
@@ -70,7 +70,7 @@ class TrCommon(object):
         -------
         t0 : real, scalar
             a measure of departure from the Gaussian model calculated as
-            trapz(xn,(xn-g(x))**2.) where int. limits is given by X.
+            trapz((xn-g(x))**2., xn) where int. limits is given by X.
         """
         if x is None:
             xn = linspace(xnmin, xnmax, n)
@@ -79,7 +79,7 @@ class TrCommon(object):
             xn = (x-self.mean)/self.sigma
 
         yn = (self._dat2gauss(x)-self.ymean)/self.ysigma
-        t0 = trapz(xn,(xn-yn)**2.)
+        t0 = trapz((xn-yn)**2., xn)
         return t0
 
     def gauss2dat(self, y, *yi):
@@ -166,14 +166,18 @@ class TrData(WafoData, TrCommon):
     True
     """
     def __init__(self, *args, **kwds):
-        super(TrData, self).__init__(*args, **kwds)
-        self.labels.title = 'Transform'
-        self.labels.ylab = 'g(x)'
-        self.labels.xlab = 'x'
-        self.ymean = kwds.get('ymean', 0e0)
-        self.ysigma = kwds.get('ysigma', 1e0)
-        self.mean = kwds.get('mean', None)
-        self.sigma = kwds.get('sigma', None)
+        self.ymean = kwds.pop('ymean', 0e0)
+        self.ysigma = kwds.pop('ysigma', 1e0)
+        self.mean = kwds.pop('mean', None)
+        self.sigma = kwds.pop('sigma', None)
+        
+        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)
+        
         
         if self.mean is None: 
             #self.mean = np.mean(self.args) # 
@@ -182,6 +186,8 @@ class TrData(WafoData, TrCommon):
             yp = self.ymean+self.ysigma
             ym = self.ymean-self.ysigma
             self.sigma = (self.gauss2dat(yp)-self.gauss2dat(ym))/2.
+            
+        self.children = [WafoData((self.args-self.mean)/self.sigma, self.args)]
 
     def _gauss2dat(self, y, *yi):
         return tranproc(self.data, self.args, y, *yi)
diff --git a/pywafo/src/wafo/transform/models.py b/pywafo/src/wafo/transform/models.py
index 15d0b7d..815f557 100644
--- a/pywafo/src/wafo/transform/models.py
+++ b/pywafo/src/wafo/transform/models.py
@@ -61,7 +61,7 @@ class TrCommon2(TrCommon):
         -------
         t0 : real, scalar
             a measure of departure from the Gaussian model calculated as
-            trapz(xn,(xn-g(x))**2.) where int. limits is given by X.
+            trapz((xn-g(x))**2., xn) where int. limits is given by X.
         """
         if x is None:
             xn = np.linspace(xnmin, xnmax, n)
diff --git a/pywafo/src/wafo/wafodata.py b/pywafo/src/wafo/wafodata.py
index 467554f..443d2e0 100644
--- a/pywafo/src/wafo/wafodata.py
+++ b/pywafo/src/wafo/wafodata.py
@@ -81,7 +81,7 @@ class WafoData(object):
         if self.children != None:
             plotbackend.hold('on')
             tmp = []
-            child_args = args + tuple(self.plot_args_children)
+            child_args = args if len(args) else tuple(self.plot_args_children)
             child_kwds = dict()
             child_kwds.update(self.plot_kwds_children)
             child_kwds.update(**kwds)
@@ -91,7 +91,7 @@ class WafoData(object):
                     tmp.append(tmp1)
             if len(tmp) == 0:
                 tmp = None
-        main_args = args + tuple(self.plot_args)
+        main_args = args if len(args) else tuple(self.plot_args)
         main_kwds = dict()
         main_kwds.update(self.plot_kwds)
         main_kwds.update(kwds)