diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py
index 4e5b7fa..6cac986 100644
--- a/pywafo/src/wafo/objects.py
+++ b/pywafo/src/wafo/objects.py
@@ -1,5 +1,5 @@
 
-#-------------------------------------------------------------------------------
+
 # Name:        module1
 # Purpose:
 #
@@ -8,7 +8,7 @@
 # Created:     16.09.2008
 # Copyright:   (c) pab 2008
 # Licence:     <your licence>
-#-------------------------------------------------------------------------------
+
 #!/usr/bin/env python
 
 
diff --git a/pywafo/src/wafo/ppimport.py b/pywafo/src/wafo/ppimport.py
new file mode 100644
index 0000000..b11fb55
--- /dev/null
+++ b/pywafo/src/wafo/ppimport.py
@@ -0,0 +1,441 @@
+#!/usr/bin/env python
+"""
+Postpone module import to future.
+
+Python versions: 1.5.2 - 2.3.x
+Author: Pearu Peterson <pearu@cens.ioc.ee>
+Created: March 2003
+$Revision: 922 $
+$Date: 2004-11-27 14:23:27 -0700 (Sat, 27 Nov 2004) $
+"""
+__all__ = ['ppimport','ppimport_attr','ppresolve']
+
+import os
+import sys
+import types
+import traceback
+
+DEBUG=0
+
+_ppimport_is_enabled = 1
+def enable():
+    """ Enable postponed importing."""
+    global _ppimport_is_enabled
+    _ppimport_is_enabled = 1
+
+def disable():
+    """ Disable postponed importing."""
+    global _ppimport_is_enabled
+    _ppimport_is_enabled = 0
+
+class PPImportError(ImportError):
+    pass
+
+def _get_so_ext(_cache={}):
+    so_ext = _cache.get('so_ext')
+    if so_ext is None:
+        if sys.platform[:5]=='linux':
+            so_ext = '.so'
+        else:
+            try:
+                # if possible, avoid expensive get_config_vars call
+                from distutils.sysconfig import get_config_vars
+                so_ext = get_config_vars('SO')[0] or ''
+            except ImportError:
+                #XXX: implement hooks for .sl, .dll to fully support
+                #     Python 1.5.x
+                so_ext = '.so'
+        _cache['so_ext'] = so_ext
+    return so_ext
+
+def _get_frame(level=0):
+    try:
+        return sys._getframe(level+1)
+    except AttributeError:
+        # Python<=2.0 support
+        frame = sys.exc_info()[2].tb_frame
+        for i in range(level+1):
+            frame = frame.f_back
+        return frame
+
+def ppimport_attr(module, name):
+    """ ppimport(module, name) is 'postponed' getattr(module, name)
+    """
+    global _ppimport_is_enabled
+    if _ppimport_is_enabled and isinstance(module, _ModuleLoader):
+        return _AttrLoader(module, name)
+    return getattr(module, name)
+
+class _AttrLoader(object):
+    def __init__(self, module, name):
+        self.__dict__['_ppimport_attr_module'] = module
+        self.__dict__['_ppimport_attr_name'] = name
+
+    def _ppimport_attr_getter(self):
+        module = self.__dict__['_ppimport_attr_module']
+        if isinstance(module, _ModuleLoader):
+            # in case pp module was loaded by other means
+            module = sys.modules[module.__name__]
+        attr = getattr(module,
+                       self.__dict__['_ppimport_attr_name'])
+        try:
+            d = attr.__dict__
+            if d is not None:
+                self.__dict__ = d
+        except AttributeError:
+            pass
+        self.__dict__['_ppimport_attr'] = attr
+
+        return attr
+
+    def __nonzero__(self):
+        return 1
+
+    def __getattr__(self, name):
+        try:
+            attr = self.__dict__['_ppimport_attr']
+        except KeyError:
+            attr = self._ppimport_attr_getter()
+        if name=='_ppimport_attr':
+            return attr
+        return getattr(attr, name)
+
+    def __repr__(self):
+        if '_ppimport_attr' in self.__dict__:
+            return repr(self._ppimport_attr)
+        module = self.__dict__['_ppimport_attr_module']
+        name = self.__dict__['_ppimport_attr_name']
+        return "<attribute %s of %s>" % (`name`,`module`)
+
+    __str__ = __repr__
+
+    # For function and class attributes.
+    def __call__(self, *args, **kwds):
+        return self._ppimport_attr(*args,**kwds)
+
+
+
+def _is_local_module(p_dir,name,suffices):
+    base = os.path.join(p_dir,name)
+    for suffix in suffices:
+        if os.path.isfile(base+suffix):
+            if p_dir:
+                return base+suffix
+            return name+suffix
+
+def ppimport(name):
+    """ ppimport(name) -> module or module wrapper
+
+    If name has been imported before, return module. Otherwise
+    return ModuleLoader instance that transparently postpones
+    module import until the first attempt to access module name
+    attributes.
+    """
+    global _ppimport_is_enabled
+
+    level = 1
+    parent_frame = p_frame = _get_frame(level)
+    while '__name__' not in p_frame.f_locals:
+        level = level + 1
+        p_frame = _get_frame(level)
+
+    p_name = p_frame.f_locals['__name__']
+    if p_name=='__main__':
+        p_dir = ''
+        fullname = name
+    elif '__path__' in p_frame.f_locals:
+        # python package
+        p_path = p_frame.f_locals['__path__']
+        p_dir = p_path[0]
+        fullname = p_name + '.' + name
+    else:
+        # python module
+        p_file = p_frame.f_locals['__file__']
+        p_dir = os.path.dirname(p_file)
+        fullname = p_name + '.' + name
+
+    # module may be imported already
+    module = sys.modules.get(fullname)
+    if module is not None:
+        if _ppimport_is_enabled or isinstance(module, types.ModuleType):
+            return module
+        return module._ppimport_importer()
+
+    so_ext = _get_so_ext()
+    py_exts = ('.py','.pyc','.pyo')
+    so_exts = (so_ext,'module'+so_ext)
+
+    for d,n,fn,e in [\
+        # name is local python module or local extension module
+        (p_dir, name, fullname, py_exts+so_exts),
+        # name is local package
+        (os.path.join(p_dir, name), '__init__', fullname, py_exts),
+        # name is package in parent directory (scipy specific)
+        (os.path.join(os.path.dirname(p_dir), name), '__init__', name, py_exts),
+        ]:
+        location = _is_local_module(d, n, e)
+        if location is not None:
+            fullname = fn
+            break
+
+    if location is None:
+        # name is to be looked in python sys.path.
+        fullname = name
+        location = 'sys.path'
+
+    # Try once more if module is imported.
+    # This covers the case when importing from python module
+    module = sys.modules.get(fullname)
+
+    if module is not None:
+        if _ppimport_is_enabled or isinstance(module,types.ModuleType):
+            return module
+        return module._ppimport_importer()
+    # It is OK if name does not exists. The ImportError is
+    # postponed until trying to use the module.
+
+    loader = _ModuleLoader(fullname,location,p_frame=parent_frame)
+    if _ppimport_is_enabled:
+        return loader
+
+    return loader._ppimport_importer()
+
+def _get_frame_code(frame):
+    filename = frame.f_code.co_filename
+    lineno = frame.f_lineno
+    result = '%s in %s:\n' % (filename,frame.f_code.co_name)
+    if not os.path.isfile(filename):
+        return result
+    f = open(filename)
+    i = 1
+    line = f.readline()
+    while line:
+        line = f.readline()
+        i = i + 1
+        if (abs(i-lineno)<2):
+            result += '#%d: %s\n' % (i,line.rstrip())
+        if i>lineno+3:
+            break
+    f.close()
+    return result
+
+def frame_traceback(frame):
+    if not frame:
+        return
+    blocks = []
+    f = frame
+    while f:
+        blocks.insert(0,_get_frame_code(f))
+        f = f.f_back
+    print '='*50
+    print '\n'.join(blocks)
+    print '='*50
+
+class _ModuleLoader(object):
+    # Don't use it directly. Use ppimport instead.
+
+    def __init__(self,name,location,p_frame=None):
+
+        # set attributes, avoid calling __setattr__
+        self.__dict__['__name__'] = name
+        self.__dict__['__file__'] = location
+        self.__dict__['_ppimport_p_frame'] = p_frame
+
+        if location != 'sys.path':
+            from numpy.testing import Tester
+            self.__dict__['test'] = Tester(os.path.dirname(location)).test
+
+        # install loader
+        sys.modules[name] = self
+
+    def _ppimport_importer(self):
+        name = self.__name__
+
+        try:
+            module = sys.modules[name]
+        except KeyError:
+            raise ImportError,self.__dict__.get('_ppimport_exc_info')[1]
+        if module is not self:
+            exc_info = self.__dict__.get('_ppimport_exc_info')
+            if exc_info is not None:
+                raise PPImportError,\
+                      ''.join(traceback.format_exception(*exc_info))
+            else:
+                assert module is self,`(module, self)`
+
+        # uninstall loader
+        del sys.modules[name]
+
+        if DEBUG:
+            print 'Executing postponed import for %s' %(name)
+        try:
+            module = __import__(name,None,None,['*'])
+        except Exception,msg: # ImportError:
+            if DEBUG:
+                p_frame = self.__dict__.get('_ppimport_p_frame',None)
+                frame_traceback(p_frame)
+            self.__dict__['_ppimport_exc_info'] = sys.exc_info()
+            raise
+
+        assert isinstance(module,types.ModuleType),`module`
+
+        #self.__dict__ = module.__dict__
+        self.__dict__.update(module.__dict__)
+        self.__dict__['_ppimport_module'] = module
+
+        # XXX: Should we check the existence of module.test? Warn?
+        from numpy.testing import Tester
+        test = Tester(os.path.dirname(module.__file__)).test
+        return module
+
+    def __setattr__(self, name, value):
+        try:
+            module = self.__dict__['_ppimport_module']
+        except KeyError:
+            module = self._ppimport_importer()
+        return setattr(module, name, value)
+
+    def __getattr__(self, name):
+        try:
+            module = self.__dict__['_ppimport_module']
+        except KeyError:
+            module = self._ppimport_importer()
+        return getattr(module, name)
+
+    def __repr__(self):
+        global _ppimport_is_enabled
+        if not _ppimport_is_enabled:
+            try:
+                module = self.__dict__['_ppimport_module']
+            except KeyError:
+                module = self._ppimport_importer()
+            return module.__repr__()
+        if '_ppimport_module' in self.__dict__:
+            status = 'imported'
+        elif '_ppimport_exc_info' in self.__dict__:
+            status = 'import error'
+        else:
+            status = 'import postponed'
+        return '<module %s from %s [%s]>' \
+               % (`self.__name__`,`self.__file__`, status)
+
+    __str__ = __repr__
+
+def ppresolve(a,ignore_failure=None):
+    """ Return resolved object a.
+
+    a can be module name, postponed module, postponed modules
+    attribute, string representing module attribute, or any
+    Python object.
+    """
+    global _ppimport_is_enabled
+    if _ppimport_is_enabled:
+        disable()
+        a = ppresolve(a,ignore_failure=ignore_failure)
+        enable()
+        return a
+    if type(a) is type(''):
+        ns = a.split('.')
+        if ignore_failure:
+            try:
+                a = ppimport(ns[0])
+            except:
+                return a
+        else:
+            a = ppimport(ns[0])
+        b = [ns[0]]
+        del ns[0]
+        while ns:
+            if hasattr(a,'_ppimport_importer') or \
+                   hasattr(a,'_ppimport_module'):
+                a = getattr(a,'_ppimport_module',a)
+            if hasattr(a,'_ppimport_attr'):
+                a = a._ppimport_attr
+            b.append(ns[0])
+            del ns[0]
+            if ignore_failure and not hasattr(a, b[-1]):
+                a = '.'.join(ns+b)
+                b = '.'.join(b)
+                if b in sys.modules and sys.modules[b] is None:
+                    del sys.modules[b]
+                return a
+            a = getattr(a,b[-1])
+    if hasattr(a,'_ppimport_importer') or \
+           hasattr(a,'_ppimport_module'):
+        a = getattr(a,'_ppimport_module',a)
+    if hasattr(a,'_ppimport_attr'):
+        a = a._ppimport_attr
+    return a
+
+def _ppresolve_ignore_failure(a):
+    return ppresolve(a,ignore_failure=1)
+
+try:
+    import pydoc as _pydoc
+except ImportError:
+    _pydoc = None
+
+if _pydoc is not None:
+    # Redefine __call__ method of help.__class__ to
+    # support ppimport.
+    import new as _new
+
+    _old_pydoc_help_call = _pydoc.help.__class__.__call__
+    def _ppimport_pydoc_help_call(self,*args,**kwds):
+        return _old_pydoc_help_call(self, *map(_ppresolve_ignore_failure,args),
+                                    **kwds)
+    _ppimport_pydoc_help_call.__doc__ = _old_pydoc_help_call.__doc__
+    _pydoc.help.__class__.__call__ = _new.instancemethod(_ppimport_pydoc_help_call,
+                                                         None,
+                                                         _pydoc.help.__class__)
+
+    _old_pydoc_Doc_document = _pydoc.Doc.document
+    def _ppimport_pydoc_Doc_document(self,*args,**kwds):
+        args = (_ppresolve_ignore_failure(args[0]),) + args[1:]
+        return _old_pydoc_Doc_document(self,*args,**kwds)
+    _ppimport_pydoc_Doc_document.__doc__ = _old_pydoc_Doc_document.__doc__
+    _pydoc.Doc.document = _new.instancemethod(_ppimport_pydoc_Doc_document,
+                                              None,
+                                              _pydoc.Doc)
+
+    _old_pydoc_describe = _pydoc.describe
+    def _ppimport_pydoc_describe(object):
+        return _old_pydoc_describe(_ppresolve_ignore_failure(object))
+    _ppimport_pydoc_describe.__doc__ = _old_pydoc_describe.__doc__
+    _pydoc.describe = _ppimport_pydoc_describe
+
+import inspect as _inspect
+_old_inspect_getfile = _inspect.getfile
+def _ppimport_inspect_getfile(object):
+    if isinstance(object,_ModuleLoader):
+        return object.__dict__['__file__']
+    return _old_inspect_getfile(_ppresolve_ignore_failure(object))
+_ppimport_inspect_getfile.__doc__ = _old_inspect_getfile.__doc__
+_inspect.getfile = _ppimport_inspect_getfile
+
+_old_inspect_getdoc = _inspect.getdoc
+def _ppimport_inspect_getdoc(object):
+    return _old_inspect_getdoc(_ppresolve_ignore_failure(object))
+_ppimport_inspect_getdoc.__doc__ = _old_inspect_getdoc.__doc__
+_inspect.getdoc = _ppimport_inspect_getdoc
+
+_old_inspect_getsource = _inspect.getsource
+def _ppimport_inspect_getsource(object):
+    return _old_inspect_getsource(_ppresolve_ignore_failure(object))
+_ppimport_inspect_getsource.__doc__ = _old_inspect_getsource.__doc__
+_inspect.getsource = _ppimport_inspect_getsource
+
+import __builtin__ as _builtin
+_old_builtin_dir = _builtin.dir
+def _ppimport_builtin_dir(*arg):
+    if not arg:
+        p_frame = _get_frame(1)
+        g = p_frame.f_globals
+        l = p_frame.f_locals
+        l['_ppimport_old_builtin_dir'] = _old_builtin_dir
+        r = eval('_ppimport_old_builtin_dir()',g,l)
+        del r[r.index('_ppimport_old_builtin_dir')]
+        return r
+    return _old_builtin_dir(*map(_ppresolve_ignore_failure,arg))
+_ppimport_builtin_dir.__doc__ = _old_builtin_dir.__doc__
+_builtin.dir = _ppimport_builtin_dir
diff --git a/pywafo/src/wafo/spectrum/models.py b/pywafo/src/wafo/spectrum/models.py
index 139c1ad..1e02867 100644
--- a/pywafo/src/wafo/spectrum/models.py
+++ b/pywafo/src/wafo/spectrum/models.py
@@ -56,14 +56,14 @@ from dispersion_relation import w2k
 #ppimport.enable()
 #_wafospectrum = ppimport.ppimport('wafo.spectrum')
 from core import SpecData1D
