diff --git a/pywafo/src/wafo/__init__.py b/pywafo/src/wafo/__init__.py
index 876ea6f..bb654e1 100644
--- a/pywafo/src/wafo/__init__.py
+++ b/pywafo/src/wafo/__init__.py
@@ -1,7 +1,8 @@
 
 from info import __doc__
 import misc
-import data
+import data
+import demos
 import kdetools
 import objects
 import spectrum
diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py
index dc40766..85d1153 100644
--- a/pywafo/src/wafo/objects.py
+++ b/pywafo/src/wafo/objects.py
@@ -179,7 +179,8 @@ class LevelCrossings(WafoData):
         L0 = vstack((L0, r1 * L0 + sqrt(1 - r2 ** 2) * randn(1)))
         #%Simulate the process, starting in L0
         lfilter = scipy.signal.lfilter
-        L = lfilter(1, [1, a1, a2], e, lfilter([1, a1, a2], 1, L0))
+        # TODO: lfilter crashes the example
+        L = lfilter(ones(1), [1, a1, a2], e, lfilter([1, a1, a2], ones(1), L0))
 
         epsilon = 1.01
         min_L = min(L)
@@ -322,20 +323,20 @@ class LevelCrossings(WafoData):
         >>> int(S.tr.dist2gauss()*100)
         141
         >>> int(g0emp.dist2gauss()*100)
-        1544
+        380995
         >>> int(g0.dist2gauss()*100)
-        340
+        143
         >>> int(g1.dist2gauss()*100)
-        352
+        162
         >>> int(g2.dist2gauss()*100)
-        365
+        120
         
-         hold on, trplot(g1,g)                          # Check the fit
-         trplot(g2)
-        
-         See also  troptset, dat2tr, trplot, findcross, smooth
+        g0.plot() # Check the fit.
+         
+        See also
+          troptset, dat2tr, trplot, findcross, smooth
         
-         NB! the transformated data will be N(0,1)
+        NB! the transformated data will be N(0,1)
         
         Reference
         --------- 
