diff --git a/pywafo/src/wafo/source/c_codes/c_functions.c b/pywafo/src/wafo/source/c_codes/c_functions.c
index febbc6b..7725902 100644
--- a/pywafo/src/wafo/source/c_codes/c_functions.c
+++ b/pywafo/src/wafo/source/c_codes/c_functions.c
@@ -229,64 +229,64 @@ void disufq(double *rvec, double *ivec,
       kw1  = kw[ix];
       Epij = kw1; 
       for (i=0;i<m;i++,ixi++,iz1++) { 
-	rrA = rA[ixi]*rA[ixi]; ///
-	iiA = iA[ixi]*iA[ixi]; ///
-	riA = rA[ixi]*iA[ixi]; ///
+    	  rrA = rA[ixi]*rA[ixi]; ///
+    	  iiA = iA[ixi]*iA[ixi]; ///
+    	  riA = rA[ixi]*iA[ixi]; ///
 	
-	/// Sum frequency effects along the diagonal
-	tmp1 = kfact*(rrA-iiA)*Epij;
-	tmp2 = kfact*2.0*riA*Epij;
-	rvec[iz1] += tmp1;
-	ivec[iz1] += tmp2;
+    	  /// Sum frequency effects along the diagonal
+    	  tmp1 = kfact*(rrA-iiA)*Epij;
+    	  tmp2 = kfact*2.0*riA*Epij;
+    	  rvec[iz1] += tmp1;
+    	  ivec[iz1] += tmp2;
 	
-	//rvec[iz2] += tmp1;
-	//ivec[iz2] -= tmp2;
-	//iz2++;
+    	  //rvec[iz2] += tmp1;
+    	  //ivec[iz2] -= tmp2;
+    	  //iz2++;
 
-	/// Difference frequency effects are zero along the diagonal
-	/// and are thus not contributing to the mean.
+    	  /// Difference frequency effects are zero along the diagonal
+    	  /// and are thus not contributing to the mean.
       }
       for (jy = ix+1;jy<nmax;jy++){
-	kw2  = kw[jy];
-	Epij = 0.5*(kw2 + kw1); 
-	Edij = -0.5*(kw2 - kw1);
-	//printf("Edij = %f Epij = %f \n", Edij,Epij);
+    	  kw2  = kw[jy];
+    	  Epij = 0.5*(kw2 + kw1);
+    	  Edij = -0.5*(kw2 - kw1);
+    	  //printf("Edij = %f Epij = %f \n", Edij,Epij);
 	
-	ixi = ix*m;
-	jyi = jy*m;
-	iz1 = ixi+jyi;
-	iv1 = jyi-ixi;
-	//iz2 = (n*m-iz1);
-	//iv2 = (n*m-iv1);
-	for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { 
+    	  ixi = ix*m;
+    	  jyi = jy*m;
+    	  iz1 = ixi+jyi;
+    	  iv1 = jyi-ixi;
+    	  //iz2 = (n*m-iz1);
+    	  //iv2 = (n*m-iv1);
+    	  for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
 	  
-	  rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
-	  iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
-	  riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
-	  irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
+    		  rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
+    		  iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
+    		  riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
+    		  irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
 	  
-	 /* Sum frequency effects */ 
-	  tmp1 = kfact*2.0*(rrA-iiA)*Epij;
-	  tmp2 = kfact*2.0*(riA+irA)*Epij;
-	  rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1;
-	  ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2;
-	  //rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] +=  tmp1;
-	  //ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -=  tmp2;
-	  // iz2++;
+    		  /* Sum frequency effects */
+    		  tmp1 = kfact*2.0*(rrA-iiA)*Epij;
+    		  tmp2 = kfact*2.0*(riA+irA)*Epij;
+    		  rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1;
+    		  ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2;
+    		  //rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] +=  tmp1;
+    		  //ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -=  tmp2;
+    		  // iz2++;
 
-	  /* Difference frequency effects */
-	  tmp1 = kfact*2.0*(rrA+iiA)*Edij;
-	  tmp2 = kfact*2.0*(riA-irA)*Edij;
+    		  /* Difference frequency effects */
+    		  tmp1 = kfact*2.0*(rrA+iiA)*Edij;
+    		  tmp2 = kfact*2.0*(riA-irA)*Edij;
 	  
-	  rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
-	  ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
+    		  rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
+    		  ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
 	  
-	  //rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
-	  //ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
-	  //iv2++;
-	}
+    		  //rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
+    		  //ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
+    		  //iv2++;
+    	  }
       }