-sech = lambda x: 1.0/cosh(x)
+sech = lambda x: 1.0 / cosh(x)
 
 eps = finfo(float).eps
 
 
-__all__ = ['Bretschneider','Jonswap','Torsethaugen','Wallop','McCormick','OchiHubble',
-            'Tmaspec','jonswap_peakfact','jonswap_seastate','spreading',
-            'w2k','k2w','phi1']
+__all__ = ['Bretschneider', 'Jonswap', 'Torsethaugen', 'Wallop', 'McCormick', 'OchiHubble',
+            'Tmaspec', 'jonswap_peakfact', 'jonswap_seastate', 'spreading',
+            'w2k', 'k2w', 'phi1']
 
 # Model spectra
 
@@ -121,17 +121,17 @@ def _gengamspec(wn, N=5, M=4):
     S = zeros_like(w)
 
     ##for w>0 # avoid division by zero
-    k = flatnonzero(w>0.0)
-    if k.size>0:
-        B  = N/M
-        C  = (N-1.0)/M
+    k = flatnonzero(w > 0.0)
+    if k.size > 0:
+        B = N / M
+        C = (N - 1.0) / M
 
         #     # A = Normalizing factor related to Bretschneider form
         #     A    = B**C*M/gamma(C)
         #     S(k) = A*wn(k)**(-N)*exp(-B*wn(k)**(-M))
         logwn = log(w.take(k))
-        logA = (C*log(B)+log(M)-sp.gammaln(C))
-        S.put(k,exp(logA-N*logwn-B*exp(-M*logwn)))
+        logA = (C * log(B) + log(M) - sp.gammaln(C))
+        S.put(k, exp(logA - N * logwn - B * exp(-M * logwn)))
     return S
 
 class ModelSpectrum(object):
@@ -160,9 +160,9 @@ class ModelSpectrum(object):
 
         if w is None:
             if wc is None:
-                wc = 33./self.Tp
-            w = linspace(0,wc,nw)
-        S = SpecData1D(self.__call__(w),w)
+                wc = 33. / self.Tp
+            w = linspace(0, wc, nw)
+        S = SpecData1D(self.__call__(w), w)
         try:
             h = self.h
             S.h = h
@@ -176,13 +176,13 @@ class ModelSpectrum(object):
         ''' Check if seastate is valid
         '''
 
-        if self.Hm0<0:
+        if self.Hm0 < 0:
             raise ValueError('Hm0 can not be negative!')
 
-        if self.Tp<=0:
+        if self.Tp <= 0:
             raise ValueError('Tp must be positve!')
 
-        if self.Hm0==0.0:
+        if self.Hm0 == 0.0:
             warnings.warn('Hm0 is zero!')
 
         self._chk_extra_param()
@@ -253,19 +253,19 @@ class Bretschneider(ModelSpectrum):
             self.chk_seastate()
 
 
-    def __call__(self,wi):
+    def __call__(self, wi):
         ''' Return Bretschnieder spectrum
         '''
         w = atleast_1d(wi)
-        if self.Hm0>0:
-            wp = 2*pi/self.Tp
-            wn = w/wp
-            S = (self.Hm0/4.0)**2/wp * _gengamspec(wn,self.N,self.M)
+        if self.Hm0 > 0:
+            wp = 2 * pi / self.Tp
+            wn = w / wp
+            S = (self.Hm0 / 4.0) ** 2 / wp * _gengamspec(wn, self.N, self.M)
         else:
             S = zeros_like(w)
         return S
 
-def jonswap_peakfact(Hm0,Tp):
+def jonswap_peakfact(Hm0, Tp):
     ''' Jonswap peakedness factor, gamma, given Hm0 and Tp
 
     Parameters
@@ -316,17 +316,17 @@ def jonswap_peakfact(Hm0,Tp):
     --------
     jonswap
     '''
-    Hm0,Tp = atleast_1d(Hm0,Tp)
+    Hm0, Tp = atleast_1d(Hm0, Tp)
 
-    x = Tp/sqrt(Hm0)
+    x = Tp / sqrt(Hm0)
 
     gam = ones_like(x)
 
-    k1 = flatnonzero(x<=5.14285714285714)
-    if k1.size>0: # #limiting gamma to [1 7]
+    k1 = flatnonzero(x <= 5.14285714285714)
+    if k1.size > 0: # #limiting gamma to [1 7]
         xk = x.take(k1)