@@ -345,22 +346,14 @@ class LevelCrossings(WafoData):
         Marine structures, Design, Construction and Safety, 
         Vol 10, pp 13--47
         ''' 
-        
-        
-        # Tested on: Matlab 5.3, 5.2, 5.1
-        # History:
-        # by pab 29.12.2000
-        # based on lc2tr, but the inversion is faster.
-        # by IR and PJ
+ 
         if mean is None:
             mean = self.mean
         if sigma is None:
             sigma = self.stdev
         
         opt = DotDict(chkder=True, plotflag=False, csm=0.9, gsm=.05,
-            param=(-5, 5, 513), delay=2, linextrap=True, ntr=inf, ne=7, cvar=1, gvar=1)
-        # If just 'defaults' passed in, return the default options in g
-        
+            param=(-5, 5, 513), delay=2, linextrap=True, ntr=10000, ne=7, gvar=1)
         opt.update(options)
         param = opt.param
         Ne = opt.ne
@@ -380,21 +373,11 @@ class LevelCrossings(WafoData):
             gvar = opt.gvar * ones(ncr)
         else:
             gvar = interp1d(linspace(0, 1, ng) , opt.gvar, kind='linear')(linspace(0, 1, ncr))  
+          
         
-        ng = len(atleast_1d(opt.cvar))
-        if ng == 1:
-            cvar = opt.cvar * ones(ncr)
-        else:
-            cvar = interp1d(linspace(0, 1, ng), opt.cvar, kind='linear')(linspace(0, 1, ncr))  
-        
-        
-       
         uu = linspace(*param)
-        
         g1 = sigma * uu + mean
         
-        g22 = lc2.copy() 
-        
         if Ner > 0: # Compute correction factors
             cor1 = trapz(lc2[0:Ner + 1], lc1[0:Ner + 1])
             cor2 = trapz(lc2[-Ner - 1::], lc1[-Ner - 1::])
@@ -402,27 +385,11 @@ class LevelCrossings(WafoData):
             cor1 = 0
             cor2 = 0
         
-    
         lc22 = hstack((0, cumtrapz(lc2, lc1) + cor1))
         lc22 = (lc22 + 0.5) / (lc22[-1] + cor2 + 1)
         lc11 = (lc1 - mean) / sigma
         
-        # find the mode
-        imin = abs(lc22 - 0.15).argmin()
-        imax = abs(lc22 - 0.85).argmin()
-        
-        inde = slice(imin, imax + 1)
-        lc222 = SmoothSpline(lc11[inde], g22[inde], opt.csm, opt.linextrap, cvar[inde])(lc11[inde])
-        
-        #tmp = smooth(cros(inde,1),g2(inde,2),opt.csm,cros(inde,1),def,cvar(inde));
-        
-        imax = lc222.argmax()
-        u0 = lc22[inde][imax]
-        #u0 = interp1q(cros(:,2),cros(:,1),.5)
-        
-        
-        #lc22 = ndtri(lc22) - u0 #
-        lc22 = invnorm(lc22) - u0
+        lc22 = invnorm(lc22)  #- ymean
         
         g2 = TrData(lc22.copy(), lc1.copy(), mean=mean, sigma=sigma)
         g2.setplotter('step')
@@ -432,9 +399,9 @@ class LevelCrossings(WafoData):
         # to be linear outside the edges or choosing a lower value for csm2.
         
         inds = slice(Ne, ncr - Ne) # indices to points we are smoothing over
-        scros2 = SmoothSpline(lc11[inds], lc22[inds], opt.gsm, opt.linextrap, gvar[inds])(uu)
+        slc22 = SmoothSpline(lc11[inds], lc22[inds], opt.gsm, opt.linextrap, gvar[inds])(uu)
         
-        g = TrData(scros2, g1, mean=mean, sigma=sigma)  #*sa; #multiply with stdev 
+        g = TrData(slc22.copy(), g1.copy(), mean=mean, sigma=sigma)  #*sa; #multiply with stdev 
         
         if opt.chkder:
             for ix in range(5):
@@ -739,7 +706,7 @@ class TimeSeries(WafoData):
     >>> h = rf.plot()
     
     >>> S = ts.tospecdata()
-    The default L is set to 68
+    The default L is set to 325
     
     >>> tp = ts.turning_points()
     >>> mm = tp.cycle_pairs()
@@ -1215,29 +1182,20 @@ class TimeSeries(WafoData):
         >>> g1, g1emp = ts.trdata(method='m', gvar=0.5 ) # Equal weight on all points
         >>> g2, g2emp = ts.trdata(method='n', gvar=[3.5, 0.5, 3.5])  # Less weight on the ends
         >>> int(S.tr.dist2gauss()*100)
-        593
+        141
         >>> int(g0emp.dist2gauss()*100)
-        439
+        217949
         >>> int(g0.dist2gauss()*100)
-        432
+        93
         >>> int(g1.dist2gauss()*100)
-        234
+        66
         >>> int(g2.dist2gauss()*100)
-        437
-        
-         Hm0 = 7;
-         S = jonswap([],Hm0); g=ochitr([],[Hm0/4]); 
-         S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
-         xs = spec2sdat(S,2^13);
-         g0 = dat2tr(xs,[],'plot','iter');             % Monitor the development
-         g1 = dat2tr(xs,'mnon','gvar', .5 );           % More weight on all points
-         g2 = dat2tr(xs,'nonl','gvar', [3.5 .5 3.5]);  % Less weight on the ends
-         hold on, trplot(g1,g)                                   % Check the fit
-         trplot(g2)
-        
+        84
+         
         See also
         --------
-        troptset, lc2tr, cdf2tr, trplot
+        LevelCrossings.trdata 
+        wafo.transform.models
            
         References
         ----------
@@ -1259,7 +1217,7 @@ class TimeSeries(WafoData):
         #    'cvar',1,'gvar',1,'multip',0);
         
         opt = DotDict(chkder=True, plotflag=False, csm=.95, gsm=.05,
-            param=[-5, 5, 513], delay=2, ntr=inf, linextrap=True, ne=7, cvar=1, gvar=1,
+            param=[-5, 5, 513], delay=2, ntr=1000, linextrap=True, ne=7, cvar=1, gvar=1,
             multip=False, crossdef='uM')
         opt.update(**options)
         
diff --git a/pywafo/src/wafo/spectrum/models.py b/pywafo/src/wafo/spectrum/models.py
index 5bfa945..75e0858 100644
--- a/pywafo/src/wafo/spectrum/models.py
+++ b/pywafo/src/wafo/spectrum/models.py
@@ -97,9 +97,10 @@ def _gengamspec(wn, N=5, M=4):
 
     Examples
     --------
+    >>> import wafo.spectrum.models as wsm
     >>> import numpy as np
     >>> wn = np.linspace(0,4,5)
-    >>> _gengamspec(wn, N=6, M=2)
+    >>> wsm._gengamspec(wn, N=6, M=2)
     array([ 0.        ,  1.16765216,  0.17309961,  0.02305179,  0.00474686])
 
     See also
@@ -232,7 +233,8 @@ class Bretschneider(ModelSpectrum):
 
     Examples
     --------
-    >>> S = Bretschneider(Hm0=6.5,Tp=10)
+    >>> import wafo.spectrum.models as wsm
+    >>> S = wsm.Bretschneider(Hm0=6.5,Tp=10)
     >>> S((0,1,2,3))
     array([ 0.        ,  1.69350993,  0.06352698,  0.00844783])
 
@@ -289,21 +291,22 @@ def jonswap_peakfact(Hm0, Tp):
 
     Examples
     --------
+    >>> import wafo.spectrum.models as wsm
     >>> import pylab as plb
     >>> Tp,Hs = plb.meshgrid(range(4,8),range(2,6))
-    >>> gam = jonswap_peakfact(Hs,Tp)
+    >>> gam = wsm.jonswap_peakfact(Hs,Tp)
 
     >>> Hm0 = plb.linspace(1,20)
     >>> Tp = Hm0
     >>> [T,H] = plb.meshgrid(Tp,Hm0)
-    >>> gam = jonswap_peakfact(H,T)
+    >>> gam = wsm.jonswap_peakfact(H,T)
     >>> v = plb.arange(0,8)
     >>> h = plb.contourf(Tp,Hm0,gam,v);h=plb.colorbar()
 
     >>> Hm0 = plb.arange(1,11)
     >>> Tp  = plb.linspace(2,16)
     >>> T,H = plb.meshgrid(Tp,Hm0)
-    >>> gam = jonswap_peakfact(H,T)
+    >>> gam = wsm.jonswap_peakfact(H,T)
     >>> h = plb.plot(Tp,gam.T)
     >>> h = plb.xlabel('Tp [s]')
     >>> h = plb.ylabel('Peakedness parameter')
@@ -356,8 +359,9 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'):
 
     Example
     --------
+    >>> import wafo.spectrum.models as wsm
     >>> fetch = 10000; u10 = 10
-    >>> ss = jonswap_seastate(u10, fetch, output='dict')
+    >>> ss = wsm.jonswap_seastate(u10, fetch, output='dict')
     >>> ss
     {'Ag': 0.016257903375341734,
      'Hm0': 0.51083679198275533,
@@ -365,13 +369,13 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'):
      'gamma': 2.4824142635861119,
      'sigmaA': 0.075317331395172021,
      'sigmaB': 0.091912084512251344}
-    >>> S = Jonswap(**ss)
+    >>> S = wsm.Jonswap(**ss)
     >>> S.Hm0
     0.51083679198275533
 
     # Alternatively
-    >>> ss1 = jonswap_seastate(u10, fetch, output='list')
-    >>> S1 = Jonswap(*ss1)
+    >>> ss1 = wsm.jonswap_seastate(u10, fetch, output='list')
+    >>> S1 = wsm.Jonswap(*ss1)
     >>> S1.Hm0
     0.51083679198275533
 
@@ -484,12 +488,14 @@ class Jonswap(ModelSpectrum):
     Examples
     ---------
     >>> import pylab as plb
-    >>> S = Jonswap(Hm0=7, Tp=11,gamma=1)
+    >>> import wafo.spectrum.models as wsm
+    >>> S = wsm.Jonswap(Hm0=7, Tp=11,gamma=1)
     >>> w = plb.linspace(0,5)
     >>> h = plb.plot(w,S(w))
 
-    >>> S2 = Bretschneider(Hm0=7, Tp=11)
-    >>> assert(all(abs(S(w)-S2(w))<1.e-7),'JONSWAP with gamma=1 differs from Bretscneider, should be equal!')
+    >>> S2 = wsm.Bretschneider(Hm0=7, Tp=11)
+    >>> all(abs(S(w)-S2(w))<1.e-7)
+    True
     >>> plb.close('all')
 
     See also
@@ -652,10 +658,10 @@ def phi1(wi, h, g=9.81):
     Example:
     -------
     Transform a JONSWAP spectrum to a spectrum for waterdepth = 30 m
-
-    >>> S = Jonswap()
+    >>> import wafo.spectrum.models as wsm
+    >>> S = wsm.Jonswap()
     >>> w = range(3.0)
-    >>> S(w)*phi1(w,30.0)
+    >>> S(w)*wsm.phi1(w,30.0)
     array([ 0.        ,  1.0358056 ,  0.03796281])
 
 
@@ -723,9 +729,10 @@ class Tmaspec(Jonswap):
 
     Example
     --------
+    >>> import wafo.spectrum.models as wsm
     >>> import pylab as plb
     >>> w = plb.linspace(0,2.5)
-    >>> S = Tmaspec(h=10,gamma=1) # Bretschneider spectrum Hm0=7, Tp=11
+    >>> S = wsm.Tmaspec(h=10,gamma=1) # Bretschneider spectrum Hm0=7, Tp=11
     >>> o=plb.plot(w,S(w))
     >>> o=plb.plot(w,S(w,h=21))
     >>> o=plb.plot(w,S(w,h=42))
@@ -812,9 +819,10 @@ class Torsethaugen(ModelSpectrum):
 
     Example
     -------
+    >>> import wafo.spectrum.models as wsm
     >>> import pylab as plb
     >>> w = plb.linspace(0,4)
-    >>> S = Torsethaugen(Hm0=6, Tp=8)
+    >>> S = wsm.Torsethaugen(Hm0=6, Tp=8)
     >>> h=plb.plot(w,S(w),w,S.wind(w),w,S.swell(w))
 
     See also
@@ -1034,7 +1042,8 @@ class McCormick(Bretschneider):
 
     Example:
     --------
-    >>> S = McCormick(Hm0=6.5,Tp=10)
+    >>> import wafo.spectrum.models as wsm
+    >>> S = wsm.McCormick(Hm0=6.5,Tp=10)
     >>> S(range(4))
     array([ 0.        ,  1.87865908,  0.15050447,  0.02994663])
 
@@ -1107,7 +1116,8 @@ class OchiHubble(ModelSpectrum):
 
     Examples
     --------
-    >>> S = OchiHubble(par=2)
+    >>> import wafo.spectrum.models as wsm
+    >>> S = wsm.OchiHubble(par=2)
     >>> S(range(4))
     array([ 0.        ,  0.90155636,  0.04185445,  0.00583207])
 
@@ -1239,7 +1249,8 @@ class Wallop(Bretschneider):
 
     Example:
     --------
-    >>> S = Wallop(Hm0=6.5, Tp=10)
+    >>> import wafo.spectrum.models as wsm
+    >>> S = wsm.Wallop(Hm0=6.5, Tp=10)
     >>> S(range(4))
     array([  0.00000000e+00,   9.36921871e-01,   2.76991078e-03,
              7.72996150e-05])
@@ -1360,8 +1371,9 @@ class Spreading(object):
 
     Examples
     --------
+    >>> import wafo.spectrum.models as wsm
     >>> import pylab as plb
-    >>> D = Spreading('cos2s',s_a=10.0)
+    >>> D = wsm.Spreading('cos2s',s_a=10.0)
 
     >>> w = plb.linspace(0,3,257)
     >>> theta = plb.linspace(-pi,pi,129)
@@ -1369,13 +1381,13 @@ class Spreading(object):
 
        # Make frequency dependent direction
     >>> theta0 = lambda w: w*plb.pi/6.0
-    >>> D2 = Spreading('cos2s',theta0=theta0)
+    >>> D2 = wsm.Spreading('cos2s',theta0=theta0)
     >>> t = plb.contour(D2(theta,w)[0])
 
     # Plot all spreading functions
     >>> alltypes = ('cos2s','box','mises','poisson','sech2','wrap_norm')
     >>> for ix in range(len(alltypes)):
-    ...     D3 = Spreading(alltypes[ix])
+    ...     D3 = wsm.Spreading(alltypes[ix])
     ...     t = plb.figure(ix)
     ...     t = plb.contour(D3(theta,w)[0])
     ...     t = plb.title(alltypes[ix])
@@ -1927,10 +1939,11 @@ class Spreading(object):
         
          Example 
          -------
-         >>> S = Jonswap().tospecdata()
-         >>> D = Spreading('cos2s')
+         >>> import wafo.spectrum.models as wsm
+         >>> S = wsm.Jonswap().tospecdata()
+         >>> D = wsm.Spreading('cos2s')
          >>> SD = D.tospecdata2d(S)
-         >>> SD.plot() 
+         >>> h = SD.plot() 
           
          See also  spreading, rotspec, jonswap, torsethaugen  
         '''  
diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py
index 4b4988d..55faf42 100644
--- a/pywafo/src/wafo/stats/estimation.py
+++ b/pywafo/src/wafo/stats/estimation.py
@@ -837,9 +837,10 @@ class FitDistribution(rv_frozen):
         plotbackend.subplot(2, 2, 2)
         self.plotepdf()
         plotbackend.subplot(2, 2, 3)
-        self.plotresprb()
-        plotbackend.subplot(2, 2, 4)
         self.plotresq()
+        plotbackend.subplot(2, 2, 4)
+        self.plotresprb()
+        
         fixstr = ''
         if not self.par_fix == None:
             numfix = len(self.i_fixed)
diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py
index 392856e..d05a98b 100644
--- a/pywafo/src/wafo/transform/core.py
+++ b/pywafo/src/wafo/transform/core.py
@@ -145,15 +145,16 @@ class TrData(WafoData, TrCommon):
     -------
     Construct a linear transformation model
     >>> import numpy as np
+    >>> import wafo.transorm as wt
     >>> sigma = 5; mean = 1
     >>> u = np.linspace(-5,5); x = sigma*u+mean; y = u
-    >>> g = TrData(y,x)
+    >>> g = wt.TrData(y,x)
     >>> g.mean
     array([ 1.])
     >>> g.sigma
     array([ 5.])
     
