diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py
index bbfa013..e660f30 100644
--- a/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py
+++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py
@@ -64,7 +64,8 @@ dtyex = S1.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)
 dtyest = Sest.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3)
 
 T, index = ts.wave_periods(vh=0, pdef='d2u')
-wm.histgrm(T, n=25, odd=True, scale=True)
+bins = wm.good_bins(T, num_bins=25, odd=True)
+wm.plot_histgrm(T, bins=bins, scale=True)
 
 dtyex.plot()
 dtyest.plot('-.')
diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py
index 63474b7..930ec3f 100644
--- a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py
+++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py
@@ -37,7 +37,7 @@ t = linspace(0.01,8,200);
 ftc = wk.TKDE(Tc, L2=0, inc=128)
 
 plot(t,ftc.eval_grid(t), t, ftc.eval_grid_fast(t),'-.') 
-wm.histgrm(Tc,scale=True)
+wm.plot_histgrm(Tc,normed=True)
 title('Kernel Density Estimates')
 axis([0, 8, 0, 0.5])
 show() 
diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py
index 90d84d6..1209b03 100644
--- a/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py
+++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py
@@ -39,363 +39,368 @@ xx_sea = wd.sea()
 ts = wo.mat2timeseries(xx_sea)
 tp = ts.turning_points()
 mM = tp.cycle_pairs(kind='min2max')
-lc = mM.level_crossings()
-lci = lc.copy()
-T_sea = lci.args[-1]-lci.args[0]
-lci.data = lci.data/T_sea
-lci.labels.ylab='Crossing intensity'
-subplot(2,2,1)
-lci.plot()
-subplot(2,2,2)
-lci.setplotter(plotmethod='step')
+lc = mM.level_crossings(intensity=True)
+T_sea = ts.args[-1]-ts.args[0]
+
+subplot(1,2,1)
+lc.plot()
+subplot(1,2,2)
+lc.setplotter(plotmethod='step')
+lc.plot()
 show() 
  
  
 m_sea = ts.data.mean() 
-f0_sea = interp(m_sea, lci.args,lci.data)
+f0_sea = interp(m_sea, lc.args,lc.data)
 extr_sea = len(tp.data)/(2*T_sea)
 alfa_sea = f0_sea/extr_sea
 print('alfa = %g ' % alfa_sea )
 
 #! Section 4.3.2 Extraction of rainflow cycles
 #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#!#! Min-max and rainflow cycle plots
-
-mM_rfc=tp.cycle_pairs(h=0.3)
-mM_sea=tp2mm(tp_sea);
-clf
+#! Min-max and rainflow cycle plots
+#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+mM_rfc = tp.cycle_pairs(h=0.3)
+ 
+clf()
 subplot(122), 
 mM.plot() 
 title('min-max cycle count')
 subplot(121), 
 mM_rfc.plot()
 title('Rainflow cycle count')
+show()
 
-#!#! Min-max and rainflow cycle distributions
-ampmM_sea=cc2amp(mM_sea);
-ampRFC_sea=cc2amp(RFC_sea);
-clf
-subplot(221), hist(ampmM_sea,25);
+#! Min-max and rainflow cycle distributions
+#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+import wafo.misc as wm
+ampmM_sea = mM.amplitudes()
+ampRFC_sea = mM_rfc.amplitudes()
+clf()
+subplot(121) 
+wm.plot_histgrm(ampmM_sea,25)
+ylim = gca().get_ylim()
 title('min-max amplitude distribution')
-subplot(222), hist(ampRFC_sea,25);
+subplot(122)
+wm.plot_histgrm(ampRFC_sea,25)
+gca().set_ylim(ylim)
 title('Rainflow amplitude distribution')
-wafostamp([],'(ER)')
-disp('Block 4'),pause(pstate)
+show()
 
 #!#! Section 4.3.3 Simulation of rainflow cycles
 #!#! Simulation of cycles in a Markov model
-n=41; param_m=[-1 1 n]; param_D=[1 n n];
+n=41; param_m=[-1, 1, n]; param_D=[1, n, n];
 u_markov=levels(param_m);
-G_markov=mktestmat(param_m,[-0.2 0.2],0.15,1);
+G_markov=mktestmat(param_m,[-0.2, 0.2],0.15,1);
 T_markov=5000;