-   }
+    }
   }
   else{ /* Finite water depth */
     for (ix = nmin-1;ix<nmax;ix++) {
diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py
index 5cf338e..5fb9176 100644
--- a/pywafo/src/wafo/spectrum/core.py
+++ b/pywafo/src/wafo/spectrum/core.py
@@ -1496,15 +1496,11 @@ class SpecData1D(WafoData):
         ...     res = fun(x2[:,1::],axis=0)
         ...     m = res.mean()
         ...     sa = res.std()
-        ...     trueval, m, sa
+        ...     #trueval, m, sa
         ...     np.abs(m-trueval)<sa
-        (0, 1.2309597785531439e-19, 3.2562166904166163e-18)
         True
-        (array([ 1.74952003]), 1.7500502911292997, 0.022461280887884415)
         array([ True], dtype=bool)
-        (0.0, -0.00024040615507690482, 0.0087615749174770451)
         True
-        (0.0, 0.0018609047569154713, 0.049873521257196997)
         True
 
         waveplot(x1,'r',x2,'g',1,1)
@@ -1726,7 +1722,29 @@ class SpecData1D(WafoData):
         ...     res = fun(x2[:,1::], axis=0)
         ...     m = res.mean()
         ...     sa = res.std()
-        ...     trueval, m, sa
+        ...     #trueval, m, sa
+        ...     np.abs(m-trueval)<sa
+        True
+        True
+        True
+        True
+        
+        >>> x = []
+        >>> for i in range(20):
+        ...     x2, x1 = S.sim_nl(ns=20000,cases=1)
+        ...     x.append(x2[:,1::])
+        >>> x2 = np.hstack(x)
+        >>> truth1 = [0,np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk')
+        >>> truth1[-1] = truth1[-1]-3
+        >>> truth1
+        [0, 1.7495200310090628, 0.18673120577479821, 0.06198852126241805]
+         
+        >>> funs = [np.mean,np.std,st.skew,st.kurtosis]
+        >>> for fun,trueval in zip(funs,truth1):
+        ...     res = fun(x2, axis=0)
+        ...     m = res.mean()
+        ...     sa = res.std()
+        ...     #trueval, m, sa
         ...     np.abs(m-trueval)<sa
         True
         True
@@ -1889,10 +1907,20 @@ class SpecData1D(WafoData):
 ##        % 1'st order + 2'nd order component.
 ##        x2(:,2:end) =x(:,2:end)+ real(x2s(1:np,:))+real(x2d(1:np,:))
 ##        else
-        amp = np.array(amp.T).ravel()
-        rvec, ivec = c_library.disufq(amp.real, amp.imag, w, kw, water_depth, g, nmin, nmax, cases, ns)
-
-        svec = rvec + 1J * ivec
+        if False:
+            # TODO: disufq does not work for cases>1
+            amp = np.array(amp.T).ravel()
+            rvec, ivec = c_library.disufq(amp.real, amp.imag, w, kw, water_depth,
+                                      g, nmin, nmax, cases, ns)
+            svec = rvec + 1J * ivec
+        else:
+            amp = amp.T
+            svec=[]
+            for i in range(cases):
+                rvec, ivec = c_library.disufq(amp[i].real, amp[i].imag, w, kw, water_depth,
+                                      g, nmin, nmax, 1, ns)
+                svec.append(rvec + 1J * ivec)
+            svec = np.hstack(svec)
         svec.shape = (cases, ns)
         x2o = fft(svec, axis=1).T # 2'nd order component
 
@@ -3022,7 +3050,7 @@ class SpecData2D(WafoData):
         Example:
         >>> import wafo.spectrum.models as sm
         >>> D = sm.Spreading()
-        >>> SD = D.tospecdata2d(sm.Jonswap().tospecdata())
+        >>> SD = D.tospecdata2d(sm.Jonswap().tospecdata(),nt=101)
         >>> m,mtext = SD.moment(nr=2,vari='xyt')
         >>> np.round(m,3),mtext
         (array([ 3.061,  0.132, -0.   ,  2.13 ,  0.011,  0.008,  1.677, -0.   ,
diff --git a/pywafo/src/wafo/spectrum/models.py b/pywafo/src/wafo/spectrum/models.py
index 07b30cb..bd88885 100644
--- a/pywafo/src/wafo/spectrum/models.py
+++ b/pywafo/src/wafo/spectrum/models.py
@@ -447,9 +447,9 @@ class Jonswap(ModelSpectrum):
     N      : scalar defining decay of high frequency part.   (default 5)
     M      : scalar defining spectral width around the peak. (default 4)
     method : String defining method used to estimate Ag when gamma>1
-            'integrate' : Ag = 1/gaussq(Gf*ggamspec(wn,N,M),0,wnc) (default)
-            'parametric': Ag = (1+f1(N,M)*log(gamma)**f2(N,M))/gamma
-            'custom'    : Ag = Ag
+            'integration': Ag = 1/gaussq(Gf*ggamspec(wn,N,M),0,wnc) (default)
+            'parametric' : Ag = (1+f1(N,M)*log(gamma)**f2(N,M))/gamma
+            'custom'     : Ag = Ag
     wnc    : wc/wp normalized cut off frequency used when calculating Ag
                 by integration (default 6)
     Parameters