-        D  = 0.036-0.0056*xk # # approx 5.061*Hm0**2/Tp**4*(1-0.287*log(gam))
-        gam.put(k1,minimum(exp(3.484*( 1.0-0.1975*D*xk**4.0 ) ),7.0)) # # gamma
+        D = 0.036 - 0.0056 * xk # # approx 5.061*Hm0**2/Tp**4*(1-0.287*log(gam))
+        gam.put(k1, minimum(exp(3.484 * (1.0 - 0.1975 * D * xk ** 4.0)), 7.0)) # # gamma
 
     return gam
 def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'):
@@ -401,35 +401,35 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'):
 
     # The following formulas are from Lewis and Allos 1990:
 
-    zeta    = g*fetch/(u10**2) # dimensionless fetch, Table 1
+    zeta = g * fetch / (u10 ** 2) # dimensionless fetch, Table 1
     #zeta = min(zeta, 2.414655013429281e+004)
     if method.startswith('h'):
-        if method[-1]=='3': # Hasselman et.al (1973)
-            A = 0.076*zeta**(-0.22)
-            ny= 3.5*zeta**(-0.33)         # dimensionless peakfrequency, Table 1
-            epsilon1 = 9.91e-8*zeta**1.1 # dimensionless surface variance, Table 1
+        if method[-1] == '3': # Hasselman et.al (1973)
+            A = 0.076 * zeta ** (-0.22)
+            ny = 3.5 * zeta ** (-0.33)         # dimensionless peakfrequency, Table 1
+            epsilon1 = 9.91e-8 * zeta ** 1.1 # dimensionless surface variance, Table 1
         else: # Hasselman et.al (1976)
-            A = 0.0662*zeta**(-0.2)
-            ny = 2.84*zeta**(-0.3)  # dimensionless peakfrequency, Table 1
-            epsilon1 = 1.6e-7*zeta      # dimensionless surface variance, Eq.4
+            A = 0.0662 * zeta ** (-0.2)
+            ny = 2.84 * zeta ** (-0.3)  # dimensionless peakfrequency, Table 1
+            epsilon1 = 1.6e-7 * zeta      # dimensionless surface variance, Eq.4
 
-        sa  = 0.07
-        sb  = 0.09
+        sa = 0.07
+        sb = 0.09
         gam = 3.3
     else:
-        A = 0.074*zeta**(-0.22)     # Eq. 10
-        ny = 3.57*zeta**(-0.33)     # dimensionless peakfrequency, Eq. 11
-        epsilon1 = 3.512e-4*A*ny**(-4.)*zeta**(-0.1)  # dimensionless surface variance, Eq.12
-        sa  = 0.05468*ny**(-0.32)      # Eq. 13
-        sb  = 0.078314*ny**(-0.16)     # Eq. 14
-        gam = maximum(17.54*zeta**(-0.28384),1)     # Eq. 15
-
-    Tp = u10/(ny*g)                          # Table 1
-    Hm0 = 4*sqrt(epsilon1)*u10**2./g            # Table 1
-    if output[0]=='l':
-        return Hm0,Tp,gam,sa,sb,A
+        A = 0.074 * zeta ** (-0.22)     # Eq. 10
+        ny = 3.57 * zeta ** (-0.33)     # dimensionless peakfrequency, Eq. 11
+        epsilon1 = 3.512e-4 * A * ny ** (-4.) * zeta ** (-0.1)  # dimensionless surface variance, Eq.12
+        sa = 0.05468 * ny ** (-0.32)      # Eq. 13
+        sb = 0.078314 * ny ** (-0.16)     # Eq. 14
+        gam = maximum(17.54 * zeta ** (-0.28384), 1)     # Eq. 15
+
+    Tp = u10 / (ny * g)                          # Table 1
+    Hm0 = 4 * sqrt(epsilon1) * u10 ** 2. / g            # Table 1
+    if output[0] == 'l':
+        return Hm0, Tp, gam, sa, sb, A
     else:
-        return dict(Hm0=Hm0,Tp=Tp,gamma=gam,sigmaA=sa,sigmaB=sb,Ag=A)
+        return dict(Hm0=Hm0, Tp=Tp, gamma=gam, sigmaA=sa, sigmaB=sb, Ag=A)
 
 class Jonswap(ModelSpectrum):
     ''' 
@@ -527,8 +527,8 @@ class Jonswap(ModelSpectrum):
         self.method = method
         self.wnc = wnc
 
-        if self.gamma==None or not isfinite(self.gamma) or self.gamma<1:
-            self.gamma = jonswap_peakfact(Hm0,Tp)
+        if self.gamma == None or not isfinite(self.gamma) or self.gamma < 1:
+            self.gamma = jonswap_peakfact(Hm0, Tp)
 
         self._preCalculateAg()
 
@@ -539,7 +539,7 @@ class Jonswap(ModelSpectrum):
         Tp = self.Tp
         Hm0 = self.Hm0
         gam = self.gamma
-        outsideJonswapRange  = Tp>5*sqrt(Hm0) or Tp<3.6*sqrt(Hm0)
+        outsideJonswapRange = Tp > 5 * sqrt(Hm0) or Tp < 3.6 * sqrt(Hm0)
         if outsideJonswapRange:
             txt0 = '''
             Hm0,Tp is outside the JONSWAP range.
@@ -547,7 +547,7 @@ class Jonswap(ModelSpectrum):
             '''
             warnings.warn(txt0)
 
-        if gam<1 or 7 < gam:
+        if gam < 1 or 7 < gam:
             txt = '''
             The peakedness factor, gamma, is possibly too large. 
             The validity of the spectral density is questionable.
@@ -555,30 +555,30 @@ class Jonswap(ModelSpectrum):
             warnings.warn(txt)
 
 
-    def _localspec(self,wn):
+    def _localspec(self, wn):
         Gf = self.peak_e_factor(wn)
-        return Gf*_gengamspec(wn,self.N, self.M)
+        return Gf * _gengamspec(wn, self.N, self.M)
 
     def _preCalculateAg(self):
         ''' PRECALCULATEAG Precalculate normalization.
         '''
-        if self.gamma==1:
+        if self.gamma == 1:
             self.Ag = 1.0
             self.method = 'parametric'
         elif self.Ag != None:
             self.method = 'custom'
-            if self.Ag<=0:
+            if self.Ag <= 0:
                 raise ValueError('Ag must be larger than 0!')
-        elif self.method[0]=='i':
+        elif self.method[0] == 'i':
             # normalizing by integration
             self.method = 'integration'
-            if self.wnc<1.0:
+            if self.wnc < 1.0:
                 raise ValueError('Normalized cutoff frequency, wnc, must be larger than one!')
             area1, unused_err1 = integrate.quad(self._localspec, 0, 1)
             area2, unused_err2 = integrate.quad(self._localspec, 1, self.wnc)
             area = area1 + area2
-            self.Ag = 1.0/area
-        elif self.method[1]=='p':
+            self.Ag = 1.0 / area
+        elif self.method[1] == 'p':
             self.method = 'parametric'
             ##   # Original normalization
             ##   # NOTE: that  Hm0**2/16 generally is not equal to intS(w)dw
@@ -587,11 +587,11 @@ class Jonswap(ModelSpectrum):
             N = self.N
             M = self.M
             gammai = self.gamma
-            parametersOK = (3<=N and N<=50) or  (2<=M and M <=9.5) and (1<= gammai and gammai<=20)
+            parametersOK = (3 <= N and N <= 50) or  (2 <= M and M <= 9.5) and (1 <= gammai and gammai <= 20)
             if parametersOK:
-                f1NM = 4.1*(N-2*M**0.28+5.3)**(-1.45*M**0.1+0.96)
-                f2NM = (2.2*M**(-3.3) + 0.57)*N**(-0.58*M**0.37+0.53)-1.04*M**(-1.9)+0.94
-                self.Ag = (1+ f1NM*log(gammai)**f2NM)/gammai
+                f1NM = 4.1 * (N - 2 * M ** 0.28 + 5.3) ** (-1.45 * M ** 0.1 + 0.96)
+                f2NM = (2.2 * M ** (-3.3) + 0.57) * N ** (-0.58 * M ** 0.37 + 0.53) - 1.04 * M ** (-1.9) + 0.94
+                self.Ag = (1 + f1NM * log(gammai) ** f2NM) / gammai
 
             ###    elseif N == 5 && M == 4,
             ###      options.Ag = (1+1.0*log(gammai).**1.16)/gammai
@@ -602,35 +602,35 @@ class Jonswap(ModelSpectrum):
             else:
                 raise ValueError('Not knowing the normalization because N, M or peakedness parameter is out of bounds!')
 
-            if self.sigmaA!=0.07 or self.sigmaB!=0.09:
+            if self.sigmaA != 0.07 or self.sigmaB != 0.09:
                 warnings.warn('Use integration to calculate Ag when sigmaA~=0.07 or sigmaB~=0.09')
 
 
 
-    def peak_e_factor(self,wn):
+    def peak_e_factor(self, wn):
         ''' PEAKENHANCEMENTFACTOR
         '''
-        w = maximum(atleast_1d(wn),0.0)
-        sab = where(w>1,self.sigmaB,self.sigmaA)
+        w = maximum(atleast_1d(wn), 0.0)
+        sab = where(w > 1, self.sigmaB, self.sigmaA)
 
-        wnm12 = 0.5*((w-1.0)/sab)**2.0
-        Gf    = self.gamma**(exp(-wnm12))
+        wnm12 = 0.5 * ((w - 1.0) / sab) ** 2.0
+        Gf = self.gamma ** (exp(-wnm12))
         return Gf
 
-    def  __call__(self,wi):
+    def  __call__(self, wi):
         ''' JONSWAP spectral density
         '''
         w = atleast_1d(wi)
-        if (self.Hm0>0.0):
+        if (self.Hm0 > 0.0):
 
-            N   = self.N
-            M   = self.M
-            wp  = 2*pi/self.Tp
-            wn  = w/wp
-            Ag  = self.Ag
+            N = self.N
+            M = self.M
+            wp = 2 * pi / self.Tp
+            wn = w / wp
+            Ag = self.Ag
             Hm0 = self.Hm0