-xxD_markov=mctpsim({G_markov []},T_markov);
-xx_markov=[(1:T_markov)' u_markov(xxD_markov)'];
-clf
-plot(xx_markov(1:50,1),xx_markov(1:50,2))
-title('Markov chain of turning points')
-wafostamp([],'(ER)')
-disp('Block 5'),pause(pstate)
-
-
-#!#! Rainflow cycles in a transformed Gaussian model
-#!#! Hermite transformed wave data and rainflow filtered turning points, h = 0.2.
-me = mean(xx_sea(:,2));
-sa = std(xx_sea(:,2));
-Hm0_sea = 4*sa;
-Tp_sea = 1/max(lc_sea(:,2));
-spec = jonswap([],[Hm0_sea Tp_sea]);
-
-[sk, ku] = spec2skew(spec);
-spec.tr = hermitetr([],[sa sk ku me]);
-param_h = [-1.5 2 51];
-spec_norm = spec;
-spec_norm.S = spec_norm.S/sa^2;
-xx_herm = spec2sdat(spec_norm,[2^15 1],0.1);
-#! ????? PJ, JR 11-Apr-2001
-#! NOTE, in the simulation program spec2sdat
-#!the spectrum must be normalized to variance 1 
-#! ?????
-h = 0.2;
-[dtp,u_herm,xx_herm_1]=dat2dtp(param_h,xx_herm,h);
-clf
-plot(xx_herm(:,1),xx_herm(:,2),'k','LineWidth',2); hold on;
-plot(xx_herm_1(:,1),xx_herm_1(:,2),'k--','Linewidth',2);
-axis([0 50 -1 1]), hold off;
-title('Rainflow filtered wave data')
-wafostamp([],'(ER)')
-disp('Block 6'),pause(pstate)
-
-#!#! Rainflow cycles and rainflow filtered rainflow cycles in the transformed Gaussian process.
-tp_herm=dat2tp(xx_herm);
-RFC_herm=tp2rfc(tp_herm);
-mM_herm=tp2mm(tp_herm);
-h=0.2;
-[dtp,u,tp_herm_1]=dat2dtp(param_h,xx_herm,h);
-RFC_herm_1 = tp2rfc(tp_herm_1);
-clf
-subplot(121), ccplot(RFC_herm)
-title('h=0')
-subplot(122), ccplot(RFC_herm_1)
-title('h=0.2')
-if (printing==1), print -deps ../bilder/fatigue_8.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 7'),pause(pstate)
-
-#!#! Section 4.3.4 Calculating the rainflow matrix
-
-
-Grfc_markov=mctp2rfm({G_markov []});
-clf
-subplot(121), cmatplot(u_markov,u_markov,G_markov), axis('square')
-subplot(122), cmatplot(u_markov,u_markov,Grfc_markov), axis('square')
-wafostamp([],'(ER)')
-disp('Block 8'),pause(pstate)
-
-#!#! 
-clf
-cmatplot(u_markov,u_markov,{G_markov Grfc_markov},3) 
-wafostamp([],'(ER)')
-disp('Block 9'),pause(pstate)	
-
-#!#! Min-max-matrix and theoretical rainflow matrix for test Markov sequence.
-cmatplot(u_markov,u_markov,{G_markov Grfc_markov},4)
-subplot(121), axis('square'), title('min2max transition matrix')
-subplot(122), axis('square'), title('Rainflow matrix')
-if (printing==1), print -deps ../bilder/fatigue_9.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 10'),pause(pstate)
-
-#!#! Observed and theoretical rainflow matrix for test Markov sequence.
-n=length(u_markov);
-Frfc_markov=dtp2rfm(xxD_markov,n);
-clf
-cmatplot(u_markov,u_markov,{Frfc_markov Grfc_markov*T_markov/2},3) 
-subplot(121), axis('square'), title('Observed rainflow matrix')
-subplot(122), axis('square'), title('Theoretical rainflow matrix')
-if (printing==1), print -deps ../bilder/fatigue_10.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 11'),pause(pstate)
-
-#!#! Smoothed observed and calculated rainflow matrix for test Markov sequence.
-tp_markov=dat2tp(xx_markov);
-RFC_markov=tp2rfc(tp_markov);
-h=1;
-Frfc_markov_smooth=cc2cmat(param_m,RFC_markov,[],1,h);
-clf
-cmatplot(u_markov,u_markov,{Frfc_markov_smooth Grfc_markov*T_markov/2},4)
-subplot(121), axis('square'), title('Smoothed observed rainflow matrix')
-subplot(122), axis('square'), title('Theoretical rainflow matrix')
-if (printing==1), print -deps ../bilder/fatigue_11.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 12'),pause(pstate)
-
-#!#! Rainflow matrix from spectrum
-clf
-#!GmM3_herm=spec2mmtpdf(spec,[],'Mm',[],[],2);
-GmM3_herm=spec2cmat(spec,[],'Mm',[],param_h,2);
-pdfplot(GmM3_herm)
-wafostamp([],'(ER)')
-disp('Block 13'),pause(pstate)
-
-
-#!#! Min-max matrix and theoretical rainflow matrix for Hermite-transformed Gaussian waves.
-Grfc_herm=mctp2rfm({GmM3_herm.f []});
-u_herm=levels(param_h);
-clf
-cmatplot(u_herm,u_herm,{GmM3_herm.f Grfc_herm},4)
-subplot(121), axis('square'), title('min-max matrix')
-subplot(122), axis('square'), title('Theoretical rainflow matrix')
-if (printing==1), print -deps ../bilder/fatigue_12.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 14'),pause(pstate)
-
-#!#!
-clf
-Grfc_direct_herm=spec2cmat(spec,[],'rfc',[],[],2);
-subplot(121), pdfplot(GmM3_herm), axis('square'), hold on
-subplot(122), pdfplot(Grfc_direct_herm), axis('square'), hold off
-if (printing==1), print -deps ../bilder/fig_mmrfcjfr.eps
-end
-wafostamp([],'(ER)')
-disp('Block 15'),pause(pstate)
-
-
-#!#! Observed smoothed and theoretical min-max matrix, 
-#!#! (and observed smoothed and theoretical rainflow matrix for Hermite-transformed Gaussian waves).
-tp_herm=dat2tp(xx_herm);
-RFC_herm=tp2rfc(tp_herm);
-mM_herm=tp2mm(tp_herm);
-h=0.2;
-FmM_herm_smooth=cc2cmat(param_h,mM_herm,[],1,h);
-Frfc_herm_smooth=cc2cmat(param_h,RFC_herm,[],1,h);
-T_herm=xx_herm(end,1)-xx_herm(1,1);
-clf
-cmatplot(u_herm,u_herm,{FmM_herm_smooth GmM3_herm.f*length(mM_herm) ; ...
-      Frfc_herm_smooth Grfc_herm*length(RFC_herm)},4)
-subplot(221), axis('square'), title('Observed smoothed min-max matrix')
-subplot(222), axis('square'), title('Theoretical min-max matrix')
-subplot(223), axis('square'), title('Observed smoothed rainflow matrix')
-subplot(224), axis('square'), title('Theoretical rainflow matrix')
-if (printing==1), print -deps ../bilder/fatigue_13.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 16'),pause(pstate)
-   
-#!#! Section 4.3.5 Simulation from crossings and rainflow structure
-
-#!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve)
-#!#! for simulated process with irregularity factor 0.25.
-clf
-cross_herm=dat2lc(xx_herm);
-alpha1=0.25;
-alpha2=0.75;
-xx_herm_sim1=lc2sdat(cross_herm,500,alpha1);
-cross_herm_sim1=dat2lc(xx_herm_sim1);
-subplot(211)
-plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2)))
-hold on
-stairs(cross_herm_sim1(:,1),...
-    cross_herm_sim1(:,2)/max(cross_herm_sim1(:,2)))
-hold off
-title('Crossing intensity, \alpha = 0.25')
-subplot(212)
-plot(xx_herm_sim1(:,1),xx_herm_sim1(:,2))
-title('Simulated load, \alpha = 0.25')
-if (printing==1), print -deps ../bilder/fatigue_14_25.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 16'),pause(pstate)
-
-#!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve)
-#!#! for simulated process with irregularity factor 0.75.
-xx_herm_sim2=lc2sdat(cross_herm,500,alpha2);
-cross_herm_sim2=dat2lc(xx_herm_sim2);
-subplot(211)
-plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2)))
-hold on
-stairs(cross_herm_sim2(:,1),...
-    cross_herm_sim2(:,2)/max(cross_herm_sim2(:,2)))
-hold off
-title('Crossing intensity, \alpha = 0.75')
-subplot(212)
-plot(xx_herm_sim2(:,1),xx_herm_sim2(:,2))
-title('Simulated load, \alpha = 0.75')
-if (printing==1), print -deps ../bilder/fatigue_14_75.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 17'),pause(pstate)
-
-#!#! Section 4.4 Fatigue damage and fatigue life distribution
-#!#! Section 4.4.1 Introduction
-beta=3.2; gam=5.5E-10; T_sea=xx_sea(end,1)-xx_sea(1,1);
-d_beta=cc2dam(RFC_sea,beta)/T_sea;
-time_fail=1/gam/d_beta/3600      #!in hours of the specific storm
-disp('Block 18'),pause(pstate)
-
-#!#! Section 4.4.2 Level crossings
-#!#! Crossing intensity as calculated from the Markov matrix (solid curve) and from the observed rainflow matrix (dashed curve).
-clf
-mu_markov=cmat2lc(param_m,Grfc_markov);
-muObs_markov=cmat2lc(param_m,Frfc_markov/(T_markov/2));
-clf
-plot(mu_markov(:,1),mu_markov(:,2),muObs_markov(:,1),muObs_markov(:,2),'--')
-title('Theoretical and observed crossing intensity ')
-if (printing==1), print -deps ../bilder/fatigue_15.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 19'),pause(pstate)
-
-#!#! Section 4.4.3 Damage
-#!#! Distribution of damage from different RFC cycles, from calculated theoretical and from observed rainflow matrix.
-beta = 4;
-Dam_markov = cmat2dam(param_m,Grfc_markov,beta)
-DamObs1_markov = cc2dam(RFC_markov,beta)/(T_markov/2)
-DamObs2_markov = cmat2dam(param_m,Frfc_markov,beta)/(T_markov/2)
-disp('Block 20'),pause(pstate)
-
-Dmat_markov = cmat2dmat(param_m,Grfc_markov,beta);
-DmatObs_markov = cmat2dmat(param_m,Frfc_markov,beta)/(T_markov/2); 
-clf
-subplot(121), cmatplot(u_markov,u_markov,Dmat_markov,4)
-title('Theoretical damage matrix') 
-subplot(122), cmatplot(u_markov,u_markov,DmatObs_markov,4)
-title('Observed damage matrix') 
-if (printing==1), print -deps ../bilder/fatigue_16.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 21'),pause(pstate)
-
-
-#!#!
-#!Damplus_markov = lc2dplus(mu_markov,beta)
-pause(pstate)
-
-#!#! Section 4.4.4 Estimation of S-N curve
-
-#!#! Load SN-data and plot in log-log scale.
-SN = load('sn.dat');
-s = SN(:,1);
-N = SN(:,2);
-clf
-loglog(N,s,'o'), axis([0 14e5 10 30])
-#!if (printing==1), print -deps ../bilder/fatigue_?.eps end
-wafostamp([],'(ER)')
-disp('Block 22'),pause(pstate)
-
-
-#!#! Check of S-N-model on normal probability paper.
-
-normplot(reshape(log(N),8,5))
-if (printing==1), print -deps ../bilder/fatigue_17.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 23'),pause(pstate)
-
-#!#! Estimation of S-N-model on linear  scale.
-clf
-[e0,beta0,s20] = snplot(s,N,12);
-title('S-N-data with estimated N(s)','FontSize',20)
-set(gca,'FontSize',20)
-if (printing==1), print -deps ../bilder/fatigue_18a.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 24'),pause(pstate)
-
-#!#! Estimation of S-N-model on log-log scale.
-clf
-[e0,beta0,s20] = snplot(s,N,14);
-title('S-N-data with estimated N(s)','FontSize',20)
-set(gca,'FontSize',20)
-if (printing==1), print -deps ../bilder/fatigue_18b.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 25'),pause(pstate)
-
-#!#! Section 4.4.5 From S-N curve to fatigue life distribution
-#!#! Damage intensity as function of $\beta$
-beta = 3:0.1:8;
-DRFC = cc2dam(RFC_sea,beta);
-dRFC = DRFC/T_sea;
-plot(beta,dRFC), axis([3 8 0 0.25])
-title('Damage intensity as function of \beta')
-if (printing==1), print -deps ../bilder/fatigue_19.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 26'),pause(pstate)
-
-#!#! Fatigue life distribution with sea load.
-dam0 = cc2dam(RFC_sea,beta0)/T_sea;
-[t0,F0] = ftf(e0,dam0,s20,0.5,1);
-[t1,F1] = ftf(e0,dam0,s20,0,1);
-[t2,F2] = ftf(e0,dam0,s20,5,1);
-plot(t0,F0,t1,F1,t2,F2)
-title('Fatigue life distribution function')
-if (printing==1), print -deps ../bilder/fatigue_20.eps 
-end
-wafostamp([],'(ER)')
-disp('Block 27, last block')
\ No newline at end of file
+#xxD_markov=mctpsim({G_markov [,]},T_markov);
+#xx_markov=[(1:T_markov)' u_markov(xxD_markov)'];
+#clf
+#plot(xx_markov(1:50,1),xx_markov(1:50,2))
+#title('Markov chain of turning points')
+#wafostamp([],'(ER)')
+#disp('Block 5'),pause(pstate)
+#
+#
+##!#! Rainflow cycles in a transformed Gaussian model
+##!#! Hermite transformed wave data and rainflow filtered turning points, h = 0.2.
+#me = mean(xx_sea(:,2));
+#sa = std(xx_sea(:,2));
+#Hm0_sea = 4*sa;
+#Tp_sea = 1/max(lc_sea(:,2));
+#spec = jonswap([],[Hm0_sea Tp_sea]);
+#
+#[sk, ku] = spec2skew(spec);
+#spec.tr = hermitetr([],[sa sk ku me]);
+#param_h = [-1.5 2 51];
+#spec_norm = spec;
+#spec_norm.S = spec_norm.S/sa^2;
+#xx_herm = spec2sdat(spec_norm,[2^15 1],0.1);
+##! ????? PJ, JR 11-Apr-2001
+##! NOTE, in the simulation program spec2sdat
+##!the spectrum must be normalized to variance 1 
+##! ?????
+#h = 0.2;
+#[dtp,u_herm,xx_herm_1]=dat2dtp(param_h,xx_herm,h);
+#clf
+#plot(xx_herm(:,1),xx_herm(:,2),'k','LineWidth',2); hold on;
+#plot(xx_herm_1(:,1),xx_herm_1(:,2),'k--','Linewidth',2);
+#axis([0 50 -1 1]), hold off;
+#title('Rainflow filtered wave data')
+#wafostamp([],'(ER)')
+#disp('Block 6'),pause(pstate)
+#
+##!#! Rainflow cycles and rainflow filtered rainflow cycles in the transformed Gaussian process.
+#tp_herm=dat2tp(xx_herm);
+#RFC_herm=tp2rfc(tp_herm);
+#mM_herm=tp2mm(tp_herm);
+#h=0.2;
+#[dtp,u,tp_herm_1]=dat2dtp(param_h,xx_herm,h);
+#RFC_herm_1 = tp2rfc(tp_herm_1);
+#clf
+#subplot(121), ccplot(RFC_herm)
+#title('h=0')
+#subplot(122), ccplot(RFC_herm_1)
+#title('h=0.2')
+#if (printing==1), print -deps ../bilder/fatigue_8.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 7'),pause(pstate)
+#
+##!#! Section 4.3.4 Calculating the rainflow matrix
+#
+#
+#Grfc_markov=mctp2rfm({G_markov []});
+#clf
+#subplot(121), cmatplot(u_markov,u_markov,G_markov), axis('square')
+#subplot(122), cmatplot(u_markov,u_markov,Grfc_markov), axis('square')
+#wafostamp([],'(ER)')
+#disp('Block 8'),pause(pstate)
+#
+##!#! 
+#clf
+#cmatplot(u_markov,u_markov,{G_markov Grfc_markov},3) 
+#wafostamp([],'(ER)')
+#disp('Block 9'),pause(pstate)	
+#
+##!#! Min-max-matrix and theoretical rainflow matrix for test Markov sequence.
+#cmatplot(u_markov,u_markov,{G_markov Grfc_markov},4)
+#subplot(121), axis('square'), title('min2max transition matrix')
+#subplot(122), axis('square'), title('Rainflow matrix')
+#if (printing==1), print -deps ../bilder/fatigue_9.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 10'),pause(pstate)
+#
+##!#! Observed and theoretical rainflow matrix for test Markov sequence.
+#n=length(u_markov);
+#Frfc_markov=dtp2rfm(xxD_markov,n);
+#clf
+#cmatplot(u_markov,u_markov,{Frfc_markov Grfc_markov*T_markov/2},3) 
+#subplot(121), axis('square'), title('Observed rainflow matrix')
+#subplot(122), axis('square'), title('Theoretical rainflow matrix')
+#if (printing==1), print -deps ../bilder/fatigue_10.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 11'),pause(pstate)
+#
+##!#! Smoothed observed and calculated rainflow matrix for test Markov sequence.
+#tp_markov=dat2tp(xx_markov);
+#RFC_markov=tp2rfc(tp_markov);
+#h=1;
+#Frfc_markov_smooth=cc2cmat(param_m,RFC_markov,[],1,h);
+#clf
+#cmatplot(u_markov,u_markov,{Frfc_markov_smooth Grfc_markov*T_markov/2},4)
+#subplot(121), axis('square'), title('Smoothed observed rainflow matrix')
+#subplot(122), axis('square'), title('Theoretical rainflow matrix')
+#if (printing==1), print -deps ../bilder/fatigue_11.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 12'),pause(pstate)
+#
+##!#! Rainflow matrix from spectrum
+#clf
+##!GmM3_herm=spec2mmtpdf(spec,[],'Mm',[],[],2);
+#GmM3_herm=spec2cmat(spec,[],'Mm',[],param_h,2);
+#pdfplot(GmM3_herm)
+#wafostamp([],'(ER)')
+#disp('Block 13'),pause(pstate)
+#
+#
+##!#! Min-max matrix and theoretical rainflow matrix for Hermite-transformed Gaussian waves.
+#Grfc_herm=mctp2rfm({GmM3_herm.f []});
+#u_herm=levels(param_h);
+#clf
+#cmatplot(u_herm,u_herm,{GmM3_herm.f Grfc_herm},4)
+#subplot(121), axis('square'), title('min-max matrix')
+#subplot(122), axis('square'), title('Theoretical rainflow matrix')
+#if (printing==1), print -deps ../bilder/fatigue_12.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 14'),pause(pstate)
+#
+##!#!
+#clf
+#Grfc_direct_herm=spec2cmat(spec,[],'rfc',[],[],2);
+#subplot(121), pdfplot(GmM3_herm), axis('square'), hold on
+#subplot(122), pdfplot(Grfc_direct_herm), axis('square'), hold off
+#if (printing==1), print -deps ../bilder/fig_mmrfcjfr.eps
+#end
+#wafostamp([],'(ER)')
+#disp('Block 15'),pause(pstate)
+#
+#
+##!#! Observed smoothed and theoretical min-max matrix, 
+##!#! (and observed smoothed and theoretical rainflow matrix for Hermite-transformed Gaussian waves).
+#tp_herm=dat2tp(xx_herm);
+#RFC_herm=tp2rfc(tp_herm);
+#mM_herm=tp2mm(tp_herm);
+#h=0.2;
+#FmM_herm_smooth=cc2cmat(param_h,mM_herm,[],1,h);
+#Frfc_herm_smooth=cc2cmat(param_h,RFC_herm,[],1,h);
+#T_herm=xx_herm(end,1)-xx_herm(1,1);
+#clf
+#cmatplot(u_herm,u_herm,{FmM_herm_smooth GmM3_herm.f*length(mM_herm) ; ...
+#      Frfc_herm_smooth Grfc_herm*length(RFC_herm)},4)
+#subplot(221), axis('square'), title('Observed smoothed min-max matrix')
+#subplot(222), axis('square'), title('Theoretical min-max matrix')
+#subplot(223), axis('square'), title('Observed smoothed rainflow matrix')
+#subplot(224), axis('square'), title('Theoretical rainflow matrix')
+#if (printing==1), print -deps ../bilder/fatigue_13.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 16'),pause(pstate)
+#   
+##!#! Section 4.3.5 Simulation from crossings and rainflow structure
+#
+##!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve)
+##!#! for simulated process with irregularity factor 0.25.
+#clf
+#cross_herm=dat2lc(xx_herm);
+#alpha1=0.25;
+#alpha2=0.75;
+#xx_herm_sim1=lc2sdat(cross_herm,500,alpha1);
+#cross_herm_sim1=dat2lc(xx_herm_sim1);
+#subplot(211)
+#plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2)))
+#hold on
+#stairs(cross_herm_sim1(:,1),...
+#    cross_herm_sim1(:,2)/max(cross_herm_sim1(:,2)))
+#hold off
+#title('Crossing intensity, \alpha = 0.25')
+#subplot(212)
+#plot(xx_herm_sim1(:,1),xx_herm_sim1(:,2))
+#title('Simulated load, \alpha = 0.25')
+#if (printing==1), print -deps ../bilder/fatigue_14_25.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 16'),pause(pstate)
+#
+##!#! Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve)
+##!#! for simulated process with irregularity factor 0.75.
+#xx_herm_sim2=lc2sdat(cross_herm,500,alpha2);
+#cross_herm_sim2=dat2lc(xx_herm_sim2);
+#subplot(211)
+#plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2)))
+#hold on
+#stairs(cross_herm_sim2(:,1),...
+#    cross_herm_sim2(:,2)/max(cross_herm_sim2(:,2)))
+#hold off
+#title('Crossing intensity, \alpha = 0.75')
+#subplot(212)
+#plot(xx_herm_sim2(:,1),xx_herm_sim2(:,2))
+#title('Simulated load, \alpha = 0.75')
+#if (printing==1), print -deps ../bilder/fatigue_14_75.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 17'),pause(pstate)
+#
+##!#! Section 4.4 Fatigue damage and fatigue life distribution
+##!#! Section 4.4.1 Introduction
+#beta=3.2; gam=5.5E-10; T_sea=xx_sea(end,1)-xx_sea(1,1);
+#d_beta=cc2dam(RFC_sea,beta)/T_sea;
+#time_fail=1/gam/d_beta/3600      #!in hours of the specific storm
+#disp('Block 18'),pause(pstate)
+#
+##!#! Section 4.4.2 Level crossings
+##!#! Crossing intensity as calculated from the Markov matrix (solid curve) and from the observed rainflow matrix (dashed curve).
+#clf
+#mu_markov=cmat2lc(param_m,Grfc_markov);
+#muObs_markov=cmat2lc(param_m,Frfc_markov/(T_markov/2));
+#clf
+#plot(mu_markov(:,1),mu_markov(:,2),muObs_markov(:,1),muObs_markov(:,2),'--')
+#title('Theoretical and observed crossing intensity ')
+#if (printing==1), print -deps ../bilder/fatigue_15.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 19'),pause(pstate)
+#
+##!#! Section 4.4.3 Damage
+##!#! Distribution of damage from different RFC cycles, from calculated theoretical and from observed rainflow matrix.
+#beta = 4;
+#Dam_markov = cmat2dam(param_m,Grfc_markov,beta)
+#DamObs1_markov = cc2dam(RFC_markov,beta)/(T_markov/2)
+#DamObs2_markov = cmat2dam(param_m,Frfc_markov,beta)/(T_markov/2)
+#disp('Block 20'),pause(pstate)
+#
+#Dmat_markov = cmat2dmat(param_m,Grfc_markov,beta);
+#DmatObs_markov = cmat2dmat(param_m,Frfc_markov,beta)/(T_markov/2); 
+#clf
+#subplot(121), cmatplot(u_markov,u_markov,Dmat_markov,4)
+#title('Theoretical damage matrix') 
+#subplot(122), cmatplot(u_markov,u_markov,DmatObs_markov,4)
+#title('Observed damage matrix') 
+#if (printing==1), print -deps ../bilder/fatigue_16.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 21'),pause(pstate)
+#
+#
+##!#!
+##!Damplus_markov = lc2dplus(mu_markov,beta)
+#pause(pstate)
+#
+##!#! Section 4.4.4 Estimation of S-N curve
+#
+##!#! Load SN-data and plot in log-log scale.
+#SN = load('sn.dat');
+#s = SN(:,1);
+#N = SN(:,2);
+#clf
+#loglog(N,s,'o'), axis([0 14e5 10 30])
+##!if (printing==1), print -deps ../bilder/fatigue_?.eps end
+#wafostamp([],'(ER)')
+#disp('Block 22'),pause(pstate)
+#
+#
+##!#! Check of S-N-model on normal probability paper.
+#
+#normplot(reshape(log(N),8,5))
+#if (printing==1), print -deps ../bilder/fatigue_17.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 23'),pause(pstate)
+#
+##!#! Estimation of S-N-model on linear  scale.
+#clf
+#[e0,beta0,s20] = snplot(s,N,12);
+#title('S-N-data with estimated N(s)','FontSize',20)
+#set(gca,'FontSize',20)
+#if (printing==1), print -deps ../bilder/fatigue_18a.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 24'),pause(pstate)
+#
+##!#! Estimation of S-N-model on log-log scale.
+#clf
+#[e0,beta0,s20] = snplot(s,N,14);
+#title('S-N-data with estimated N(s)','FontSize',20)
+#set(gca,'FontSize',20)
+#if (printing==1), print -deps ../bilder/fatigue_18b.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 25'),pause(pstate)
+#
+##!#! Section 4.4.5 From S-N curve to fatigue life distribution
+##!#! Damage intensity as function of $\beta$
+#beta = 3:0.1:8;
+#DRFC = cc2dam(RFC_sea,beta);
+#dRFC = DRFC/T_sea;
+#plot(beta,dRFC), axis([3 8 0 0.25])
+#title('Damage intensity as function of \beta')
+#if (printing==1), print -deps ../bilder/fatigue_19.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 26'),pause(pstate)
+#
+##!#! Fatigue life distribution with sea load.
+#dam0 = cc2dam(RFC_sea,beta0)/T_sea;
+#[t0,F0] = ftf(e0,dam0,s20,0.5,1);
+#[t1,F1] = ftf(e0,dam0,s20,0,1);
+#[t2,F2] = ftf(e0,dam0,s20,5,1);
+#plot(t0,F0,t1,F1,t2,F2)
+#title('Fatigue life distribution function')
+#if (printing==1), print -deps ../bilder/fatigue_20.eps 
+#end
+#wafostamp([],'(ER)')
+#disp('Block 27, last block')
\ No newline at end of file
diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py
index 62751a8..ad77de8 100644
--- a/pywafo/src/wafo/misc.py
+++ b/pywafo/src/wafo/misc.py
@@ -12,6 +12,7 @@ from numpy import (abs, amax, any, logical_and, arange, linspace, atleast_1d,atl
                    finfo, inf, pi, interp, isnan, isscalar, zeros, ones,
                    r_, sign, unique, hstack, vstack, nonzero, where, extract)
 from scipy.special import gammaln