-    >>> g = TrData(y,x,mean=1,sigma=5)
+    >>> g = wt.TrData(y,x,mean=1,sigma=5)
     >>> g.mean
     1
     >>> g.sigma
diff --git a/pywafo/src/wafo/transform/models.py b/pywafo/src/wafo/transform/models.py
index 815f557..d59afd5 100644
--- a/pywafo/src/wafo/transform/models.py
+++ b/pywafo/src/wafo/transform/models.py
@@ -26,19 +26,20 @@ from core import TrCommon, TrData
 __all__ = ['TrHermite', 'TrLinear', 'TrOchi']
 
 _example = '''
+    >>> import numpy as np
+    >>> import wafo.spectrum.models as sm
+    >>> import wafo.transform.models as tm
     >>> std = 7./4
-    >>> g = <generic>(sigma=std, ysigma=std)
+    >>> g = tm.<generic>(sigma=std, ysigma=std)
     
     Simulate a Transformed Gaussian process:
-    >>> import numpy as np
-    >>> import wafo.spectrum.models as sm
     >>> Sj = sm.Jonswap(Hm0=4*std, Tp=11)
     >>> w = np.linspace(0,4,256)
-    >>> S = Sj.tospecdata(w) #Make spectrum object from numerical values
+    >>> S = Sj.tospecdata(w) # Make spectrum object from numerical values
     >>> ys = S.sim(ns=15000) # Simulated in the Gaussian world
     
     >>> me, va, sk, ku = S.stats_nl(moments='mvsk')
-    >>> g2 = <generic>(mean=me, var=va, skew=sk, kurt=ku, ysigma=std)
+    >>> g2 = tm.<generic>(mean=me, var=va, skew=sk, kurt=ku, ysigma=std)
     >>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world 
     '''
 class TrCommon2(TrCommon):