-            Gf  = self.peak_e_factor(wn)
-            S   = ((Hm0/4.0)**2/wp*Ag)*Gf*_gengamspec(wn,N,M)
+            Gf = self.peak_e_factor(wn)
+            S = ((Hm0 / 4.0) ** 2 / wp * Ag) * Gf * _gengamspec(wn, N, M)
         else:
             S = zeros_like(w)
         return S
@@ -674,14 +674,14 @@ def phi1(wi, h, g=9.81):
     if h == inf: # % special case infinite water depth
         return ones_like(w)
 
-    k1 = w2k(w,0,inf,g=g)[0]
-    dw1 = 2.0*w/g  # % dw/dk|h=inf
-    k2 = w2k(w,0,h,g=g)[0]
+    k1 = w2k(w, 0, inf, g=g)[0]
+    dw1 = 2.0 * w / g  # % dw/dk|h=inf
+    k2 = w2k(w, 0, h, g=g)[0]
 
-    k2h = k2*h
-    dw2 = where(k1!=0,dw1/(tanh(k2h)+k2h/cosh(k2h)**2.0),0) # dw/dk|h=h0
+    k2h = k2 * h
+    dw2 = where(k1 != 0, dw1 / (tanh(k2h) + k2h / cosh(k2h) ** 2.0), 0) # dw/dk|h=h0
 
-    return where(k1!=0,(k1/k2)**3.0*dw2/dw1,0)
+    return where(k1 != 0, (k1 / k2) ** 3.0 * dw2 / dw1, 0)
 
 class Tmaspec(Jonswap):
     ''' JONSWAP spectrum for finite water depth
@@ -755,24 +755,24 @@ class Tmaspec(Jonswap):
     Zeitschrift.
 
     '''
-    def __init__(self, Hm0=7.0, Tp=11.0, gamma=None, sigmaA=0.07, sigmaB=0.09, 
-                 Ag=None, N=5, M=4, method='integration', wnc=6.0, 
+    def __init__(self, Hm0=7.0, Tp=11.0, gamma=None, sigmaA=0.07, sigmaB=0.09,
+                 Ag=None, N=5, M=4, method='integration', wnc=6.0,
                  chk_seastate=True, h=42, g=9.81):
         self.g = g
         self.h = h
-        super(Tmaspec, self).__init__(Hm0,Tp,gamma,sigmaA,sigmaB,Ag,N,M,method,wnc,chk_seastate)
+        super(Tmaspec, self).__init__(Hm0, Tp, gamma, sigmaA, sigmaB, Ag, N, M, method, wnc, chk_seastate)
         self.type = 'TMA'
 
-    def phi(self,w,h=None,g=None):
-        if h==None:
+    def phi(self, w, h=None, g=None):
+        if h == None:
             h = self.h
-        if g==None:
+        if g == None:
             g = self.g
-        return phi1(w,h,g)
+        return phi1(w, h, g)
 
-    def __call__(self,w,h=None,g=None):
+    def __call__(self, w, h=None, g=None):
         jonswap = super(Tmaspec, self).__call__(w)
-        return jonswap*self.phi(w,h,g)
+        return jonswap * self.phi(w, h, g)
 
 class Torsethaugen(ModelSpectrum):
     ''' 
@@ -863,20 +863,20 @@ class Torsethaugen(ModelSpectrum):
 
         self._init_spec()
 
-    def __call__(self,w):
+    def __call__(self, w):
         ''' TORSETHAUGEN spectral density
         '''
-        return self.wind(w)+self.swell(w)
+        return self.wind(w) + self.swell(w)
 
     def _chk_extra_param(self):
         Hm0 = self.Hm0
         Tp = self.Tp
-        if Hm0>11 or Hm0>max((Tp/3.6)**2, (Tp-2)*12/11):
+        if Hm0 > 11 or Hm0 > max((Tp / 3.6) ** 2, (Tp - 2) * 12 / 11):
             txt0 = '''Hm0 is outside the valid range.
                     The validity of the spectral density is questionable'''
             warnings.warn(txt0)
 
-        if Tp>20 or Tp<3:
+        if Tp > 20 or Tp < 3:
             txt1 = '''Tp is outside the valid range.
                     The validity of the spectral density is questionable'''
             warnings.warn(txt1)
@@ -916,33 +916,33 @@ class Torsethaugen(ModelSpectrum):
         # the parameters on the line below can be dependent on geographical location
         A10 = 0.7; A1 = 0.5; A20 = 0.6; A2 = 0.3; A3 = 6
 
-        Tf = Af*(Hm0)**(1.0/3.0)
-        Tl = AL*sqrt(Hm0)   # lower limit
+        Tf = Af * (Hm0) ** (1.0 / 3.0)
+        Tl = AL * sqrt(Hm0)   # lower limit
         Tu = Au             # upper limit
 
         #Non-dimensional scales
         # New call pab April 2005
-        El = min(max((Tf-Tp)/(Tf-Tl),0),1) #wind sea
-        Eu = min(max((Tp-Tf)/(Tu-Tf),0),1) #Swell
+        El = min(max((Tf - Tp) / (Tf - Tl), 0), 1) #wind sea
+        Eu = min(max((Tp - Tf) / (Tu - Tf), 0), 1) #Swell
 
 
 
-        if Tp<Tf: # Wind dominated seas
+        if Tp < Tf: # Wind dominated seas
             # Primary peak (wind dominated)
-            Nw  = K0*sqrt(Hm0)+K00             # high frequency exponent
-            Mw  = M0                           # spectral width exponent
-            Rpw = min((1-A10)*exp(-(El/A1)**2)+A10,1)
-            Hpw = Rpw*Hm0                      # significant waveheight wind
+            Nw = K0 * sqrt(Hm0) + K00             # high frequency exponent
+            Mw = M0                           # spectral width exponent
+            Rpw = min((1 - A10) * exp(-(El / A1) ** 2) + A10, 1)
+            Hpw = Rpw * Hm0                      # significant waveheight wind
             Tpw = Tp                           # primary peak period
             # peak enhancement factor
-            gammaw = KG*(1+KG0*exp(-Hm0/KG1))*(2*pi/gravity1*Rpw*Hm0/(Tp**2))**r
-            gammaw = max(gammaw,1)
+            gammaw = KG * (1 + KG0 * exp(-Hm0 / KG1)) * (2 * pi / gravity1 * Rpw * Hm0 / (Tp ** 2)) ** r
+            gammaw = max(gammaw, 1)
             # Secondary peak (swell)
-            Ns  = Nw                # high frequency exponent
-            Ms  = Mw                # spectral width exponent
-            Rps = sqrt(1.0-Rpw**2.0)
-            Hps = Rps*Hm0           # significant waveheight swell
-            Tps = Tf+B1
+            Ns = Nw                # high frequency exponent
+            Ms = Mw                # spectral width exponent
+            Rps = sqrt(1.0 - Rpw ** 2.0)
+            Hps = Rps * Hm0           # significant waveheight swell
+            Tps = Tf + B1
             gammas = 1.0
 
             if monitor:
@@ -953,29 +953,29 @@ class Torsethaugen(ModelSpectrum):
         else: #swell dominated seas
 
             # Primary peak (swell)
-            Ns  = K0*sqrt(Hm0)+K00             #high frequency exponent
-            Ms  = M0                           #spectral width exponent
-            Rps = min((1-A20)*exp(-(Eu/A2)**2)+A20,1)
-            Hps = Rps*Hm0                      # significant waveheight swell
+            Ns = K0 * sqrt(Hm0) + K00             #high frequency exponent
+            Ms = M0                           #spectral width exponent
+            Rps = min((1 - A20) * exp(-(Eu / A2) ** 2) + A20, 1)
+            Hps = Rps * Hm0                      # significant waveheight swell
             Tps = Tp                           # primary peak period
             # peak enhancement factor
-            gammas = KG*(1+KG0*exp(-Hm0/KG1))*(2*pi/gravity1*Hm0/(Tf**2))**r*(1+A3*Eu)
-            gammas = max(gammas,1)
+            gammas = KG * (1 + KG0 * exp(-Hm0 / KG1)) * (2 * pi / gravity1 * Hm0 / (Tf ** 2)) ** r * (1 + A3 * Eu)
+            gammas = max(gammas, 1)
 
             # Secondary peak (wind)
-            Nw   = Ns                       # high frequency exponent
-            Mw   = M0*(1-B2*exp(-Hm0/B3))   # spectral width exponent
-            Rpw  = sqrt(1-Rps**2)
-            Hpw  = Rpw*Hm0                  # significant waveheight wind
-
-            C = (Nw-1)/Mw
-            B = Nw/Mw
-            G0w    = B**C*Mw/sp.gamma(C)#normalizing factor
+            Nw = Ns                       # high frequency exponent
+            Mw = M0 * (1 - B2 * exp(-Hm0 / B3))   # spectral width exponent
+            Rpw = sqrt(1 - Rps ** 2)
+            Hpw = Rpw * Hm0                  # significant waveheight wind
+
+            C = (Nw - 1) / Mw
+            B = Nw / Mw
+            G0w = B ** C * Mw / sp.gamma(C)#normalizing factor
             #G0w = exp(C*log(B)+log(Mw)-gammaln(C))
             #G0w  = Mw/((B)**(-C)*gamma(C))
 
-            if Hpw>0:
-                Tpw  = (16*S0*(1-exp(-Hm0/S1))*(0.4)**Nw/( G0w*Hpw**2))**(-1.0/(Nw-1.0))
+            if Hpw > 0:
+                Tpw = (16 * S0 * (1 - exp(-Hm0 / S1)) * (0.4) ** Nw / (G0w * Hpw ** 2)) ** (-1.0 / (Nw - 1.0))
             else:
                 Tpw = inf
 
@@ -989,13 +989,13 @@ class Torsethaugen(ModelSpectrum):
 
 
         if monitor:
-            if (3.6*sqrt(Hm0)<= Tp & Tp<=5*sqrt(Hm0)):
+            if (3.6 * sqrt(Hm0) <= Tp & Tp <= 5 * sqrt(Hm0)):
                 print('     Jonswap range')
 
             print('Hm0 = %g' % Hm0)
             print('Ns, Ms = %g, %g  Nw, Mw = %g, %g' % (Ns, Ms, Nw, Mw))
-            print('gammas = %g gammaw = ' % (gammas,gammaw))
-            print('Rps = %g Rpw = %g' % (Rps,Rpw))
+            print('gammas = %g gammaw = ' % (gammas, gammaw))
+            print('Rps = %g Rpw = %g' % (Rps, Rpw))
             print('Hps = %g Hpw = %g' % (Hps, Hpw))
             print('Tps = %g Tpw = %g' % (Tps, Tpw))
 
@@ -1003,10 +1003,10 @@ class Torsethaugen(ModelSpectrum):
         #G0s=Ms/((Ns/Ms)**(-(Ns-1)/Ms)*gamma((Ns-1)/Ms )) #normalizing factor
 
         # Wind part
-        self.wind  = Jonswap(Hm0=Hpw, Tp=Tpw, gamma=gammaw, N=Nw, M=Mw, 
+        self.wind = Jonswap(Hm0=Hpw, Tp=Tpw, gamma=gammaw, N=Nw, M=Mw,
                              method=self.method, chk_seastate=False)
         # Swell part
-        self.swell = Jonswap(Hm0=Hps, Tp=Tps, gamma=gammas, N=Ns, M=Ms, 
+        self.swell = Jonswap(Hm0=Hps, Tp=Tps, gamma=gammas, N=Ns, M=Ms,
                              method=self.method, chk_seastate=False)
 
 class McCormick(Bretschneider):
@@ -1060,23 +1060,23 @@ class McCormick(Bretschneider):
         self.type = 'McCormick'
         self.Hm0 = Hm0
         self.Tp = Tp
-        if Tz==None:
-            Tz = 0.8143*Tp
+        if Tz == None:
+            Tz = 0.8143 * Tp
 
         self.Tz = Tz
         if chk_seastate:
             self.chk_seastate()
 
-        if M==None and self.Hm0>0:
-            self._TpdTz = Tp/Tz
-            M = 1.0/optimize.fminbound(self._localoptfun,0.01,5)
+        if M == None and self.Hm0 > 0:
+            self._TpdTz = Tp / Tz
+            M = 1.0 / optimize.fminbound(self._localoptfun, 0.01, 5)
         self.M = M
-        self.N = M+1.0
+        self.N = M + 1.0
 
-    def _localoptfun(self,x):
+    def _localoptfun(self, x):
         #LOCALOPTFUN Local function to optimize.
-        y = 1.0+x
-        return (y**(x)/sp.gamma(y)-self._TpdTz)**2.0
+        y = 1.0 + x
+        return (y ** (x) / sp.gamma(y) - self._TpdTz) ** 2.0
 
 class OchiHubble(ModelSpectrum):
     ''' OchiHubble bimodal spectral density model.
@@ -1141,8 +1141,8 @@ class OchiHubble(ModelSpectrum):
             self.chk_seastate()
         self._init_spec()
 
-    def __call__(self,w):
-        return self.wind(w)+self.swell(w)
+    def __call__(self, w):
+        return self.wind(w) + self.swell(w)
 
 
     def _init_spec(self):
@@ -1156,7 +1156,7 @@ class OchiHubble(ModelSpectrum):
                         [0.65, 0.76],
                         [0.90, 0.44],
                         [0.77, 0.64],
-                        [0.73,0.68],
+                        [0.73, 0.68],
                         [0.92, 0.39]])
         wa = array([ [0.7, 1.15],
                         [0.93, 1.5],
@@ -1180,31 +1180,31 @@ class OchiHubble(ModelSpectrum):
                         [0.039, 0.000],
                         [0.046, 0.039],
                         [0.046, 0.039]])