+from scipy.integrate import trapz, simps
 import types
 import warnings
 from wafo import plotbackend
@@ -30,7 +31,7 @@ __all__ = ['is_numlike','JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_sel
     'findoutliers', 'common_shape', 'argsreduce',
     'stirlerr', 'getshipchar', 'betaloge', 'gravity', 'nextpow2',
     'discretize', 'discretize2', 'pol2cart', 'cart2pol',  'meshgrid', 'ndgrid', 
-    'trangood', 'tranproc', 'histgrm', 'num2pistr']
+    'trangood', 'tranproc', 'plot_histgrm', 'num2pistr']
 
 
 def is_numlike(obj):
@@ -59,7 +60,7 @@ class JITImport(object):
         except:
             if self._module is None:
                 self._module = __import__(self._module_name, None, None, ['*'])
-                assert(isinstance(self._module, types.ModuleType), 'module')
+                #assert(isinstance(self._module, types.ModuleType), 'module')
                 return getattr(self._module, attr)
             else:
                 raise
@@ -1942,44 +1943,46 @@ def tranproc(x, f, x0, *xi):
                     if N > 4:
                         warnings.warn('Transformation of derivatives of order>4 not supported.')
     return y #y0,y1,y2,y3,y4
-
-def histgrm(data, n=None, odd=False, scale=False, lintype='b-'):
-    '''
-    Plot histogram
-             
+def good_bins(data=None, range=None, num_bins=None, num_data=None, odd=False, loose=True):
+    ''' Return good bins for histogram
+    
     Parameters