@@ -111,7 +112,9 @@ class TrHermite(TrCommon2):
     --------
     """ + _example.replace('<generic>', 'TrHermite') + """
     >>> g.dist2gauss()
-     3.9858776379926808
+    0.88230868748851499
+    >>> g2.dist2gauss()
+    1.1411663205144991
     
     See also
     --------  
@@ -337,6 +340,8 @@ class TrLinear(TrCommon2):
     """ + _example.replace('<generic>', 'TrLinear') + """
     >>> g.dist2gauss()
     0.0
+    >>> g2.dist2gauss()
+    3.8594770921678001e-31
     
     See also
     --------  
@@ -399,7 +404,9 @@ class TrOchi(TrCommon2):
     -------
     """ + _example.replace('<generic>', 'TrOchi') + """
     >>> g.dist2gauss()
-    5.9322684525265501
+    1.410698801056657
+    >>> g2.dist2gauss()
+    1.988807188766706
     
     See also
     --------
@@ -486,7 +493,8 @@ class TrOchi(TrCommon2):
         mean = self.mean
         sigma = self.sigma
         xn = atleast_1d(x)
-        xn = (xn - mean) / sigma
+        shape0 = xn.shape
+        xn = (xn.ravel() - mean) / sigma
         igp, = where(0 <= xn)
         igm, = where(xn < 0)
         