-        Lpar = array([[3.00, 1.54,-0.062],
-                        [3.00, 2.77,-0.112],
-                        [2.55, 1.82,-0.089],
-                        [2.65, 3.90,-0.085],
-                        [2.60, 0.53,-0.069],
-                        [1.35, 2.48,-0.102],
-                        [4.95, 2.48,-0.102],
-                        [1.80, 2.95,-0.105],
-                        [4.50, 1.95,-0.082],
-                        [6.40, 1.78,-0.069],
-                        [0.70, 1.78,-0.069]])
+        Lpar = array([[3.00, 1.54, -0.062],
+                        [3.00, 2.77, -0.112],
+                        [2.55, 1.82, -0.089],
+                        [2.65, 3.90, -0.085],
+                        [2.60, 0.53, -0.069],
+                        [1.35, 2.48, -0.102],
+                        [4.95, 2.48, -0.102],
+                        [1.80, 2.95, -0.105],
+                        [4.50, 1.95, -0.082],
+                        [6.40, 1.78, -0.069],
+                        [0.70, 1.78, -0.069]])
         Hm0 = self.Hm0
         Lpari = Lpar[self.par]
-        Li = hstack((Lpari[0],Lpari[1]*exp(Lpari[2]*Hm0)))
+        Li = hstack((Lpari[0], Lpari[1] * exp(Lpari[2] * Hm0)))
 
-        Hm0i = hp[self.par]*Hm0
-        Tpi  = 2*pi*exp(wb[self.par]*Hm0)/wa[self.par]
-        Ni   = 4*Li+1
-        Mi   = [4, 4]
+        Hm0i = hp[self.par] * Hm0
+        Tpi = 2 * pi * exp(wb[self.par] * Hm0) / wa[self.par]
+        Ni = 4 * Li + 1
+        Mi = [4, 4]
 
         self.swell = Bretschneider(Hm0=Hm0i[0], Tp=Tpi[0], N=Ni[0], M=Mi[0])
-        self.wind  = Bretschneider(Hm0=Hm0i[1], Tp=Tpi[1], N=Ni[1], M=Mi[1])
+        self.wind = Bretschneider(Hm0=Hm0i[1], Tp=Tpi[1], N=Ni[1], M=Mi[1])
 
     def _chk_extra_param(self):
-        if self.par<0 or 10<self.par:
+        if self.par < 0 or 10 < self.par:
             raise ValueError('Par must be an integer from 0 to 10!')
 
 class Wallop(Bretschneider):
@@ -1266,10 +1266,10 @@ class Wallop(Bretschneider):
         self.Tp = Tp
         self.M = 4
         if N is None:            
-            wp = 2.*pi/Tp
-            kp = w2k(wp,0,inf)[0] # wavenumber at peak frequency
-            Lp = 2.*pi/kp # wave length at the peak frequency
-            N = abs((log(2.*pi**2.)+2*log(Hm0/4)-2.0*log(Lp))/log(2))
+            wp = 2. * pi / Tp
+            kp = w2k(wp, 0, inf)[0] # wavenumber at peak frequency
+            Lp = 2. * pi / kp # wave length at the peak frequency
+            N = abs((log(2. * pi ** 2.) + 2 * log(Hm0 / 4) - 2.0 * log(Lp)) / log(2))
 
         self.N = N
 
@@ -1459,7 +1459,7 @@ class Spreading(object):
 ##% By es, jr 1999.11.25
 
 
-    def __init__(self,type='cos2s',theta0=0,method='mitsuyasu',s_a=15.,s_b=15.,m_a=5.,m_b=-2.5,wn_lo=0.0,wn_c=1.,wn_up=inf):
+    def __init__(self, type='cos2s', theta0=0, method='mitsuyasu', s_a=15., s_b=15., m_a=5., m_b= -2.5, wn_lo=0.0, wn_c=1., wn_up=inf):
 
         self.type = type
         self.theta0 = theta0
@@ -1475,40 +1475,40 @@ class Spreading(object):
         methods = dict(n=None, m='mitsuyasu', d='donelan', b='banner')
         methodslist = (None, 'mitsuyasu', 'donelan', 'banner')
 
-        if isinstance(self.method,str):
+        if isinstance(self.method, str):
             if not self.method[0] in methods:
                 raise ValueError('Unknown method')
             self.method = methods[self.method[0]]
         elif self.method == None:
             pass
         else:
-            if method<0 or 3<method:
+            if method < 0 or 3 < method:
                 method = 2
             self.method = methodslist[method]
 
         self._spreadfun = dict(c=self.cos2s, b=self.box, m=self.mises,
                         p=self.poisson, s=self.sech2, w=self.wrap_norm)
         self._fourierdispatch = dict(b=self.fourier2a, m=self.fourier2k,
-                                    p=self.fourier2x, s=self.fourier2b, 
+                                    p=self.fourier2x, s=self.fourier2b,
                                     w=self.fourier2d)
 
-    def __call__(self,*args):
+    def __call__(self, *args):
         spreadfun = self._spreadfun[self.type[0]]
         return spreadfun(*args)
 
-    def  chk_input(self,theta,w=1,wc=1): # [s_par,TH,phi0,Nt] =
+    def  chk_input(self, theta, w=1, wc=1): # [s_par,TH,phi0,Nt] =
         ''' CHK_INPUT
 
         CALL [s_par,TH,phi0,Nt] = inputchk(theta,w,wc)
         '''
 
-        wn = atleast_1d(w/wc)
+        wn = atleast_1d(w / wc)
         theta = theta.ravel()
         Nt = len(theta)
 
         # Make sure theta is from -pi to pi
         phi0 = 0.0
-        theta = mod(theta+pi,2*pi)-pi
+        theta = mod(theta + pi, 2 * pi) - pi
 
 
         if hasattr(self.theta0, '__call__'):
@@ -1517,32 +1517,32 @@ class Spreading(object):
             th0 = atleast_1d(self.theta0).flatten()
 
         Nt0 = th0.size
-        Nw  = w.size
-        isFreqDepDir = (Nt0==Nw)
+        Nw = w.size
+        isFreqDepDir = (Nt0 == Nw)
         if isFreqDepDir:
             # frequency dependent spreading and/or
             # frequency dependent direction
             # make sure -pi<=TH<pi
-            TH = mod(theta[:,newaxis]-th0[newaxis,:]+pi,2*pi)-pi
-        elif Nt0!=1:
+            TH = mod(theta[:, newaxis] - th0[newaxis, :] + pi, 2 * pi) - pi
+        elif Nt0 != 1:
             raise ValueError('The length of theta0 must equal to 1 or the length of w')
         else:
-            TH = mod(theta-th0+pi,2*pi)-pi  # make sure -pi<=TH<pi
-            if self.method!=None: # frequency dependent spreading
-                TH = TH[:,newaxis]
+            TH = mod(theta - th0 + pi, 2 * pi) - pi  # make sure -pi<=TH<pi
+            if self.method != None: # frequency dependent spreading
+                TH = TH[:, newaxis]
 
         # Get spreading parameter
         #~~~~~~~~~~~~~~~~
         s = self.spread_par_s(wn)
 