-    -----------
+    ----------
     data : array-like
         the data
-    n : scalar integer
-        approximate number of bins wanted  (default depending on length(data))  
+    range : (float, float)
+        minimum and maximum range of bins (default data.min(), data.max())
+    num_bins : scalar integer
+        approximate number of bins wanted  (default depending on num_data=len(data))  
     odd : bool
         placement of bins (0 or 1) (default 0)
-    scale : bool
-        argument for scaling (default 0)
-        scale = 1 yields the area 1 under the histogram
-    lintype : specify color and lintype, see PLOT for possibilities.
-    
-    Returns
+    loose : bool
+        if True add extra space to min and max
+        if False the bins are made tight to the min and max
+        
+    Example
     -------
-    binwidth = the width of each bin
-    
-     Example:
-       R=rndgumb(2,2,1,100);
-       histgrm(R,20,0,1)
-       hold on
-       x=linspace(-3,16,200);
-       plot(x,pdfgumb(x,2,2),'r')
-       hold off
+    >>> import wafo.misc as wm
+    >>> wm.good_bins(range=(0,5), num_bins=6)
+    array([-1.,  0.,  1.,  2.,  3.,  4.,  5.,  6.])
+    >>> wm.good_bins(range=(0,5), num_bins=6, loose=False)
+    array([ 0.,  1.,  2.,  3.,  4.,  5.])
+    >>> wm.good_bins(range=(0,5), num_bins=6, odd=True)
+    array([-1.5, -0.5,  0.5,  1.5,  2.5,  3.5,  4.5,  5.5,  6.5])
+    >>> wm.good_bins(range=(0,5), num_bins=6, odd=True, loose=False)
+    array([-0.5,  0.5,  1.5,  2.5,  3.5,  4.5,  5.5])
     '''
     
-    x = np.atleast_1d(data)
-    if n is None:
-        n = np.ceil(4 * np.sqrt(np.sqrt(len(x))))
+    if data:
+        x = np.atleast_1d(data)
+        num_data = len(x)
     
-    mn = x.min()
-    mx = x.max()
-    d = (mx - mn) / n * 2
+    mn, mx = range if range else (x.min(), x.max()) 
+    
+    if num_bins is None:
+        num_bins = np.ceil(4 * np.sqrt(np.sqrt(num_data)))
+    
+    d = float(mx - mn) / num_bins * 2
     e = np.floor(np.log(d) / np.log(10));
     m = np.floor(d / 10 ** e)
     if m > 5:
@@ -1988,11 +1991,65 @@ def histgrm(data, n=None, odd=False, scale=False, lintype='b-'):
         m = 2
     
     d = m * 10 ** e
-    mn = (np.floor(mn / d) - 1) * d - odd * d / 2
-    mx = (np.ceil(mx / d) + 1) * d + odd * d / 2
-    limits = np.arange(mn, mx, d)
+    mn = (np.floor(mn / d) - loose) * d - odd * d / 2
+    mx = (np.ceil(mx / d) + loose) * d + odd * d / 2
+    limits = np.arange(mn, mx+d/2, d)
+    return limits
+
+def plot_histgrm(data, bins=None, range=None, normed=False, weights=None, lintype='b-'):
+    '''
+    Plot histogram
+             
+    Parameters
+    -----------
+    data : array-like
+        the data
+    bins : int or sequence of scalars, optional
+        If an int, it defines the number of equal-width
+        bins in the given range (4 * sqrt(sqrt(len(data)), by default). 
+        If a sequence, it defines the bin edges, including the 
+        rightmost edge, allowing for non-uniform bin widths.
+    range : (float, float), optional
+        The lower and upper range of the bins.  If not provided, range
+        is simply ``(data.min(), data.max())``.  Values outside the range are
+        ignored.
+    normed : bool, optional
+        If False, the result will contain the number of samples in each bin.  
+        If True, the result is the value of the probability *density* function 
+        at the bin, normalized such that the *integral* over the range is 1. 
+    weights : array_like, optional
+        An array of weights, of the same shape as `data`.  Each value in `data`
+        only contributes its associated weight towards the bin count
+        (instead of 1).  If `normed` is True, the weights are normalized,
+        so that the integral of the density over the range remains 1
+    lintype : specify color and lintype, see PLOT for possibilities.
+    
+    Returns
+    -------
+    h : list
+        of plot-objects 
     
-    bin, limits = np.histogram(data, bins=limits, normed=scale) #, new=True)
+    Example
+    -------
+    >>> import pylab as plb
+    >>> import wafo.misc as wm
+    >>> import wafo.stats as ws
+    >>> R = ws.weibull_min.rvs(2,loc=0,scale=2, size=100)
+    >>> h0 = wm.plot_histgrm(R, 20, normed=True)
+    >>> x = linspace(-3,16,200)
+    >>> h1 = plb.plot(x,ws.weibull_min.pdf(x,2,0,2),'r')
+            
+    See also
+    --------
+    wafo.misc.good_bins
+    numpy.histogram
+    '''
+    
+    x = np.atleast_1d(data)
+    if bins is None:
+        bins = np.ceil(4 * np.sqrt(np.sqrt(len(x))))
+     
+    bin, limits = np.histogram(data, bins=bins, normed=normed, weights=weights) #, new=True)
     limits.shape = (-1, 1)
     xx = limits.repeat(3, axis=1)
     xx.shape = (-1,)
@@ -2003,10 +2060,8 @@ def histgrm(data, n=None, odd=False, scale=False, lintype='b-'):
     yy[:, 0] = 0.0 # histogram
     yy.shape = (-1,)
     yy = np.hstack((yy, 0.0))
-    plotbackend.plotbackend.plot(xx, yy, lintype, limits, limits * 0)
-    binwidth = d
-    return binwidth
-
+    return plotbackend.plotbackend.plot(xx, yy, lintype, limits, limits * 0)
+     
 def num2pistr(x, n=3):
     ''' 
     Convert a scalar to a text string in fractions of pi
@@ -2041,29 +2096,28 @@ def num2pistr(x, n=3):
         xtxt = format % x
     return xtxt     
 
-def fourier(x,t,T=None,M=None,N=None, method='trapz'):
+def fourier(data,t=None,T=None,m=None,n=None, method='trapz'):
     '''
     Returns Fourier coefficients.
-    
-      CALL:  [a,b] = fourier(t,x,T,M);
-    
-        
+ 
     Parameters
     ----------
-    x : array-like
+    data : array-like
         vector or matrix of row vectors with data points shape p x n.
     t : array-like
         vector with n values indexed from 1 to N.
     T : real scalar
         primitive period of signal, i.e., smallest period. (default T = t[-1]-t[0] 
-    M : scalar integer
+    m : scalar integer
         defines no of harmonics desired (default M = N)
-    N : scalar integer 
+    n : scalar integer 
         no of data points (default len(t))
+    method : string
+        integration method used
     
     Returns
     -------
-    a,b  = Fourier coefficients size M x P
+    a,b  = Fourier coefficients size m x p
     
     FOURIER finds the coefficients for a Fourier series representation
     of the signal x(t) (given in digital form).  It is assumed the signal
@@ -2083,15 +2137,20 @@ def fourier(x,t,T=None,M=None,N=None, method='trapz'):
      
     Example
     -------
-    >>> t = linspace(0,4*T) 
-    >>> x = sin(t);
-    >>> a, b = fourier(t, x, T=2*pi, M=5)
+    >>> import wafo.misc as wm
+    >>> import numpy as np 
+    >>> T = 2*np.pi
+    >>> t = np.linspace(0,4*T) 
+    >>> x = np.sin(t)
+    >>> a, b = wm.fourier(x, t, T=T, m=5)
+    >>> (np.round(a.ravel()), np.round(b.ravel()))
+    (array([ 0., -0.,  0., -0.,  0.]), array([ 0.,  4., -0., -0.,  0.]))
     
     See also
     --------
     fft
     '''
-    x = np.atleast_2d(x)
+    x = np.atleast_2d(data)
     p, n = x.shape
     if t is None:
         t = np.arange(n)
@@ -2099,112 +2158,48 @@ def fourier(x,t,T=None,M=None,N=None, method='trapz'):
         t = np.atleast_1d(t)
     
     n = len(t) if n is None else n
-    m = n if n is none else m
+    m = n if n is None else m
     T = t[-1]-t[0] if T is None else T
     
+    if method.startswith('trapz'):
+        intfun = trapz     
+    elif method.startswith('simp'):
+        intfun = simps
+
+    # Define the vectors for computing the Fourier coefficients
+    t.shape = (1,-1) 
+    a = zeros((m,p))
+    b = zeros((m,p))
+    a[0] = intfun(x,t, axis=-1)
      
-    
-#    switch 0
-#      case -1,
-#           % Define the vectors for computing the Fourier coefficients
-#      %
-#      a = zeros(M,P);
-#      b = zeros(M,P);
-#      a(1,:) = simpson(x);
-#    
-#      %
-#      % Compute M-1 more coefficients
-#      tmp  = 2*pi*t(:,ones(1,P))/T;
-#      % tmp =  2*pi*(0:N-1).'/(N-1); 
-#      for n1 = 1:M-1,
-#        n = n1+1;
-#        a(n,:) = simpson(x.*cos(n1*tmp));
-#        b(n,:) = simpson(x.*sin(n1*tmp));
-#      end
-#      
-#      a = 2*a/N;
-#      b = 2*b/N;
-#      
-#      case 0,
-#         %
-#      a = zeros(M,P);
-#      b = zeros(M,P);
-#      a(1,:) = trapz(t,x);
-#    
-#      %
-#      % Compute M-1 more coefficients
-#      tmp  = 2*pi*t(:,ones(1,P))/T;
-#      % tmp =  2*pi*(0:N-1).'/(N-1); 
-#      for n1 = 1:M-1,
-#        n = n1+1;
-#        a(n,:) = trapz(t,x.*cos(n1*tmp));
-#        b(n,:) = trapz(t,x.*sin(n1*tmp));
-#      end
-#      a = a/pi;
-#      b = b/pi;
-#      
-#      case 1,
-#      % Define the vectors for computing the Fourier coefficients
-#      %
-#      a = zeros(M,P);
-#      b = zeros(M,P);
-#      
-#      %
-#      % Compute the dc-level (the a(0) component).
-#      %
-#      % Note: the index has to begin with "1".
-#      %
-#      
-#      a(1,:) = sum(x);
-#    
-#      %
-#      % Compute M-1 more coefficients
-#      tmp  = 2*pi*t(:,ones(1,P))/T;
-#      % tmp =  2*pi*(0:N-1).'/(N-1); 
-#      for n1 = 1:M-1,
-#        n = n1+1;
-#        a(n,:) = sum(x.*cos(n1*tmp));
-#        b(n,:) = sum(x.*sin(n1*tmp));
-#      end
-#      a = 2*a/N;
-#      b = 2*b/N;
-#    case 2,
-#       % Define the vectors for computing the Fourier coefficients
-#      %
-#      a = zeros(M,P);
-#      b = zeros(M,P);
-#      a(1,:) = trapz(x);
-#    
-#      %
-#      % Compute M-1 more coefficients
-#      tmp  = 2*pi*t(:,ones(1,P))/T;
-#      % tmp =  2*pi*(0:N-1).'/(N-1); 
-#      for n1 = 1:M-1,
-#        n = n1+1;
-#        a(n,:) = trapz(x.*cos(n1*tmp));
-#        b(n,:) = trapz(x.*sin(n1*tmp));
-#      end
+    # Compute M-1 more coefficients
+    tmp  = 2*pi*t/T
+    #% tmp =  2*pi*(0:N-1).'/(N-1); 
+    for i in range(1,m):
+        a[i] = intfun(x*cos(i*tmp),t, axis=-1)
+        b[i] = intfun(x*sin(i*tmp),t, axis=-1)
+      
+    a = a/pi
+    b = b/pi
+     
+    # Alternative:  faster for large M, but gives different results than above.
+#    nper = diff(t([1 end]))/T; %No of periods given
+#    if nper == round(nper):
+#        N1 = n/nper
+#    else:
+#        N1 = n
 #      
-#      a = 2*a/N;
-#      b = 2*b/N;
-#    case 3
-#      % Alternative:  faster for large M, but gives different results than above.
-#      nper = diff(t([1 end]))/T; %No of periods given
-#      if nper == round(nper), 
-#        N1 = N/nper;
-#      else
-#        N1 = N;
-#      end
+#
 #    
-#      % Fourier coefficients by fft
-#      Fcof1 = 2*ifft(x(1:N1,:),[],1);
-#      Pcor = [1; exp(sqrt(-1)*(1:M-1).'*t(1))]; % correction term to get
+#    # Fourier coefficients by fft
+#    Fcof1 = 2*ifft(x(1:N1,:),[],1);
+#    Pcor = [1; exp(sqrt(-1)*(1:M-1).'*t(1))]; % correction term to get
 #                                                  % the correct integration limits
-#      Fcof = Fcof1(1:M,:).*Pcor(:,ones(1,P));
-#      a = real(Fcof(1:M,:));
-#      b = imag(Fcof(1:M,:));
-#    end
-#    return
+#    Fcof = Fcof1(1:M,:).*Pcor(:,ones(1,P));
+#    a = real(Fcof(1:M,:));
+#    b = imag(Fcof(1:M,:));
+    
+    return a, b
 
  
 
diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py
index f46aade..9f63c66 100644
--- a/pywafo/src/wafo/objects.py
+++ b/pywafo/src/wafo/objects.py
@@ -100,7 +100,8 @@ class LevelCrossings(WafoData):
                 logcmin = logcros[icmax]
                 logcros = sqrt(2 * abs(logcros - logcmin))
                 logcros[0:icmax + 1] = 2 * logcros[icmax] - logcros[0:icmax + 1]
-                p = polyfit(self.args[10:-9], logcros[10:-9], 1) #least square fit
+                ncr = 10
+                p = polyfit(self.args[ncr:-ncr], logcros[ncr:-ncr], 1) #least square fit
                 if self.sigma is None:
                     self.sigma = 1.0 / p[0] #estimated standard deviation of x
                 if self.mean is None:
@@ -449,9 +450,10 @@ class CyclePairs(WafoData):
     >>> h1 = mm.plot(marker='x')
     '''
     def __init__(self, *args, **kwds):
-        self.kind = kwds.pop('kind', 'max2min')
+        self.kind = kwds.pop('kind', 'min2max')
         self.sigma = kwds.pop('sigma', None)
         self.mean = kwds.pop('mean', None)
+        self.time = kwds.pop('time', 1)
         
         options = dict(title=self.kind + ' cycle pairs',
                             xlab='min', ylab='max',
@@ -505,8 +507,8 @@ class CyclePairs(WafoData):
         amp = abs(self.amplitudes())
         return atleast_1d([K * np.sum(amp ** betai) for betai in beta])
 
-    def level_crossings(self, kind='uM'):
-        """ Return number of upcrossings from a cycle count.
+    def level_crossings(self, kind='uM', intensity=False):
+        """ Return level crossing spectrum from a cycle count.
 
         Parameters
         ----------
@@ -516,6 +518,9 @@ class CyclePairs(WafoData):
             1,'uM' : upcrossings and maxima (default).
             2,'umM': upcrossings, minima, and maxima.
             3,'um' : upcrossings and minima.
+        intensity : bool
+            True if level crossing intensity spectrum
+            False if level crossing count spectrum
         Return
         ------
         lc : level crossing object
@@ -594,7 +599,11 @@ class CyclePairs(WafoData):
             dcount = cumsum(extr[1, 0:nx]) - extr[3, 0:nx]
         elif defnr == 3: ## This are upcrossings + minima + maxima
             dcount = cumsum(extr[1, 0:nx]) + extr[2, 0:nx]
-        return LevelCrossings(dcount, levels, mean=self.mean, sigma=self.sigma)
+        ylab='Count'
+        if intensity:
+            dcount = dcount/self.time
+            ylab = 'Intensity [count/sec]'
+        return LevelCrossings(dcount, levels, mean=self.mean, sigma=self.sigma, ylab=ylab)
 
 class TurningPoints(WafoData):
     '''
@@ -705,7 +714,7 @@ class TurningPoints(WafoData):
         """
         if h>0:
             ind = findrfc(self.data, h, method=method)
-            data = self.data(ind)
+            data = self.data[ind]
         else:
             data = self.data
         if data[0] > data[1]:
@@ -724,7 +733,8 @@ class TurningPoints(WafoData):
             kind = 'max2min'
             M = data[iM:-1:2]
             m = data[iM + 1::2]
-        return CyclePairs(M, m, kind=kind, mean=self.mean, sigma=self.sigma)
+        time = self.args[-1]-self.args[0]
+        return CyclePairs(M, m, kind=kind, mean=self.mean, sigma=self.sigma, time=time)
 
 def mat2timeseries(x):
     """