@@ -498,6 +506,7 @@ class TrOchi(TrCommon2):
         if gb != 0:
             np.put(g, igm, (-expm1(-gb * xn[igm])) / gb)
 
+        g.shape = shape0
         return (g - mean2) * self.ysigma / sigma2 + self.ymean
 
     def _gauss2dat(self, y, *yi):
@@ -509,7 +518,7 @@ class TrOchi(TrCommon2):
         sigma = self.sigma
         
         yn = (atleast_1d(y) - self.ymean) / self.ysigma
-        xn = sigma2 * yn + mean2
+        xn = sigma2 * yn.ravel() + mean2
 
         igp, = where(0 <= xn)
         igm, = where(xn < 0)
@@ -519,7 +528,8 @@ class TrOchi(TrCommon2):
 
         if gb != 0:
             np.put(xn, igm, -log1p(-gb * xn[igm]) / gb)
-
+        
+        xn.shape = yn.shape
         return  sigma * xn + mean
 
 def main():
diff --git a/pywafo/src/wafo/wafodata.py b/pywafo/src/wafo/wafodata.py
index 443d2e0..b456d41 100644
--- a/pywafo/src/wafo/wafodata.py
+++ b/pywafo/src/wafo/wafodata.py
@@ -2,6 +2,7 @@ import warnings
 from plotbackend import plotbackend
 from time import gmtime, strftime
 import numpy as np
+from scipy.integrate.quadrature import cumtrapz #@UnresolvedImport
 
 __all__ = ['WafoData', 'AxisLabels']
 
@@ -66,19 +67,21 @@ class WafoData(object):
         self.data = data
         self.args = args
         self.date = now()
-        self.plotter = None
+        self.plotter = kwds.pop('plotter', None)
         self.children = None