-        if self.type[0]=='c': # cos2s
+        if self.type[0] == 'c': # cos2s
             s_par = s
         else:
-            r1 = abs(s/(s+1)) #% First Fourier coefficient of the directional spreading function.
+            r1 = abs(s / (s + 1)) #% First Fourier coefficient of the directional spreading function.
             #% Find distribution parameter from first Fourier coefficient.
             s_par = self.fourier2distpar(r1)
-        if self.method!=None:
-            s_par = s_par[newaxis,:]
+        if self.method != None:
+            s_par = s_par[newaxis, :]
         return s_par, TH, phi0, Nt
 
 
@@ -1569,12 +1569,12 @@ class Spreading(object):
 ##            S = S[newaxis,:]
         gammaln = sp.gammaln
 
-        D = (exp(gammaln(S+1)-gammaln(S+1.0/2.0))/(2*sqrt(pi)))*cos(TH/2.0)**(2.0*S)
+        D = (exp(gammaln(S + 1) - gammaln(S + 1.0 / 2.0)) / (2 * sqrt(pi))) * cos(TH / 2.0) ** (2.0 * S)
         return D, phi0
 
 
 
-    def  poisson(self,theta,w=1,wc=1): # [D,phi0] =
+    def  poisson(self, theta, w=1, wc=1): # [D,phi0] =
         ''' POISSON spreading function
 
         CALL  [D, phi0] = poisson(theta,w,wc)
@@ -1589,17 +1589,17 @@ class Spreading(object):
         phi0 = Parameter defining the actual principal direction of D.
 
         '''
-        [X, TH, phi0, unused_Nt] = self.chk_input(theta,w,wc)
+        [X, TH, phi0, unused_Nt] = self.chk_input(theta, w, wc)
 
 ##        if not self.method[0]=='n':
 ##            X = X[newaxis,:]
 
 
-        D    = (1-X**2.)/(1.-(2.*cos(TH)-X)*X)/(2.*pi)
+        D = (1 - X ** 2.) / (1. - (2. * cos(TH) - X) * X) / (2. * pi)
         return D, phi0
 
 
-    def wrap_norm(self,theta,w=1,wc=1):
+    def wrap_norm(self, theta, w=1, wc=1):
         ''' Wrapped Normal spreading function
 
         CALL  [D, phi0] = wnormal(theta,w,wc)
@@ -1616,27 +1616,27 @@ class Spreading(object):
 
 
 
-        [par,TH,phi0,Nt] = self.chk_input(theta,w,wc)
+        [par, TH, phi0, Nt] = self.chk_input(theta, w, wc)
 
-        D1 = par**2./2.
+        D1 = par ** 2. / 2.
 ##        if not self.method(1)=='n':
 ##            D1 = D1[newaxis,:]
 
 
-        ix  =  arange(1,Nt)
-        ix2 = ix**2
+        ix = arange(1, Nt)
+        ix2 = ix ** 2
         Nd2 = D1.size
-        Fcof = vstack( (ones((1,Nd2))/2, exp(-ix2[:,newaxis]*D1)))/pi
+        Fcof = vstack((ones((1, Nd2)) / 2, exp(-ix2[:, newaxis] * D1))) / pi
 
-        cor = exp(1j*ix[:,newaxis]*TH[0,:])
-        Pcor = vstack( (ones((1,TH.shape[1]) ), cor ) ) # % correction term to get
+        cor = exp(1j * ix[:, newaxis] * TH[0, :])
+        Pcor = vstack((ones((1, TH.shape[1])), cor)) # % correction term to get
         #% the correct integration limits
-        Fcof = Fcof*Pcor.conj()
-        D = real(fft(Fcof,axis=0))
-        D[D<0]=0
+        Fcof = Fcof * Pcor.conj()
+        D = real(fft(Fcof, axis=0))
+        D[D < 0] = 0
         return D, phi0
 
-    def sech2(self,theta,w=1,wc=1):
+    def sech2(self, theta, w=1, wc=1):
         '''Sech**2 spreading function
 
         CALL  [D, phi0] = sech2(theta,w,wc)
@@ -1652,17 +1652,17 @@ class Spreading(object):
         '''
 
 
-        [B, TH, phi0, unused_Nt] = self.chk_input(theta,w,wc)
-        NB = tanh(pi*B) #% Normalization factor.
-        NB = where(NB==0,1.0,NB) # Avoid division by zero
+        [B, TH, phi0, unused_Nt] = self.chk_input(theta, w, wc)
+        NB = tanh(pi * B) #% Normalization factor.
+        NB = where(NB == 0, 1.0, NB) # Avoid division by zero
 ##        if not self.method[0]=='n':
 ##            B = B[newaxis,:]
 ##            NB = NB[newaxis,:]
 
-        D = 0.5*B*sech(B*TH)**2./NB
+        D = 0.5 * B * sech(B * TH) ** 2. / NB
         return D, phi0
 
-    def mises(self,theta,w=1,wc=1):
+    def mises(self, theta, w=1, wc=1):
         '''Mises spreading function
 
         CALL  [D, phi0] = mises(theta,w,wc)
@@ -1677,16 +1677,16 @@ class Spreading(object):
         phi0 = Parameter defining the actual principal direction of D.
         '''
 
-        [K, TH, phi0, unused_Nt] = self.chk_input(theta,w,wc)
+        [K, TH, phi0, unused_Nt] = self.chk_input(theta, w, wc)
 ##        if not self.method[0]=='n':
 ##            K = K[newaxis,:]
 
-        D = exp(K*(cos(TH)-1.))/(2*pi*sp.ive(0,K))
+        D = exp(K * (cos(TH) - 1.)) / (2 * pi * sp.ive(0, K))
         return D, phi0
 
 
 
-    def box(self,theta,w=1,wc=1):
+    def box(self, theta, w=1, wc=1):
         ''' Box car spreading function
 
         CALL  [D, phi0] = box(theta,w,wc)
@@ -1703,17 +1703,17 @@ class Spreading(object):
 
 
 
-        [A,TH,phi0,unused_Nt] = self.chk_input(theta,w,wc)
+        [A, TH, phi0, unused_Nt] = self.chk_input(theta, w, wc)
 
 ##        if not self.method[0]=='n':
 ##            A = A[newaxis,:]
 
-        D = ((-A<=TH) & (TH<=A))/(2.*A)
+        D = ((-A <= TH) & (TH <= A)) / (2. * A)
         return D, phi0
 
 # Local sub functions
 
-    def fourier2distpar(self,r1):
+    def fourier2distpar(self, r1):
         ''' Fourier coefficients to distribution parameter
 
         Parameters
@@ -1743,21 +1743,21 @@ class Spreading(object):
         fourierfun = self._fourierdispatch.get(self.type[0])
         return fourierfun(r1)
 
-    def fourier2x(self,r1):
+    def fourier2x(self, r1):
         ''' Returns the solution of r1 = x.
         '''
         X = r1
-        if any(X>=1):
+        if any(X >= 1):
             raise ValueError('POISSON spreading: X value must be less than 1')
         return X
 
-    def fourier2a(self,r1):
+    def fourier2a(self, r1):
         ''' Returns the solution of R1 = sin(A)/A.
         '''
-        A0 = flipud(linspace(0,pi+0.1,1025))
-        funA = interp1d(sinc(A0/pi),A0)
+        A0 = flipud(linspace(0, pi + 0.1, 1025))
+        funA = interp1d(sinc(A0 / pi), A0)
         A0 = funA(r1.ravel())
-        A  = asarray(A0)
+        A = asarray(A0)
 
         # Newton-Raphson
         da = ones_like(r1)
@@ -1766,18 +1766,18 @@ class Spreading(object):
         ix = flatnonzero(A)
         for unused_iy in range(max_count):
             Ai = A[ix]
-            da[ix] = (sin(Ai) -Ai*r1[ix])/(cos(Ai)-r1[ix])
+            da[ix] = (sin(Ai) - Ai * r1[ix]) / (cos(Ai) - r1[ix])
             Ai = Ai - da[ix]
             # Make sure that the current guess is larger than zero and less than pi
             #x(ix) = xi + 0.1*(dx(ix) - 9*xi).*(xi<=0) + 0.38*(dx(ix)-6.2*xi +6.2).*(xi>pi) 
             # Make sure that the current guess is larger than zero.
-            A[ix] = Ai + 0.5*(da[ix] - Ai)*(Ai<=0.0)
+            A[ix] = Ai + 0.5 * (da[ix] - Ai) * (Ai <= 0.0)
 
-            ix = flatnonzero((abs(da) > sqrt(eps)*abs(A))*(abs(da) > sqrt(eps)))
-            if ix.size==0:
-                if any(A>pi):
+            ix = flatnonzero((abs(da) > sqrt(eps) * abs(A)) * (abs(da) > sqrt(eps)))
+            if ix.size == 0:
+                if any(A > pi):
                     raise ValueError('BOX-CAR spreading: The A value must be less than pi')
-                return A.clip(min=1e-16,max=pi)
+                return A.clip(min=1e-16, max=pi)
 
 
         warnings.warn('Newton raphson method did not converge.')
@@ -1787,52 +1787,52 @@ class Spreading(object):
         ''' 
         Returns the solution of R1 = besseli(1,K)/besseli(0,K),
         '''
-        K0 = hstack((linspace(0,10,513), linspace(10.00001,100)))
-        fun0 = lambda x : sp.ive(1,x)/sp.ive(0,x)
-        funK = interp1d(fun0(K0),K0)
+        K0 = hstack((linspace(0, 10, 513), linspace(10.00001, 100)))
+        fun0 = lambda x : sp.ive(1, x) / sp.ive(0, x)
+        funK = interp1d(fun0(K0), K0)
         K0 = funK(r1.ravel())
         k1 = flatnonzero(isnan(K0))
-        if (k1.size>0):
+        if (k1.size > 0):
             K0[k1] = 0.0
             K0[k1] = K0.max()
 
-        ix0 = flatnonzero(r1!=0.0)
+        ix0 = flatnonzero(r1 != 0.0)
         K = zeros_like(r1)
-        fun = lambda x : fun0(x)-r1[ix]
+        fun = lambda x : fun0(x) - r1[ix]
         for ix in ix0:
-            K[ix] = optimize.fsolve(fun,K0[ix])
+            K[ix] = optimize.fsolve(fun, K0[ix])
         return K
 
-    def fourier2b(self,r1):
+    def fourier2b(self, r1):
         ''' Returns the solution of R1 = pi/(2*B*sinh(pi/(2*B)).
         '''
-        B0 = hstack((linspace(eps,5,513), linspace(5.0001,100)))
-        funB = interp1d(self._r1ofsech2(B0),B0)
+        B0 = hstack((linspace(eps, 5, 513), linspace(5.0001, 100)))
+        funB = interp1d(self._r1ofsech2(B0), B0)
 
         B0 = funB(r1.ravel())
         k1 = flatnonzero(isnan(B0))
-        if (k1.size>0):
+        if (k1.size > 0):
             B0[k1] = 0.0
             B0[k1] = max(B0)
 
-        ix0 = flatnonzero(r1!=0.0)
+        ix0 = flatnonzero(r1 != 0.0)
         B = zeros_like(r1)
-        fun = lambda x : 0.5*pi/(sinh(.5*pi/x))-x*r1[ix]
+        fun = lambda x : 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix]
         for ix in ix0:
-            B[ix] = abs(optimize.fsolve(fun,B0[ix]))
+            B[ix] = abs(optimize.fsolve(fun, B0[ix]))
         return B
 
-    def fourier2d(self,r1):
+    def fourier2d(self, r1):
         ''' Returns the solution of R1 = exp(-D**2/2).
         '''
-        r = clip(r1,0.,1.0)
-        return where(r<=0,inf,sqrt(-2.0*log(r)))
+        r = clip(r1, 0., 1.0)
+        return where(r <= 0, inf, sqrt(-2.0 * log(r)))
 
 
 
 
 
-    def  spread_par_s(self,wn):
+    def  spread_par_s(self, wn):
         ''' Return spread parameter, S, of COS2S function
 
         Parameters
@@ -1844,50 +1844,50 @@ class Spreading(object):
         S   : ndarray
             spread parameter of COS2S functions
         '''
-        if self.method==None:
+        if self.method == None:
             # no frequency dependent spreading,
             # but possible frequency dependent direction
-            s =  atleast_1d(self.s_a)
+            s = atleast_1d(self.s_a)
         else:
             wn_lo = self.wn_lo
             wn_up = self.wn_up
             wn_c = self.wn_c
 
-            spa  = self.s_a
-            spb  = self.s_b
-            ma   = self.m_a
-            mb   = self.m_b
+            spa = self.s_a
+            spb = self.s_b
+            ma = self.m_a
+            mb = self.m_b
 
 
             # Mitsuyasu et. al and Hasselman et. al parametrization   of
             # frequency dependent spreading
-            s = where(wn<=wn_c,spa*wn**ma,spb*wn**mb)
-            s[wn<=wn_lo] = 0.0
+            s = where(wn <= wn_c, spa * wn ** ma, spb * wn ** mb)
+            s[wn <= wn_lo] = 0.0
 
-            k = flatnonzero(wn_up<wn)
-            if k.size>0:
-                if self.method[0]=='d':
+            k = flatnonzero(wn_up < wn)
+            if k.size > 0:
+                if self.method[0] == 'd':
                     # Donelan et. al. parametrization for B in SECH-2
-                    s[k] = spb*(wn_up)**mb
+                    s[k] = spb * (wn_up) ** mb
 
                     # Convert to S-paramater in COS-2S distribution
                     r1 = self.r1ofsech2(s)
-                    s  = r1/(1.-r1)
+                    s = r1 / (1. - r1)
 
-                elif self.method[0]=='b':
+                elif self.method[0] == 'b':
                     # Banner parametrization  for B in SECH-2
-                    s3m   = spb*(wn_up)**mb
-                    s3p   = self._donelan(wn_up)
-                    scale = s3m/s3p #% Scale so that parametrization will be continous
-                    s[k]  = scale*self.donelan(wn[k])
+                    s3m = spb * (wn_up) ** mb
+                    s3p = self._donelan(wn_up)
+                    scale = s3m / s3p #% Scale so that parametrization will be continous
+                    s[k] = scale * self.donelan(wn[k])
                     r1 = self.r1ofsech2(s)
 
                     #% Convert to S-paramater in COS-2S distribution
-                    s  = r1/(1.-r1)
+                    s = r1 / (1. - r1)
                 else:
                     s[k] = 0.0
 
-        if any(s<0):
+        if any(s < 0):
             raise ValueError('The COS2S spread parameter, S(w), value must be larger than 0')
         return s
 
@@ -1895,16 +1895,16 @@ class Spreading(object):
     def _donelan(self, wn):
         ''' High frequency decay of B of sech2 paramater
         '''
-        return  10.0**(-0.4+0.8393*exp(-0.567*log(wn**2)))
+        return  10.0 ** (-0.4 + 0.8393 * exp(-0.567 * log(wn ** 2)))
 
     def _r1ofsech2(self, B):
         ''' R1OFSECH2   Computes R1 = pi./(2*B.*sinh(pi./(2*B)))
         '''
         realmax = finfo(float).max
-        tiny = 1./realmax
-        x = clip(2.*B,tiny,realmax)
-        xk = pi/x
-        return where(x<100.,xk/sinh(xk),-2.*xk/(exp(xk)*expm1(-2.*xk)))
+        tiny = 1. / realmax
+        x = clip(2. * B, tiny, realmax)
+        xk = pi / x
+        return where(x < 100., xk / sinh(xk), -2. * xk / (exp(xk) * expm1(-2. * xk)))
 
 
 
@@ -1912,17 +1912,17 @@ def test_some_spectra():
     S = Jonswap()
 
     w = arange(3.0)
-    S(w)*phi1(w,30.0)
+    S(w) * phi1(w, 30.0)
     S1 = S.tospecdata(w)
     S1.plot()
 
 
     import pylab as plb
-    w = plb.linspace(0,2.5)
-    S = Tmaspec(h=10,gamma=1) # Bretschneider spectrum Hm0=7, Tp=11
-    plb.plot(w,S(w))
-    plb.plot(w,S(w,h=21))
-    plb.plot(w,S(w,h=42))
+    w = plb.linspace(0, 2.5)
+    S = Tmaspec(h=10, gamma=1) # Bretschneider spectrum Hm0=7, Tp=11
+    plb.plot(w, S(w))
+    plb.plot(w, S(w, h=21))
+    plb.plot(w, S(w, h=42))
     plb.show()
     plb.close('all')
 
@@ -1930,39 +1930,39 @@ def test_some_spectra():
 
     import pylab as plb
     #w = plb.linspace(0,3)
-    w,th = plb.ogrid[0:4,0:6]
-    k,k2 = w2k(w,th)
-    k1,k12 = w2k(w,th,h=20)
-    plb.plot(w,k,w,k2)
+    w, th = plb.ogrid[0:4, 0:6]
+    k, k2 = w2k(w, th)
+    k1, k12 = w2k(w, th, h=20)
+    plb.plot(w, k, w, k2)
 
     plb.show()
 
     plb.close('all')
-    w = plb.linspace(0,2,100)
+    w = plb.linspace(0, 2, 100)
     S = Torsethaugen(Hm0=6, Tp=8)
-    plb.plot(w,S(w),w,S.wind(w),w,S.swell(w))
+    plb.plot(w, S(w), w, S.wind(w), w, S.swell(w))
 
-    S1 = Jonswap(Hm0=7, Tp=11,gamma=1)
-    w = plb.linspace(0,2,100)
-    plb.plot(w,S1(w))
+    S1 = Jonswap(Hm0=7, Tp=11, gamma=1)
+    w = plb.linspace(0, 2, 100)
+    plb.plot(w, S1(w))
     plb.show()
     plb.close('all')
 
 
-    Hm0 = plb.arange(1,11)
-    Tp  = plb.linspace(2,16)
-    T,H = plb.meshgrid(Tp,Hm0)
-    gam = jonswap_peakfact(H,T)
-    plb.plot(Tp,gam.T)
+    Hm0 = plb.arange(1, 11)
+    Tp = plb.linspace(2, 16)
+    T, H = plb.meshgrid(Tp, Hm0)
+    gam = jonswap_peakfact(H, T)
+    plb.plot(Tp, gam.T)
     plb.xlabel('Tp [s]')
     plb.ylabel('Peakedness parameter')
 
-    Hm0 = plb.linspace(1,20)
+    Hm0 = plb.linspace(1, 20)
     Tp = Hm0
-    [T,H] = plb.meshgrid(Tp,Hm0)
-    gam = jonswap_peakfact(H,T)
-    v = plb.arange(0,8)
-    plb.contourf(Tp,Hm0,gam,v)
+    [T, H] = plb.meshgrid(Tp, Hm0)
+    gam = jonswap_peakfact(H, T)
+    v = plb.arange(0, 8)
+    plb.contourf(Tp, Hm0, gam, v)
     plb.colorbar()
     plb.show()
     plb.close('all')
@@ -1970,19 +1970,19 @@ def test_some_spectra():
 def test_spreading():
     import pylab as plb
     pi = plb.pi
-    w = plb.linspace(0,3,257)
-    theta = plb.linspace(-pi,pi,129)
-    theta0 = lambda w: w*plb.pi/6.0
-    D2 = Spreading('cos2s',theta0=theta0)
-    d1 =D2(theta,w)[0]
+    w = plb.linspace(0, 3, 257)
+    theta = plb.linspace(-pi, pi, 129)
+    theta0 = lambda w: w * plb.pi / 6.0
+    D2 = Spreading('cos2s', theta0=theta0)
+    d1 = D2(theta, w)[0]
     t = plb.contour(d1.squeeze())
 
     pi = plb.pi
-    D = Spreading('wrap_norm',s_a=10.0)
+    D = Spreading('wrap_norm', s_a=10.0)
 
-    w = plb.linspace(0,3,257)
-    theta = plb.linspace(-pi,pi,129)
-    d1 = D(theta,w)
+    w = plb.linspace(0, 3, 257)
+    theta = plb.linspace(-pi, pi, 129)
+    d1 = D(theta, w)
     plb.contour(d1[0])
     plb.show()
 
@@ -1994,4 +1994,4 @@ def main():
         doctest.testmod()
 
 if __name__ == '__main__':
-    main()
\ No newline at end of file
+    main()
diff --git a/pywafo/src/wafo/stats/__init__.py b/pywafo/src/wafo/stats/__init__.py
index 7588461..f08000f 100644
--- a/pywafo/src/wafo/stats/__init__.py
+++ b/pywafo/src/wafo/stats/__init__.py
@@ -6,7 +6,7 @@ Statistics package in WAFO Toolbox.
 """
 from scipy.stats import *
 from wafo.stats.core import *
-#from wafo.stats.distributions import *
+from wafo.stats.distributions import *
 #from wafo.spectrum.core import SpecData1D
 #import wafo.spectrum.models
 #import wafo.spectrum.dispersion_relation
diff --git a/pywafo/src/wafo/stats/distributions.py b/pywafo/src/wafo/stats/distributions.py
index a3b4b01..42f50a5 100644
--- a/pywafo/src/wafo/stats/distributions.py
+++ b/pywafo/src/wafo/stats/distributions.py
@@ -163,7 +163,7 @@ class rv_frozen(object):
     def __init__(self, dist, *args, **kwds):
         self.dist = dist
         loc0, scale0 = map(kwds.get, ['loc', 'scale'])
-        if isinstance(dist, rv_continuous):
+        if hasattr(dist, 'fix_loc_scale'): #isinstance(dist, rv_continuous):
             args, loc0, scale0 = dist.fix_loc_scale(args, loc0, scale0)
             self.par = args + (loc0, scale0)
         else: # rv_discrete
@@ -569,7 +569,7 @@ class rv_generic(object):
     and rv_continuous.
 
     """
-    def fix_loc_scale(self, args, loc, scale=1):
+    def _fix_loc_scale(self, args, loc, scale=1):
         N = len(args)
         if N > self.numargs:
             if N == self.numargs + 1 and loc is None:
@@ -585,8 +585,8 @@ class rv_generic(object):
             loc = 0.0
         return args, loc, scale
 
-    def fix_loc(self, args, loc):
-        args, loc, scale = self.fix_loc_scale(args, loc)
+    def _fix_loc(self, args, loc):
+        args, loc, scale = self._fix_loc_scale(args, loc)
         return args, loc
 
     # These are actually called, and should not be overwritten if you
@@ -617,7 +617,7 @@ class rv_generic(object):
         loc, scale, size, discrete = map(kwds.get, kwd_names,
                                          [None] * len(kwd_names))
 
-        args, loc, scale = self.fix_loc_scale(args, loc, scale)
+        args, loc, scale = self._fix_loc_scale(args, loc, scale)
         cond = logical_and(self._argcheck(*args), (scale >= 0))
         if not all(cond):
             raise ValueError, "Domain error in arguments."
@@ -765,7 +765,8 @@ class rv_continuous(rv_generic):
                  shapes=None, extradoc=None):
 
         rv_generic.__init__(self)
-
+        self.fix_loc_scale = self._fix_loc_scale
+        
         if badvalue is None:
             badvalue = nan
         self.badvalue = badvalue
@@ -4518,7 +4519,7 @@ class rv_discrete(rv_generic):
                  shapes=None, extradoc=None):
 
         rv_generic.__init__(self)
-
+        self.fix_loc = self._fix_loc
         if badvalue is None:
             badvalue = nan
         self.badvalue = badvalue
@@ -4543,9 +4544,9 @@ class rv_discrete(rv_generic):
             self.pk = take(ravel(self.pk), indx, 0)
             self.a = self.xk[0]
             self.b = self.xk[-1]
-            self.P = make_dict(self.xk, self.pk)
+            self.P = dict(zip(self.xk, self.pk)) #make_dict(self.xk, self.pk)
             self.qvals = numpy.cumsum(self.pk, axis=0)
-            self.F = make_dict(self.xk, self.qvals)
+            self.F = dict(zip(self.xk, self.qvals)) #make_dict(self.xk, self.qvals)
             self.Finv = reverse_dict(self.F)
             self._ppf = new.instancemethod(sgf(_drv_ppf, otypes='d'),
                                            self, rv_discrete)
diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py
index d7c6198..b75db37 100644
--- a/pywafo/src/wafo/stats/estimation.py
+++ b/pywafo/src/wafo/stats/estimation.py
@@ -1,26 +1,23 @@
-# Functions to implement several important functions for
-#   various Continous and Discrete Probability Distributions
-#
-# Author:  Travis Oliphant  2002-2003
+''' 
+Contains FitDistribution and Profile class, which are
 
+important classes for fitting to various Continous and Discrete Probability Distributions
+
+Author:  Per A. Brodtkorb 2008
+'''
 
 from __future__ import division
 
 from wafo.plotbackend import plotbackend
-from wafo.misc import ecross, findcross, JITImport
+from wafo.misc import ecross, findcross
 from scipy.misc.ppimport import ppimport
 
 import numdifftools 
-
+from scipy import special
 from scipy.linalg import pinv2
 
 
 from scipy import optimize
-# Trick to avoid error due to circular import
-#_WAFOCOV = ppimport('wafo.covariance')
-_WAFODIST = ppimport('wafo.stats.distributions')
-#_WAFODIST = JITImport('wafo.stats.distributions')
-
 from numpy import alltrue, arange, \
      ravel, ones, sum, \
      zeros, log, sqrt, exp
@@ -43,6 +40,14 @@ arr = asarray
 
 all = alltrue
 
+def chi2isf(p,df):
+    return special.chdtri(df, p)
+def chi2sf(x, df):
+    return special.chdtrc(df, x)
+
+def norm_ppf(q):
+    return special.ndtri(q)
+
 def valarray(shape, value=nan, typecode=None):
     """Return an array of all value.
     """
@@ -87,7 +92,7 @@ class rv_frozen(object):
     def __init__(self, dist, *args, **kwds):
         self.dist = dist
         loc0, scale0 = map(kwds.get, ['loc', 'scale'])
-        if isinstance(dist, _WAFODIST.rv_continuous):
+        if hasattr(dist, 'fix_loc_scale'): #isinstance(dist, _WAFODIST.rv_continuous):
             args, loc0, scale0 = dist.fix_loc_scale(args, loc0, scale0)
             self.par = args + (loc0, scale0)
         else: # rv_discrete
@@ -247,7 +252,7 @@ class Profile(object):
         self.i_free = nonzero(isfree)
 
         self.Lmax = Lmax
-        self.alpha_Lrange = 0.5 * _WAFODIST.chi2.isf(self.alpha, 1)
+        self.alpha_Lrange = 0.5 * chi2isf(self.alpha, 1) #_WAFODIST.chi2.isf(self.alpha, 1)
         self.alpha_cross_level = Lmax - self.alpha_Lrange
         lowLevel = self.alpha_cross_level - self.alpha_Lrange / 7.0
 
@@ -344,7 +349,8 @@ class Profile(object):
                 pcov = self.fit_dist.par_cov[i_notfixed, :][:, i_notfixed]
                 pvar = sum(numpy.dot(drl, pcov) * drl)
 
-            p_crit = _WAFODIST.norm.isf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5
+            #p_crit = _WAFODIST.norm.isf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5
+            p_crit = -norm_ppf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5
             if self.pmin == None:
                 self.pmin = p_opt - 5.0 * p_crit
             if self.pmax == None:
@@ -392,7 +398,7 @@ class Profile(object):
         if alpha < self.alpha:
             raise ValueError('Unable to return CI with alpha less than %g' % self.alpha)
 
-        cross_level = self.Lmax - 0.5 * _WAFODIST.chi2.isf(alpha, 1)
+        cross_level = self.Lmax - 0.5 * chi2isf(alpha, 1) #_WAFODIST.chi2.isf(alpha, 1)
         ind = findcross(self.data, cross_level)
         N = len(ind)
         if N == 0:
@@ -536,7 +542,7 @@ class FitDistribution(rv_frozen):
             self.par_cov[:, :] = nan
 
         pvar = numpy.diag(self.par_cov)
-        zcrit = -_WAFODIST.norm.ppf(self.alpha / 2.0)
+        zcrit = -norm_ppf(self.alpha / 2.0)#_WAFODIST.norm.ppf(self.alpha / 2.0)
         self.par_lower = self.par - zcrit * sqrt(pvar)
         self.par_upper = self.par + zcrit * sqrt(pvar)
 
@@ -702,7 +708,7 @@ class FitDistribution(rv_frozen):
         Other distribution types will introduce deviations in the plot.
         '''
 
-        bin, limits = numpy.histogram(self.data, normed=True, new=True)
+        bin, limits = numpy.histogram(self.data, normed=True) #, new=True)
         limits.shape = (-1, 1)
         xx = limits.repeat(3, axis=1)
         xx.shape = (-1,)
@@ -793,7 +799,7 @@ class FitDistribution(rv_frozen):
         C1 = m - sqrt(0.5 * n * v)
         C2 = sqrt(v / (2.0 * n))
         Tn = (T + 0.5 * k * isParUnKnown - C1) / C2 # chi2 with n degrees of freedom
-        pvalue = _WAFODIST.chi2.sf(Tn, n)
+        pvalue = chi2sf(Tn, n) #_WAFODIST.chi2.sf(Tn, n)
         return pvalue
 
     def nlogps(self, theta, x):
@@ -964,6 +970,7 @@ class FitDistribution(rv_frozen):
         return - H
 
 def main():
+    _WAFODIST = ppimport('wafo.stats.distributions')
      #nbinom(10, 0.75).rvs(3)
     import matplotlib
     matplotlib.interactive(True)
@@ -986,15 +993,15 @@ def main():
     Lx.get_CI(alpha=0.2)
 
     # CI for logSF=log(SF)
-    Lpr = phat.profile(i=1, logSF=log(R), link=phat.dist.link)
+    Lpr = phat.profile(i=0, logSF=log(R), link=phat.dist.link)
     Lpr.plot()
     Lpr.get_CI(alpha=0.075)
 
-    dlaplace.stats(0.8, loc=0)
+    _WAFODIST.dlaplace.stats(0.8, loc=0)
 #    pass
-    t = planck(0.51000000000000001)
+    t = _WAFODIST.planck(0.51000000000000001)
     t.ppf(0.5)
-    t = zipf(2)
+    t = _WAFODIST.zipf(2)
     t.ppf(0.5)
     import pylab as plb
     _WAFODIST.rice.rvs(1)