-        self.plot_args_children = kwds.get('plot_args_children',[])
-        self.plot_kwds_children = kwds.get('plot_kwds_children',{})
-        self.plot_args = kwds.get('plot_args',[])
-        self.plot_kwds = kwds.get('plot_kwds',{})
+        self.plot_args_children = kwds.pop('plot_args_children',[])
+        self.plot_kwds_children = kwds.pop('plot_kwds_children',{})
+        self.plot_args = kwds.pop('plot_args',[])
+        self.plot_kwds = kwds.pop('plot_kwds',{})
         
         self.labels = AxisLabels(**kwds)
-        self.setplotter()
+        if not self.plotter:
+            self.setplotter()
 
     def plot(self, *args, **kwds):
         tmp = None
-        if self.children != None:
+        plotflag = kwds.get('plotflag', None)
+        if not plotflag and self.children != None:
             plotbackend.hold('on')
             tmp = []
             child_args = args if len(args) else tuple(self.plot_args_children)
@@ -106,12 +109,10 @@ class WafoData(object):
         newcopy.__dict__.update(self.__dict__)
         return newcopy
 
-
     def setplotter(self, plotmethod=None):
         '''
             Set plotter based on the data type data_1d, data_2d, data_3d or data_nd
         '''
-
         if isinstance(self.args, (list, tuple)): # Multidimensional data
             ndim = len(self.args)
             if ndim < 2:
@@ -140,10 +141,6 @@ class AxisLabels:
         newcopy = empty_copy(self)
         newcopy.__dict__.update(self.__dict__)
         return newcopy
-        #lbkwds = self.labels.__dict__.copy()
-        #labels = AxisLabels(**lbkwds)
-        #return labels
-
     def labelfig(self):
         try:
             h1 = plotbackend.title(self.title)
@@ -178,7 +175,6 @@ class Plotter_1d(object):
             plotmethod = 'plot'
         self.plotbackend = plotbackend
         try:
-            #self.plotfun = plotbackend.__dict__[plotmethod]
             self.plotfun = getattr(plotbackend, plotmethod)
         except:
             pass
@@ -190,9 +186,136 @@ class Plotter_1d(object):
             args1 = tuple((wdata.args)) + (wdata.data,) + args
         else:
             args1 = tuple((wdata.args,)) + (wdata.data,) + args
-        h1 = self.plotfun(*args1, **kwds)
+        plotflag = kwds.pop('plotflag', None)
+        if plotflag:
+            h1 = self._plot(plotflag, *args1, **kwds)
+        else:
+            h1 = self.plotfun(*args1, **kwds)
         h2 = wdata.labels.labelfig()
         return h1, h2
+    
+    def _plot(self, plotflag, *args1, **kwds):
+        x = args1[0]
+        data = transformdata(x,args1[1], plotflag)
+        dataCI = ()
+        h1 = plot1d(x,data, dataCI, plotflag, *args1[2:], **kwds)
+        return h1
+def plot1d(args,data,dataCI,plotflag,*varargin, **kwds):
+     
+    plottype = np.mod(plotflag,10)
+    if plottype==0: # %  No plotting
+        return []
+    elif plottype==1: 
+        H = plotbackend.plot(args,data,*varargin, **kwds)
+    elif plottype==2: 
+        H = plotbackend.step(args,data,*varargin, **kwds)
+    elif plottype==3: 
+        H = plotbackend.stem(args,data,*varargin, **kwds)
+    elif plottype==4: 
+        H = plotbackend.errorbar(args,data,dataCI[:,0]-data,dataCI[:,1]-data,*varargin, **kwds)
+    elif plottype==5: 
+        H = plotbackend.bar(args,data,*varargin, **kwds)
+    elif plottype==6:  
+        level = 0
+        if np.isfinite(level):
+            H = plotbackend.fill_between(args,data,level,*varargin, **kwds);
+        else:
+            H = plotbackend.fill_between(args,data,*varargin, **kwds);
+        
+    scale = plotscale(plotflag);
+    logXscale = any(scale=='x');
+    logYscale = any(scale=='y');
+    logZscale = any(scale=='z');
+    ax = plotbackend.gca()
+    if logXscale: 
+        plotbackend.setp(ax,xscale='log')
+    if logYscale:
+        plotbackend.setp(ax,yscale='log') 
+    if logZscale: 
+        plotbackend.setp(ax,zscale='log')
+    
+    transFlag = np.mod(plotflag//10,10)
+    logScale = logXscale or logYscale or logZscale
+    if  logScale or (transFlag ==5 and  not logScale):
+        ax = list(plotbackend.axis())
+        fmax1 = data.max()
+        if transFlag==5 and not logScale:
+            ax[3] = 11*np.log10(fmax1)
+            ax[2] = ax[3]-40
+        else:
+            ax[3] = 1.15*fmax1;
+            ax[2] = ax[3]*1e-4;
+      
+        plotbackend.axis(ax)
+    
+    if dataCI and plottype < 3:
+        plotbackend.hold('on')
+       
+        plot1d(args,dataCI,(),plotflag,'r--');
+    return H
+
+def plotscale(plotflag):
+    '''
+    Return plotscale from plotflag
+    
+     CALL scale = plotscale(plotflag)
+    
+     plotflag = integer defining plotscale.
+       Let scaleId = floor(plotflag/100). 
+       If scaleId < 8 then:
+          0 'linear' : Linear scale on all axes.
+          1 'xlog'   : Log scale on x-axis.
+          2 'ylog'   : Log scale on y-axis.
+          3 'xylog'  : Log scale on xy-axis.
+          4 'zlog'   : Log scale on z-axis.
+          5 'xzlog'  : Log scale on xz-axis.
+          6 'yzlog'  : Log scale on yz-axis.
+          7 'xyzlog' : Log scale on xyz-axis.
+      otherwise
+       if (mod(scaleId,10)>0)            : Log scale on x-axis.
+       if (mod(floor(scaleId/10),10)>0)  : Log scale on y-axis.
+       if (mod(floor(scaleId/100),10)>0) : Log scale on z-axis.
+    
+     scale    = string defining plotscale valid options are:
+           'linear', 'xlog', 'ylog', 'xylog', 'zlog', 'xzlog',
+           'yzlog',  'xyzlog' 
+    
+     Example
+     plotscale(100)  % xlog
+     plotscale(200)  % xlog
+     plotscale(1000) % ylog
+    
+     See also plotscale
+    '''
+    scaleId = plotflag//100
+    if scaleId>7:
+        logXscaleId = np.mod(scaleId,10)>0
+        logYscaleId = (np.mod(scaleId//10, 10)>0)*2
+        logZscaleId = (np.mod(scaleId//100, 10)>0)*4
+        scaleId = logYscaleId +logXscaleId+logZscaleId
+    
+    scales = ['linear','xlog','ylog','xylog','zlog','xzlog','yzlog','xyzlog']
+    
+    return scales[scaleId]
+
+def transformdata(x,f,plotflag):
+    transFlag = np.mod(plotflag//10,10)    
+    if transFlag==0:
+        data = f
+    elif transFlag==1:
+        data = 1-f
+    elif transFlag==2:
+        data = cumtrapz(f,x)
+    elif transFlag==3:
+        data = 1-cumtrapz(f,x)
+    if transFlag in (4,5):
+        if transFlag==4:
+            data = -np.log1p(-cumtrapz(f,x))
+        else:
+            if any(f<0):
+                raise ValueError('Invalid plotflag: Data or dataCI is negative, but must be positive')     
+            data = 10*np.log10(f)
+    return data
 
 class Plotter_2d(Plotter_1d):
     """
@@ -201,6 +324,7 @@ class Plotter_2d(Plotter_1d):
     plotmethod : string
         defining type of plot. Options are:
         contour (default)
+        contourf
         mesh
         surf
     """
@@ -209,7 +333,9 @@ class Plotter_2d(Plotter_1d):
         if plotmethod is None:
             plotmethod = 'contour'
         super(Plotter_2d, self).__init__(plotmethod)
-        #self.plotfun = plotbackend.__dict__[plotmethod]
+        
+    def _plot(self,*args,**kwds):
+        pass
 
        
 def main():