Added tutorial_scripts

pywafo/src/wafo/__init__.py

@@ -1,7 +1,8 @@
 
 from info import __doc__
 import misc
-import data
+import data
+import kdetools
 import objects
 import spectrum
 import transform

pywafo/src/wafo/doc/tutorial_scripts/chapter1.py

@@ -0,0 +1,188 @@
+from scipy import *
+from pylab import *
+
+# pyreport -o chapter1.html chapter1.py
+
+#! CHAPTER1 demonstrates some applications of WAFO
+#!================================================
+#!
+#! CHAPTER1 gives an overview through examples some of the capabilities of
+#! WAFO. WAFO is a toolbox of Matlab routines for statistical analysis and
+#! simulation of random waves and loads.
+#!
+#! The commands are edited for fast computation.
+
+
+#! Section 1.4 Some applications of WAFO
+#!---------------------------------------
+#! Section 1.4.1 Simulation from spectrum, estimation of spectrum 
+#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#! Simulation of the sea surface from spectrum
+#! The following code generates 200 seconds of data sampled with 10Hz from
+#! the Torsethaugen spectrum
+import wafo.spectrum.models as sm
+S = sm.Torsethaugen(Hm0=6, Tp=8);
+S1 = S.tospecdata()
+S1.plot()
+show()
+
+##
+import wafo.objects as wo
+xs = S1.sim(ns=2000, dt=0.1)
+ts = wo.mat2timeseries(xs)
+ts.plot_wave('-')
+show()
+
+#! Estimation of spectrum 
+#!~~~~~~~~~~~~~~~~~~~~~~~
+#! A common situation is that one wants to estimate the spectrum for wave
+#! measurements. The following code simulate 20 minutes signal sampled at 4Hz
+#! and compare the spectral estimate with the original Torsethaugen spectum.
+clf()
+Fs = 4;  
+xs = S1.sim(ns=fix(20 * 60 * Fs), dt=1. / Fs) 
+ts = wo.mat2timeseries(xs) 
+Sest = ts.tospecdata(NFFT=400)
+S1.plot()
+Sest.plot('--')
+axis([0, 3, 0, 5]) # This may depend on the simulation
+show()
+
+## Section 1.4.2 Probability distributions of wave characteristics.
+## Probability distribution of wave trough period
+# WAFO gives the possibility of computing the exact probability
+# distributions for a number of characteristics given a spectral density.
+# In the following example we study the trough period extracted from the
+# time series and compared with the theoretical density computed with exact
+# spectrum, S1, and the estimated spectrum, Sest.
+clf()
+import wafo.misc as wm
+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)
+
+dtyex.plot()
+dtyest.plot('-.')
+axis([0, 10, 0, 0.35])
+show()
+
+#! Section 1.4.3 Directional spectra
+#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#! Here are a few lines of code, which produce directional spectra 
+#! with frequency independent and frequency dependent spreading.
+clf()
+plotflag = 1
+Nt = 101;   # number of angles
+th0 = pi / 2; # primary direction of waves
+Sp = 15;   # spreading parameter
+
+D1 = sm.Spreading(type='cos', theta0=th0, method=None) # frequency independent
+D12 = sm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent
+
+#SD1 = mkdspec(S1, D1)
+#SD12 = mkdspec(S1, D12);
+#plotspec(SD1, plotflag), hold on, plotspec(SD12, plotflag, '-.'); hold off
+#wafostamp('', '(ER)')
+#disp('Block = 5'), pause(pstate)
+
+#! 3D Simulation of the sea surface 
+#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#! The simulations show that frequency dependent spreading leads to
+#! much more irregular surface so the orientation of waves is less
+#! transparent compared to the frequency independent case.
+#
+#! Frequency independent spreading
+#plotflag = 1; iseed = 1;
+#
+#Nx = 2 ^ 8;Ny = Nx;Nt = 1;dx = 0.5; dy = dx; dt = 0.25; fftdim = 2;
+#randn('state', iseed)
+#Y1 = seasim(SD1, Nx, Ny, Nt, dx, dy, dt, fftdim, plotflag);
+#wafostamp('', '(ER)')
+#axis('fill')
+#disp('Block = 6'), pause(pstate)
+#
+###
+## Frequency dependent spreading
+#randn('state', iseed)
+#Y12 = seasim(SD12, Nx, Ny, Nt, dx, dy, dt, fftdim, plotflag);
+#wafostamp('', '(ER)')
+#axis('fill')
+#disp('Block = 7'), pause(pstate)
+#
+#! Estimation of directional spectrum
+#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#!  The figure is not shown in the Tutorial
+#
+# Nx = 3; Ny = 2; Nt = 2 ^ 12; dx = 10; dy = 10;dt = 0.5;
+# F = seasim(SD12, Nx, Ny, Nt, dx, dy, dt, 1, 0);  
+# Z = permute(F.Z, [3 1 2]);
+# [X, Y] = meshgrid(F.x, F.y);
+# N = Nx * Ny;
+# types = repmat(sensortypeid('n'), N, 1);
+# bfs = ones(N, 1);
+# pos = [X(:), Y(:), zeros(N, 1)];
+# h = inf;
+# nfft = 128;
+# nt = 101;
+# SDe = dat2dspec([F.t Z(:, :)], [pos types, bfs], h, nfft, nt);
+#plotspec(SDe), hold on
+#plotspec(SD12, '--'), hold off
+#disp('Block = 8'), pause(pstate)
+
+#! Section 1.4.4 Fatigue, Load cycles and Markov models.
+#! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#! Switching Markow chain of turningpoints 
+#! In fatigue applications the exact sample path is not important, but
+#! only the tops and bottoms of the load, called the sequence of turning
+#! points (TP). From the turning points one can extract load cycles, from
+#! which damage calculations and fatigue life predictions can be
+#! performed.
+#!
+#! The matlab commands below computes the intensity of rainflowcycles for 
+#! the Gaussian model with spectrum S1 using the Markow approximation. 
+#! The rainflow cycles found in the simulated load signal are shown in the 
+# figure.
+#clf
+#paramu = [-6 6 61];
+#frfc = spec2cmat(S1, [], 'rfc', [], paramu);
+#pdfplot(frfc);
+#hold on
+#tp = dat2tp(xs);
+#rfc = tp2rfc(tp);
+#plot(rfc(:, 2), rfc(:, 1), '.')
+#wafostamp('', '(ER)')
+#hold off
+#disp('Block = 9'), pause(pstate)
+
+#! Section 1.4.5 Extreme value statistics
+#!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Plot of yura87 data
+clf()
+import wafo.data as wd
+xn = wd.yura87()
+#xn = load('yura87.dat'); 
+subplot(211) 
+plot(xn[::30, 0] / 3600, xn[::30, 1], '.')
+title('Water level')
+ylabel('(m)')
+
+#! Formation of 5 min maxima
+yura = xn[:85500, 1]
+yura = np.reshape(yura, (300, 285))
+maxyura = yura.max(axis=0)
+subplot(212)
+plot(xn[299:85500:300, 0] / 3600, maxyura, '.')
+xlabel('Time (h)')
+ylabel('(m)')
+title('Maximum 5 min water level')
+show()
+
+#! Estimation of GEV for yuramax
+clf()
+import wafo.stats as ws
+phat = ws.genextreme.fit2(maxyura, method='mps')
+phat.plotfitsummary()
+show()
+#disp('Block = 11, Last block')

pywafo/src/wafo/doc/tutorial_scripts/chapter2.py

@@ -0,0 +1,356 @@
+import numpy as np
+from scipy import *
+from pylab import *
+
+# pyreport -o chapter1.html chapter1.py
+
+#! CHAPTER2 Modelling random loads and stochastic waves
+#!=======================================================
+#!
+#! Chapter2 contains the commands used in Chapter 2 of the tutorial and
+#! present some tools for analysis of random functions with
+#! respect to their correlation, spectral and distributional properties.
+#! The presentation is divided into three examples: 
+#!
+#! Example1 is devoted to estimation of different parameters in the model.
+#! Example2 deals with spectral densities and
+#! Example3 presents the use of WAFO to simulate samples of a Gaussian
+#! process.
+#!
+#! Some of the commands are edited for fast computation. 
+#! Each set of commands is followed by a 'pause' command.
+#! 
+
+#!
+#! Tested on Matlab 5.3, 7.01
+#! History
+#! Revised by Georg Lindgren sept 2009 for WAFO ver 2.5 on Matlab 7.1
+#! Revised pab sept2005
+#! Added sections -> easier to evaluate using cellmode evaluation.
+#! Revised pab Dec2004
+#! Created by GL July 13, 2000
+#! from commands used in Chapter 2
+#!
+
+pstate =  'off';
+
+#! Section 2.1 Introduction and preliminary analysis
+#!====================================================
+#! Example 1: Sea data
+#!----------------------
+#! Observed crossings compared to the expected for Gaussian signals
+
+import wafo
+import wafo.objects as wo
+xx = wafo.data.sea()
+me = xx[:,1].mean()
+sa = xx[:,1].std()
+xx[:,1] -=  me
+ts = wo.mat2timeseries(xx)
+tp = ts.turning_points()
+
+
+cc = tp.cycle_pairs()
+lc = cc.level_crossings()
+lc.plot()
+show()
+
+#! Average number of upcrossings per time unit
+#!----------------------------------------------
+#! Next we compute the mean frequency as the average number of upcrossings 
+#! per time unit of the mean level (= 0); this may require interpolation in the 
+#! crossing intensity curve, as follows.  
+T = xx[:,0].max()-xx[:,0].min()
+f0 = np.interp(0,lc.args,lc.data,0)/T  #! zero up-crossing frequency 
+
+#! Turningpoints and irregularity factor
+#!----------------------------------------
+
+fm = len(tp.data)/(2*T)           #! frequency of maxima
+alfa = f0/fm                     #! approx Tm24/Tm02
+ 
+#! Visually examine data
+#!------------------------
+#! We finish this section with some remarks about the quality
+#! of the measured data. Especially sea surface measurements can be
+#! of poor quality. We shall now check the  quality of the dataset {\tt xx}. 
+#! It is always good practice to visually examine the data 
+#! before the analysis to get an impression of the quality, 
+#! non-linearities and narrow-bandedness of the data.
+#! First we shall plot the data and zoom in on a specific region. 
+#! A part of sea data is visualized with the following commands
+clf()
+ts.plot_wave('k-',tp,'*',nfig=1, nsub=1)
+
+axis([0, 2, -2, 2])
+show()
+
+#! Finding possible spurious points
+#!------------------------------------
+#! However, if the amount of data is too large for visual examinations one
+#! could use the following criteria to find possible spurious points. One
+#! must be careful using the criteria for extremevalue analysis, because
+#! it might remove extreme waves that are OK and not spurious.
+
+import wafo.misc as wm
+dt = ts.sampling_period()
+# dt = np.diff(xx[:2,0])
+dcrit = 5*dt
+ddcrit = 9.81/2*dt*dt
+zcrit = 0
+inds, indg = wm.findoutliers(xx[:,1],zcrit,dcrit,ddcrit, verbose=True)
+
+#! Section 2.2 Frequency Modeling of Load Histories
+#!----------------------------------------------------
+#! Periodogram: Raw spectrum
+#!
+clf()
+Lmax = 9500
+S = ts.tospecdata(NFFT=Lmax)
+S.plot()
+axis([0, 5, 0, 0.7])
+show()
+
+#! Calculate moments  
+#!-------------------
+mom, text= S.moment(nr=4)
+[sa, sqrt(mom[0])]
+
+#! Section 2.2.1 Random functions in Spectral Domain - Gaussian processes
+#!--------------------------------------------------------------------------
+#! Smoothing of spectral estimate 
+#1----------------------------------
+#! By decreasing Lmax the spectrum estimate becomes smoother.
+
+clf()
+Lmax0 = 200; Lmax1 = 50
+S1 = ts.tospecdata(NFFT=Lmax0)
+S2 = ts.tospecdata(NFFT=Lmax1)
+S1.plot('-.')
+S2.plot()
+show()
+
+#! Estimated autocovariance
+#!----------------------------
+#! Obviously knowing the spectrum one can compute the covariance
+#! function. The following matlab code will compute the covariance for the 
+#! unimodal spectral density S1 and compare it with estimated 
+#! covariance of the signal xx.
+clf()
+Lmax = 80
+R1 = S1.tocovdata(nr=1)   
+Rest = ts.tocovdata(lag=Lmax)
+R1.plot('.')
+Rest.plot()
+show()
+
+#! We can see in Figure below that the covariance function corresponding to 
+#! the spectral density S2 significantly differs from the one estimated 
+#! directly from data. 
+#! It can be seen in Figure above that the covariance corresponding to S1 
+#! agrees much better with the estimated covariance function
+
+clf()
+R2 = S2.tocovdata(nr=1)
+R2.plot('.')
+Rest.plot()
+show()
+
+#! Section 2.2.2 Transformed Gaussian models
+#!-------------------------------------------
+#! We begin with computing skewness and kurtosis
+#! for the data set xx and compare it with the second order wave approximation
+#! proposed by Winterstein:
+import wafo.stats as ws
+rho3 = ws.skew(xx[:,1])
+rho4 = ws.kurtosis(xx[:,1])
+
+sk, ku=S1.stats_nl(moments='sk' )
+ 
+#! Comparisons of 3 transformations
+clf()
+import wafo.transform.models as wtm
+gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku ).trdata()
+g = wtm.TrLinear(mean=me, sigma=sa ).trdata()
+
+
+#! Linear transformation
+ 
+glc, gemp = lc.trdata()
+g.plot('r')
+glc.plot('b-') #! Transf. estimated from level-crossings
+gh.plot('b-.')  #! Hermite Transf. estimated from moments
+show()
+ 
+#!  Test Gaussianity of a stochastic process.
+#!---------------------------------------------
+#! TESTGAUSSIAN simulates  e(g(u)-u) = int (g(u)-u)^2 du  for Gaussian processes 
+#!  given the spectral density, S. The result is plotted if test0 is given.
+#!  This is useful for testing if the process X(t) is Gaussian.
+#!  If 95#! of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5#! level.
+#! 
+#! As we see from the figure below: none of the simulated values of test1 is 
+#! above 1.00. Thus the data significantly departs from a Gaussian distribution. 
+clf()
+test0 = glc.dist2gauss()
+#! the following test takes time
+N = len(xx)
+test1 = S1.testgaussian(ns=N,cases=50,t0=test0)
+sum(test1>test0)<5
+show()
+
+#! Normalplot of data xx
+#!------------------------
+#! indicates that the underlying distribution has a "heavy" upper tail and a
+#! "light" lower tail. 
+clf()
+import pylab
+ws.probplot(ts.data, dist='norm', plot=pylab)
+show()
+#! Section 2.2.3 Spectral densities of sea data
+#!-----------------------------------------------
+#! Example 2: Different forms of spectra
+#!
+import wafo.spectrum.models as wsm
+clf()
+Hm0 = 7; Tp = 11;
+spec = wsm.Jonswap(Hm0=Hm0, Tp=Tp).tospecdata()
+spec.plot()
+show()
+
+#! Directional spectrum and Encountered directional spectrum
+#! Directional spectrum
+clf()
+D = wsm.Spreading('cos2s')
+Sd = D.tospecdata2d(S)
+Sd.plot()
+show()
+
+
+##!Encountered directional spectrum
+##!--------------------------------- 
+#clf()
+#Se = spec2spec(Sd,'encdir',0,10);
+#plotspec(Se), hold on
+#plotspec(Sd,1,'--'), hold off
+##!wafostamp('','(ER)')
+#disp('Block = 17'),pause(pstate)
+#
+##!#! Frequency spectra
+#clf
+#Sd1 =spec2spec(Sd,'freq');
+#Sd2 = spec2spec(Se,'enc');
+#plotspec(spec), hold on
+#plotspec(Sd1,1,'.'),
+#plotspec(Sd2),
+##!wafostamp('','(ER)')
+#hold off
+#disp('Block = 18'),pause(pstate)
+#
+##!#! Wave number spectrum
+#clf
+#Sk = spec2spec(spec,'k1d')
+#Skd = spec2spec(Sd,'k1d')
+#plotspec(Sk), hold on
+#plotspec(Skd,1,'--'), hold off
+##!wafostamp('','(ER)')
+#disp('Block = 19'),pause(pstate)
+#
+##!#! Effect of waterdepth on spectrum
+#clf
+#plotspec(spec,1,'--'), hold on
+#S20 = spec;
+#S20.S = S20.S.*phi1(S20.w,20);
+#S20.h = 20;
+#plotspec(S20),  hold off
+##!wafostamp('','(ER)')
+#disp('Block = 20'),pause(pstate)
+#
+##!#! Section 2.3 Simulation of transformed Gaussian process
+##!#! Example 3: Simulation of random sea    
+##! The reconstruct function replaces the spurious points of seasurface by
+##! simulated data on the basis of the remaining data and a transformed Gaussian
+##! process. As noted previously one must be careful using the criteria 
+##! for finding spurious points when reconstructing a dataset, because
+##! these criteria might remove the highest and steepest waves as we can see
+##! in this plot where the spurious points is indicated with a '+' sign:
+##!
+#clf
+#[y, grec] = reconstruct(xx,inds);
+#waveplot(y,'-',xx(inds,:),'+',1,1)
+#axis([0 inf -inf inf])
+##!wafostamp('','(ER)')
+#disp('Block = 21'),pause(pstate)
+#
+##! Compare transformation (grec) from reconstructed (y) 
+##! with original (glc) from (xx)
+#clf
+#trplot(g), hold on
+#plot(gemp(:,1),gemp(:,2))
+#plot(glc(:,1),glc(:,2),'-.')
+#plot(grec(:,1),grec(:,2)), hold off 
+#disp('Block = 22'),pause(pstate)
+#
+##!#!
+#clf
+#L = 200;
+#x = dat2gaus(y,grec);
+#Sx = dat2spec(x,L);
+#disp('Block = 23'),pause(pstate)
+#      
+##!#!
+#clf
+#dt = spec2dt(Sx)
+#Ny = fix(2*60/dt) #! = 2 minutes
+#Sx.tr = grec;
+#ysim = spec2sdat(Sx,Ny);
+#waveplot(ysim,'-')
+##!wafostamp('','(CR)')
+#disp('Block = 24'),pause(pstate)
+# 
+##!#! Estimated spectrum compared to Torsethaugen spectrum
+#clf
+#Tp = 1.1;
+#H0 = 4*sqrt(spec2mom(S1,1))
+#St = torsethaugen([0:0.01:5],[H0  2*pi/Tp]);
+#plotspec(S1)
+#hold on
+#plotspec(St,'-.')
+#axis([0 6 0 0.4])
+##!wafostamp('','(ER)')
+#disp('Block = 25'),pause(pstate)
+#
+##!#!
+#clf
+#Snorm = St;
+#Snorm.S = Snorm.S/sa^2;
+#dt = spec2dt(Snorm)
+#disp('Block = 26'),pause(pstate)
+#
+##!#!
+#clf
+#[Sk Su] = spec2skew(St);
+#sa = sqrt(spec2mom(St,1));
+#gh = hermitetr([],[sa sk ku me]);
+#Snorm.tr = gh;
+#disp('Block = 27'),pause(pstate)
+#
+##!#! Transformed Gaussian model compared to Gaussian model
+#clf
+#dt = 0.5;
+#ysim_t = spec2sdat(Snorm,240,dt);
+#xsim_t = dat2gaus(ysim_t,Snorm.tr);
+#disp('Block = 28'),pause(pstate)
+#
+##!#! Compare
+##! In order to compare the Gaussian and non-Gaussian models we need to scale  
+##! \verb+xsim_t+ #!{\tt xsim$_t$} 
+##! to have the same first spectral moment as 
+##! \verb+ysim_t+,  #!{\tt ysim$_t$},  Since the process xsim_t has variance one
+##! which will be done by the following commands. 
+#clf
+#xsim_t(:,2) = sa*xsim_t(:,2);
+#waveplot(xsim_t,ysim_t,5,1,sa,4.5,'r.','b')
+##!wafostamp('','(CR)')
+#disp('Block = 29, Last block'),pause(pstate)
+

pywafo/src/wafo/doc/tutorial_scripts/chapter3.py

@@ -0,0 +1,611 @@
+%% CHAPTER3  Demonstrates distributions of wave characteristics
+%
+% Chapter3 contains the commands used in Chapter3 in the tutorial.
+% 
+% Some of the commands are edited for fast computation. 
+% Each set of commands is followed by a 'pause' command.
+% 
+
+% Tested on Matlab 5.3
+% History
+% Revised pab sept2005
+%  Added sections -> easier to evaluate using cellmode evaluation.
+% Revised by pab Feb 2005
+% -updated calls to kdetools+spec2XXpdf programs
+% Created by GL July 12, 2000
+% from commands used in Chapter 3, written by IR
+%
+
+
+%% Section 3.2 Estimation of wave characteristics from data
+%% Example 1
+pstate = 'off';
+speed = 'fast'
+%speed = 'slow'
+
+xx = load('sea.dat');
+xx(:,2) = detrend(xx(:,2));
+rate = 8;
+Tcrcr = dat2wa(xx,0,'c2c','tw',rate);
+Tc = dat2wa(xx,0,'u2d','tw',rate);
+disp('Block = 1'), pause(pstate)
+
+%% Histogram of crestperiod compared to the kernel density estimate
+clf
+mean(Tc)
+max(Tc)
+t = linspace(0.01,8,200);
+kopt = kdeoptset('L2',0);
+tic
+ftc1 = kde(Tc,kopt,t);
+toc
+pdfplot(ftc1), hold on
+histgrm(Tc,[],[],1)
+axis([0 8 0 0.5])
+wafostamp([],'(ER)')
+disp('Block = 2'), pause(pstate)
+%%
+tic
+kopt.inc = 128;
+ftc2 = kdebin(Tc,kopt);
+toc
+pdfplot(ftc2,'-.')
+title('Kernel Density Estimates')
+hold off
+disp('Block = 3'), pause(pstate)
+
+%% Extreme waves - model check: the highest and steepest wave
+clf
+method = 0;
+rate = 8;
+[S, H, Ac, At, Tcf, Tcb, z_ind, yn] = ...
+   dat2steep(xx,rate,method);
+disp('Block = 4'), pause(pstate)
+
+clf
+[Smax indS]=max(S)
+[Amax indA]=max(Ac)
+spwaveplot(yn,[indA indS],'k.')
+wafostamp([],'(ER)')
+disp('Block = 5'), pause(pstate)
+
+%% Does the highest wave contradict a transformed Gaussian model?
+clf
+inds1 = (5965:5974)'; % points to remove
+Nsim = 10;
+[y1, grec1, g2, test, tobs, mu1o, mu1oStd] = ...
+   reconstruct(xx,inds1,Nsim);
+spwaveplot(y1,indA-10)
+hold on
+plot(xx(inds1,1),xx(inds1,2),'+')
+lamb = 2.;
+muLstd = tranproc(mu1o-lamb*mu1oStd,fliplr(grec1));
+muUstd = tranproc(mu1o+lamb*mu1oStd,fliplr(grec1));
+plot (y1(inds1,1), [muLstd muUstd],'b-')
+axis([1482 1498 -1 3]), 
+wafostamp([],'(ER)')
+disp('Block = 6'), 
+pause(pstate)
+
+%% Expected value (solid) compared to data removed
+clf
+plot(xx(inds1,1),xx(inds1,2),'+'), hold on
+mu = tranproc(mu1o,fliplr(grec1));
+plot(y1(inds1,1), mu), hold off
+disp('Block = 7'), pause(pstate)
+
+%% Crest height PDF
+% Transform data so that kde works better
+clf
+L2 = 0.6;
+plotnorm(Ac.^L2)
+
+%%
+%
+fac = kde(Ac,{'L2',L2},linspace(0.01,3,200));
+pdfplot(fac)
+wafostamp([],'(ER)')
+simpson(fac.x{1},fac.f)
+disp('Block = 8'), pause(pstate)
+
+%% Empirical crest height CDF
+clf
+Fac = flipud(cumtrapz(fac.x{1},flipud(fac.f)));
+Fac = [fac.x{1} 1-Fac];
+Femp = plotedf(Ac,Fac);
+axis([0 2 0 1])
+wafostamp([],'(ER)')
+disp('Block = 9'), pause(pstate)
+
+%% Empirical crest height CDF compared to a Transformed Rayleigh approximation
+facr = trraylpdf(fac.x{1},'Ac',grec1);
+Facr = cumtrapz(facr.x{1},facr.f);
+hold on
+plot(facr.x{1},Facr,'.')
+axis([1.25 2.25 0.95 1])
+wafostamp([],'(ER)')
+disp('Block = 10'), pause(pstate)
+
+%% Joint pdf of crest period and crest amplitude
+clf
+kopt2 = kdeoptset('L2',0.5,'inc',256);
+Tc = Tcf+Tcb;
+fTcAc = kdebin([Tc Ac],kopt2);
+fTcAc.labx={'Tc [s]'  'Ac [m]'}
+pdfplot(fTcAc)
+hold on
+plot(Tc,Ac,'k.')
+hold off
+wafostamp([],'(ER)')
+disp('Block = 11'), pause(pstate)
+
+%% Example 4:  Simple wave characteristics obtained from Jonswap spectrum 
+clf
+S = jonswap([],[5 10]);
+[m,  mt]= spec2mom(S,4,[],0);
+disp('Block = 12'), pause(pstate)
+
+clf
+spec2bw(S)
+[ch Sa2] = spec2char(S,[1  3])
+disp('Block = 13'), pause(pstate)
+
+%% Section 3.3.2 Explicit form approximations of wave characteristic densities
+%% Longuett-Higgins model for Tc and Ac
+clf
+t = linspace(0,15,100);
+h = linspace(0,6,100);
+flh = lh83pdf(t,h,[m(1),m(2),m(3)]);
+disp('Block = 14'), pause(pstate)
+
+%% Transformed Longuett-Higgins model for Tc and Ac
+clf
+[sk, ku ]=spec2skew(S);
+sa = sqrt(m(1));
+gh = hermitetr([],[sa sk ku 0]);
+flhg = lh83pdf(t,h,[m(1),m(2),m(3)],gh);
+disp('Block = 15'), pause(pstate)
+
+%% Cavanie model for Tc and Ac
+clf
+t = linspace(0,10,100);
+h = linspace(0,7,100);
+fcav = cav76pdf(t,h,[m(1) m(2) m(3) m(5)],[]);
+disp('Block = 16'), pause(pstate)
+
+%% Example 5 Transformed Rayleigh approximation of crest- vs trough- amplitude
+clf
+xx = load('sea.dat');
+x = xx;
+x(:,2) = detrend(x(:,2));
+SS = dat2spec2(x);
+[sk, ku, me, si ] = spec2skew(SS);
+gh = hermitetr([],[si sk ku me]);
+Hs = 4*si;
+r = (0:0.05:1.1*Hs)';
+fac_h = trraylpdf(r,'Ac',gh);
+fat_h = trraylpdf(r,'At',gh);
+h = (0:0.05:1.7*Hs)';
+facat_h = trraylpdf(h,'AcAt',gh);
+pdfplot(fac_h)
+hold on
+pdfplot(fat_h,'--')
+hold off
+wafostamp([],'(ER)')
+disp('Block = 17'), pause(pstate)
+
+%%
+clf
+TC = dat2tc(xx, me);
+tc = tp2mm(TC);
+Ac = tc(:,2);
+At = -tc(:,1);
+AcAt = Ac+At;
+disp('Block = 18'), pause(pstate)
+
+%%
+clf
+Fac_h = [fac_h.x{1} cumtrapz(fac_h.x{1},fac_h.f)];
+subplot(3,1,1)
+Fac = plotedf(Ac,Fac_h);
+hold on
+plot(r,1-exp(-8*r.^2/Hs^2),'.')
+axis([1. 2. 0.9 1])
+title('Ac CDF')
+
+Fat_h = [fat_h.x{1} cumtrapz(fat_h.x{1},fat_h.f)];
+subplot(3,1,2)
+Fat = plotedf(At,Fat_h);
+hold on
+plot(r,1-exp(-8*r.^2/Hs^2),'.')
+axis([1. 2. 0.9 1])
+title('At CDF')
+
+Facat_h = [facat_h.x{1} cumtrapz(facat_h.x{1},facat_h.f)];
+subplot(3,1,3)
+Facat = plotedf(AcAt,Facat_h);
+hold on
+plot(r,1-exp(-2*r.^2/Hs^2),'.')
+axis([1.5 3.5 0.9 1])
+title('At+Ac CDF')
+
+wafostamp([],'(ER)')
+disp('Block = 19'), pause(pstate)
+
+%% Section 3.4 Exact wave distributions in transformed Gaussian Sea
+%% Section 3.4.1 Density of crest period, crest length or encountered crest period
+clf
+S1 = torsethaugen([],[6 8],1);
+D1 = spreading(101,'cos',pi/2,[15],[],0);
+D12 = spreading(101,'cos',0,[15],S1.w,1);
+SD1 = mkdspec(S1,D1);
+SD12 = mkdspec(S1,D12);
+disp('Block = 20'), pause(pstate)
+
+%% Crest period
+clf
+tic 
+f_tc = spec2tpdf(S1,[],'Tc',[0 11 56],[],4);
+toc
+pdfplot(f_tc)
+wafostamp([],'(ER)')
+simpson(f_tc.x{1},f_tc.f)
+disp('Block = 21'), pause(pstate)
+
+%% Crest length
+
+if strncmpi(speed,'slow',1)
+  opt1 = rindoptset('speed',5,'method',3);
+  opt2 = rindoptset('speed',5,'nit',2,'method',0);
+else
+  % fast
+  opt1 = rindoptset('speed',7,'method',3);
+  opt2 = rindoptset('speed',7,'nit',2,'method',0);
+end
+
+
+clf
+if strncmpi(speed,'slow',1)
+  NITa = 5;
+else
+  disp('NIT=5 may take time, running with NIT=3 in the following')
+  NITa = 3;
+end
+%f_Lc = spec2tpdf2(S1,[],'Lc',[0 200 81],opt1); % Faster and more accurate
+f_Lc = spec2tpdf(S1,[],'Lc',[0 200 81],[],NITa);
+pdfplot(f_Lc,'-.')
+wafostamp([],'(ER)')
+disp('Block = 22'), pause(pstate)
+
+
+f_Lc_1 = spec2tpdf(S1,[],'Lc',[0 200 81],1.5,NITa);
+%f_Lc_1 = spec2tpdf2(S1,[],'Lc',[0 200 81],1.5,opt1);
+
+hold on
+pdfplot(f_Lc_1)
+wafostamp([],'(ER)')
+
+disp('Block = 23'), pause(pstate)
+%%
+clf
+simpson(f_Lc.x{1},f_Lc.f)
+simpson(f_Lc_1.x{1},f_Lc_1.f)
+      
+disp('Block = 24'), pause(pstate)
+%%
+clf
+tic
+
+f_Lc_d1 = spec2tpdf(rotspec(SD1,pi/2),[],'Lc',[0 300 121],[],NITa);
+f_Lc_d12 = spec2tpdf(SD12,[],'Lc',[0 200 81],[],NITa);
+% f_Lc_d1 = spec2tpdf2(rotspec(SD1,pi/2),[],'Lc',[0 300 121],opt1);
+% f_Lc_d12 = spec2tpdf2(SD12,[],'Lc',[0 200 81],opt1);
+toc
+pdfplot(f_Lc_d1,'-.'), hold on
+pdfplot(f_Lc_d12),     hold off
+wafostamp([],'(ER)')
+
+disp('Block = 25'), pause(pstate)
+
+%%
+
+
+clf
+opt1 = rindoptset('speed',5,'method',3);
+SD1r = rotspec(SD1,pi/2);
+if strncmpi(speed,'slow',1)
+  f_Lc_d1_5 = spec2tpdf(SD1r,[], 'Lc',[0 300 121],[],5);
+  pdfplot(f_Lc_d1_5),     hold on
+else
+  % fast
+  disp('Run the following example only if you want a check on computing time')
+  disp('Edit the command file and remove %')
+end
+f_Lc_d1_3 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],3);
+f_Lc_d1_2 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],2);
+f_Lc_d1_0 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],0);
+%f_Lc_d1_n4 = spec2tpdf2(SD1r,[],'Lc',[0 400 161],opt1);
+
+pdfplot(f_Lc_d1_3), hold on
+pdfplot(f_Lc_d1_2)
+pdfplot(f_Lc_d1_0)
+%pdfplot(f_Lc_d1_n4)
+
+%simpson(f_Lc_d1_n4.x{1},f_Lc_d1_n4.f)
+
+disp('Block = 26'), pause(pstate)
+
+%% Section 3.4.2 Density of wave period, wave length or encountered wave period
+%% Example 7: Crest period and high crest waves
+clf
+tic
+xx = load('sea.dat');
+x = xx;
+x(:,2) = detrend(x(:,2));
+SS = dat2spec(x);
+si = sqrt(spec2mom(SS,1));
+SS.tr = dat2tr(x);
+Hs = 4*si
+method = 0;
+rate = 2;
+[S, H, Ac, At, Tcf, Tcb, z_ind, yn] = dat2steep(x,rate,method);
+Tc = Tcf+Tcb;
+t = linspace(0.01,8,200);
+ftc1 = kde(Tc,{'L2',0},t);
+pdfplot(ftc1)
+hold on
+%      f_t = spec2tpdf(SS,[],'Tc',[0 8 81],0,4);
+f_t = spec2tpdf(SS,[],'Tc',[0 8 81],0,2);
+simpson(f_t.x{1},f_t.f)
+pdfplot(f_t,'-.')
+hold off
+wafostamp([],'(ER)')
+toc
+disp('Block = 27'), pause(pstate)
+
+%%
+clf
+tic
+
+if strncmpi(speed,'slow',1)
+  NIT = 4;
+else
+  NIT = 2;
+end
+%      f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],[Hs/2],4);
+tic
+f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],Hs/2,NIT);
+toc
+
+Pemp = sum(Ac>Hs/2)/sum(Ac>0)
+simpson(f_t2.x{1},f_t2.f)
+index = find(Ac>Hs/2);
+ftc1 = kde(Tc(index),{'L2',0},t);
+ftc1.f = Pemp*ftc1.f;
+pdfplot(ftc1)
+hold on
+pdfplot(f_t2,'-.')
+hold off
+wafostamp([],'(ER)')
+toc
+disp('Block = 28'), pause(pstate)
+
+%% Example 8: Wave period for high crest waves 
+%      clf
+      tic 
+      f_tcc2 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],-1);
+toc
+      simpson(f_tcc2.x{1},f_tcc2.f)
+      f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],3,5);
+%      f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],1,5);
+      simpson(f_tcc3.x{1},f_tcc3.f)
+      pdfplot(f_tcc2,'-.')
+      hold on
+      pdfplot(f_tcc3)
+      hold off
+       toc
+disp('Block = 29'), pause(pstate)
+
+%%
+clf
+[TC tc_ind v_ind] = dat2tc(yn,[],'dw');
+N = length(tc_ind);
+t_ind = tc_ind(1:2:N);
+c_ind = tc_ind(2:2:N);
+Pemp = sum(yn(t_ind,2)<-Hs/2 & yn(c_ind,2)>Hs/2)/length(t_ind)
+ind = find(yn(t_ind,2)<-Hs/2 & yn(c_ind,2)>Hs/2);
+spwaveplot(yn,ind(2:4))
+wafostamp([],'(ER)')
+disp('Block = 30'), pause(pstate)
+
+%%
+clf
+Tcc = yn(v_ind(1+2*ind),1)-yn(v_ind(1+2*(ind-1)),1);
+t = linspace(0.01,14,200);
+ftcc1 = kde(Tcc,{'kernel' 'epan','L2',0},t);
+ftcc1.f = Pemp*ftcc1.f;
+pdfplot(ftcc1,'-.')
+wafostamp([],'(ER)')
+disp('Block = 31'), pause(pstate)
+
+tic 
+f_tcc22_1 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[Hs/2],-1);
+toc
+simpson(f_tcc22_1.x{1},f_tcc22_1.f)
+hold on
+pdfplot(f_tcc22_1)
+hold off
+wafostamp([],'(ER)')
+disp('Block = 32'), pause(pstate)
+
+disp('The rest of this chapter deals with joint densities.')
+disp('Some calculations may take some time.') 
+disp('You could experiment with other NIT.')
+%return
+
+%% Section 3.4.3 Joint density of crest period and crest height
+%% Example 9. Some preliminary analysis of the data
+clf
+tic
+yy = load('gfaksr89.dat');
+SS = dat2spec(yy);
+si = sqrt(spec2mom(SS,1));
+SS.tr = dat2tr(yy);
+Hs = 4*si
+v = gaus2dat([0 0],SS.tr);
+v = v(2)
+toc
+disp('Block = 33'), pause(pstate)
+
+%%
+clf
+tic
+[TC, tc_ind, v_ind] = dat2tc(yy,v,'dw');
+N       = length(tc_ind);
+t_ind   = tc_ind(1:2:N);
+c_ind   = tc_ind(2:2:N);
+v_ind_d = v_ind(1:2:N+1);
+v_ind_u = v_ind(2:2:N+1);
+T_d     = ecross(yy(:,1),yy(:,2),v_ind_d,v);
+T_u     = ecross(yy(:,1),yy(:,2),v_ind_u,v);
+
+Tc = T_d(2:end)-T_u(1:end);
+Tt = T_u(1:end)-T_d(1:end-1);
+Tcf = yy(c_ind,1)-T_u;
+Ac = yy(c_ind,2)-v;
+At = v-yy(t_ind,2);
+toc
+disp('Block = 34'), pause(pstate)
+
+%%
+clf
+tic
+t = linspace(0.01,15,200);
+kopt3 = kdeoptset('hs',0.25,'L2',0); 
+ftc1 = kde(Tc,kopt3,t);
+ftt1 = kde(Tt,kopt3,t);
+pdfplot(ftt1,'k')
+hold on
+pdfplot(ftc1,'k-.')
+f_tc4 = spec2tpdf(SS,[],'Tc',[0 12 81],0,4,5);
+f_tc2 = spec2tpdf(SS,[],'Tc',[0 12 81],0,2,5);
+f_tc = spec2tpdf(SS,[],'Tc',[0 12 81],0,-1);
+pdfplot(f_tc,'b')
+hold off
+legend('kde(Tt)','kde(Tc)','f_{tc}')
+wafostamp([],'(ER)')
+toc
+disp('Block = 35'), pause(pstate)
+
+%% Example 10: Joint characteristics of a half wave:
+%% position and height of a crest for a wave with given period
+clf
+tic
+ind = find(4.4<Tc & Tc<4.6);
+f_AcTcf = kde([Tcf(ind) Ac(ind)],{'L2',[1 .5]});
+pdfplot(f_AcTcf)
+hold on
+plot(Tcf(ind), Ac(ind),'.');
+wafostamp([],'(ER)')
+toc
+disp('Block = 36'), pause(pstate)
+
+%%
+clf
+tic
+opt1 = rindoptset('speed',5,'method',3);
+opt2 = rindoptset('speed',5,'nit',2,'method',0);
+
+f_tcfac1 = spec2thpdf(SS,[],'TcfAc',[4.5 4.5 46],[0:0.25:8],opt1);
+f_tcfac2 = spec2thpdf(SS,[],'TcfAc',[4.5 4.5 46],[0:0.25:8],opt2);
+
+pdfplot(f_tcfac1,'-.')
+hold on
+pdfplot(f_tcfac2)
+plot(Tcf(ind), Ac(ind),'.');
+
+simpson(f_tcfac1.x{1},simpson(f_tcfac1.x{2},f_tcfac1.f,1))
+simpson(f_tcfac2.x{1},simpson(f_tcfac2.x{2},f_tcfac2.f,1))
+f_tcf4=spec2tpdf(SS,[],'Tc',[4.5 4.5 46],[0:0.25:8],6);
+f_tcf4.f(46)
+toc
+wafostamp([],'(ER)')
+disp('Block = 37'), pause(pstate)
+
+%%
+clf
+f_tcac_s = spec2thpdf(SS,[],'TcAc',[0 12 81],[Hs/2:0.1:2*Hs],opt1);
+disp('Block = 38'), pause(pstate)
+
+clf
+tic
+mom = spec2mom(SS,4,[],0);
+t = f_tcac_s.x{1};
+h = f_tcac_s.x{2};
+flh_g = lh83pdf(t',h',[mom(1),mom(2),mom(3)],SS.tr);
+clf
+ind=find(Ac>Hs/2);
+plot(Tc(ind), Ac(ind),'.');
+hold on
+pdfplot(flh_g,'k-.')
+pdfplot(f_tcac_s)
+toc
+wafostamp([],'(ER)')
+disp('Block = 39'), pause(pstate)
+
+%%
+clf
+%      f_tcac = spec2thpdf(SS,[],'TcAc',[0 12 81],[0:0.2:8],opt1);
+%      pdfplot(f_tcac)
+disp('Block = 40'), pause(pstate)
+
+%% Section 3.4.4 Joint density of crest and trough height
+%% Section 3.4.5 Min-to-max distributions – Markov method
+%% Example 11. (min-max problems with Gullfaks data)
+%% Joint density of maximum and the following minimum
+clf
+tic
+tp = dat2tp(yy);
+Mm = fliplr(tp2mm(tp));
+fmm = kde(Mm);
+f_mM = spec2mmtpdf(SS,[],'mm',[],[-7 7 51],opt2);
+
+pdfplot(f_mM,'-.')
+hold on
+pdfplot(fmm,'k-')
+hold off
+wafostamp([],'(ER)')
+toc
+disp('Block = 41'), pause(pstate)
+
+%% The joint density of ”still water separated”  maxima and minima.
+clf
+tic
+ind = find(Mm(:,1)>v & Mm(:,2)<v);
+Mmv = abs(Mm(ind,:)-v);
+fmmv = kde(Mmv);
+f_vmm = spec2mmtpdf(SS,[],'vmm',[],[-7 7 51],opt2);
+clf
+pdfplot(fmmv,'k-')
+hold on
+pdfplot(f_vmm,'-.')
+hold off
+wafostamp([],'(ER)')
+toc
+disp('Block = 42'), pause(pstate)
+
+
+%%
+clf
+tic
+facat = kde([Ac At]);
+f_acat = spec2mmtpdf(SS,[],'AcAt',[],[-7 7 51],opt2);
+clf
+pdfplot(f_acat,'-.')
+hold on
+pdfplot(facat,'k-')
+hold off
+wafostamp([],'(ER)')
+toc
+disp('Block = 43'), pause(pstate)
+

pywafo/src/wafo/doc/tutorial_scripts/chapter4.py

@@ -0,0 +1,388 @@
+%% CHAPTER4 contains the commands used in Chapter 4 of the tutorial
+%
+% CALL:  Chapter4
+% 
+% Some of the commands are edited for fast computation. 
+% Each set of commands is followed by a 'pause' command.
+% 
+% This routine also can print the figures; 
+% For printing the figures on directory ../bilder/ edit the file and put  
+%     printing=1;
+
+% Tested on Matlab 5.3
+% History
+% Revised pab sept2005
+%  Added sections -> easier to evaluate using cellmode evaluation.
+% revised pab Feb2004
+% updated call to lc2sdat
+% Created by GL July 13, 2000
+% from commands used in Chapter 4
+%
+ 
+%% Chapter 4 Fatigue load analysis and rain-flow cycles
+
+pstate = 'off';
+
+printing=0;
+%set(0,'DefaultAxesFontSize',15')
+
+%% Section 4.3.1 Crossing intensity
+xx_sea = load('sea.dat'); 
+tp_sea = dat2tp(xx_sea);
+lc_sea = tp2lc(tp_sea);
+T_sea = xx_sea(end,1)-xx_sea(1,1);
+lc_sea(:,2) = lc_sea(:,2)/T_sea;
+clf
+subplot(221), plot(lc_sea(:,1),lc_sea(:,2)) 
+title('Crossing intensity, (u, \mu(u))')
+subplot(222), semilogx(lc_sea(:,2),lc_sea(:,1)) 
+title('Crossing intensity, (log \mu(u), u)')
+wafostamp([],'(ER)')
+disp('Block 1'), pause(pstate)
+
+m_sea = mean(xx_sea(:,2));
+f0_sea = interp1(lc_sea(:,1),lc_sea(:,2),m_sea,'linear')
+extr_sea = length(tp_sea)/(2*T_sea);
+alfa_sea = f0_sea/extr_sea
+disp('Block 2'),pause(pstate)
+
+%% Section 4.3.2 Extraction of rainflow cycles
+%% Min-max and rainflow cycle plots
+RFC_sea=tp2rfc(tp_sea);
+mM_sea=tp2mm(tp_sea);
+clf
+subplot(122), ccplot(mM_sea);
+title('min-max cycle count')
+subplot(121), ccplot(RFC_sea);
+title('Rainflow cycle count')
+wafostamp([],'(ER)')
+disp('Block 3'),pause(pstate)
+
+%% Min-max and rainflow cycle distributions
+ampmM_sea=cc2amp(mM_sea);
+ampRFC_sea=cc2amp(RFC_sea);
+clf
+subplot(221), hist(ampmM_sea,25);
+title('min-max amplitude distribution')
+subplot(222), hist(ampRFC_sea,25);
+title('Rainflow amplitude distribution')
+wafostamp([],'(ER)')
+disp('Block 4'),pause(pstate)
+
+%% 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];
+u_markov=levels(param_m);
+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')

pywafo/src/wafo/doc/tutorial_scripts/chapter5.py

@@ -0,0 +1,195 @@
+%% CHAPTER5 contains the commands used in Chapter 5 of the tutorial
+%
+% CALL:  Chapter5
+% 
+% Some of the commands are edited for fast computation. 
+% Each set of commands is followed by a 'pause' command.
+% 
+
+% Tested on Matlab 5.3
+% History
+% Added Return values by GL August 2008
+% Revised pab sept2005
+% Added sections -> easier to evaluate using cellmode evaluation.
+% Created by GL July 13, 2000
+% from commands used in Chapter 5
+%
+ 
+%% Chapter 5 Extreme value analysis
+
+%% Section 5.1 Weibull and Gumbel papers
+pstate = 'off'
+
+% Significant wave-height data on Weibull paper,
+clf
+Hs = load('atlantic.dat');
+wei = plotweib(Hs)
+wafostamp([],'(ER)')
+disp('Block = 1'),pause(pstate)
+
+%%
+% Significant wave-height data on Gumbel paper,
+clf
+gum=plotgumb(Hs)
+wafostamp([],'(ER)')
+  disp('Block = 2'),pause(pstate)
+
+%%
+% Significant wave-height data on Normal probability paper,
+plotnorm(log(Hs),1,0);
+wafostamp([],'(ER)')
+disp('Block = 3'),pause(pstate)
+
+%%
+% Return values in the Gumbel distribution
+clf
+T=1:100000;
+sT=gum(2) - gum(1)*log(-log(1-1./T));
+semilogx(T,sT), hold
+N=1:length(Hs); Nmax=max(N);
+plot(Nmax./N,sort(Hs,'descend'),'.')
+title('Return values in the Gumbel model')
+xlabel('Return period')
+ylabel('Return value') 
+wafostamp([],'(ER)')
+disp('Block = 4'),pause(pstate)
+
+%% Section 5.2 Generalized Pareto and Extreme Value distributions
+%% Section 5.2.1 Generalized Extreme Value distribution
+
+% Empirical distribution of significant wave-height with estimated 
+% Generalized Extreme Value distribution,
+gev=fitgev(Hs,'plotflag',2)
+wafostamp([],'(ER)')
+disp('Block = 5a'),pause(pstate)
+
+clf
+x = linspace(0,14,200);
+plotkde(Hs,[x;pdfgev(x,gev)]')
+disp('Block = 5b'),pause(pstate)
+
+% Analysis of yura87 wave data. 
+% Wave data interpolated (spline) and organized in 5-minute intervals
+% Normalized to mean 0 and std = 1 to get stationary conditions. 
+% maximum level over each 5-minute interval analysed by GEV
+xn  = load('yura87.dat');
+XI  = 0:0.25:length(xn);
+N   = length(XI); N = N-mod(N,4*60*5); 
+YI  = interp1(xn(:,1),xn(:,2),XI(1:N),'spline');
+YI  = reshape(YI,4*60*5,N/(4*60*5)); % Each column holds 5 minutes of 
+                                     % interpolated data.
+Y5  = (YI-ones(1200,1)*mean(YI))./(ones(1200,1)*std(YI)); 
+Y5M = max(Y5);
+Y5gev = fitgev(Y5M,'method','mps','plotflag',2)
+wafostamp([],'(ER)')
+disp('Block = 6'),pause(pstate)
+
+%% Section 5.2.2 Generalized Pareto distribution
+
+% Exceedances of significant wave-height data over level 3,
+gpd3 = fitgenpar(Hs(Hs>3)-3,'plotflag',1);
+wafostamp([],'(ER)')
+
+%%
+figure
+% Exceedances of significant wave-height data over level 7,
+gpd7 = fitgenpar(Hs(Hs>7),'fixpar',[nan,nan,7],'plotflag',1);
+wafostamp([],'(ER)')
+disp('Block = 6'),pause(pstate)
+
+%% 
+%Simulates 100 values from the GEV distribution with parameters (0.3, 1, 2), then estimates the
+%parameters using two different methods and plots the estimated distribution functions together
+%with the empirical distribution.
+Rgev = rndgev(0.3,1,2,1,100);
+gp = fitgev(Rgev,'method','pwm');
+gm = fitgev(Rgev,'method','ml','start',gp.params,'plotflag',0);
+
+x=sort(Rgev);
+plotedf(Rgev,gp,{'-','r-'});
+hold on
+plot(x,cdfgev(x,gm),'--')
+hold off
+wafostamp([],'(ER)')
+disp('Block =7'),pause(pstate)
+
+%%
+% Similarly for the GPD distribution;
+Rgpd = rndgenpar(0.4,1,0,1,100);
+plotedf(Rgpd);
+hold on
+gp = fitgenpar(Rgpd,'method','pkd','plotflag',0);
+x=sort(Rgpd);
+plot(x,cdfgenpar(x,gp))
+% gm = fitgenpar(Rgpd,'method','mom','plotflag',0);
+% plot(x,cdfgenpar(x,gm),'g--')
+gw = fitgenpar(Rgpd,'method','pwm','plotflag',0);
+plot(x,cdfgenpar(x,gw),'g:')
+gml = fitgenpar(Rgpd,'method','ml','plotflag',0);
+plot(x,cdfgenpar(x,gml),'--')
+gmps = fitgenpar(Rgpd,'method','mps','plotflag',0);
+plot(x,cdfgenpar(x,gmps),'r-.')
+hold off
+wafostamp([],'(ER)')
+disp('Block = 8'),pause(pstate)
+
+%%
+% Return values for the GEV distribution
+T = logspace(1,5,10);
+[sT, sTlo, sTup] = invgev(1./T,Y5gev,'lowertail',false,'proflog',true);
+
+%T = 2:100000;
+%k=Y5gev.params(1); mu=Y5gev.params(3); sigma=Y5gev.params(2);
+%sT1 = invgev(1./T,Y5gev,'lowertail',false);
+%sT=mu + sigma/k*(1-(-log(1-1./T)).^k);
+clf
+semilogx(T,sT,T,sTlo,'r',T,sTup,'r'), hold
+N=1:length(Y5M); Nmax=max(N);
+plot(Nmax./N,sort(Y5M,'descend'),'.')
+title('Return values in the GEV model')
+xlabel('Return priod')
+ylabel('Return value') 
+grid on 
+disp('Block = 9'),pause(pstate)
+
+%% Section 5.3 POT-analysis
+
+% Estimated expected exceedance over level u as function of u.
+clf
+plotreslife(Hs,'umin',2,'umax',10,'Nu',200);
+wafostamp([],'(ER)')
+disp('Block = 10'),pause(pstate)
+
+%%
+% Estimated distribution functions of monthly maxima 
+%with the POT method (solid),
+% fitting a GEV (dashed) and the empirical distribution.
+
+% POT- method
+gpd7 = fitgenpar(Hs(Hs>7)-7,'method','pwm','plotflag',0);
+khat = gpd7.params(1);
+sigmahat = gpd7.params(2);
+muhat = length(Hs(Hs>7))/(7*3*2);
+bhat = sigmahat/muhat^khat;
+ahat = 7-(bhat-sigmahat)/khat;
+x = linspace(5,15,200);
+plot(x,cdfgev(x,khat,bhat,ahat))
+disp('Block = 11'),pause(pstate)
+
+%%
+% Since we have data to compute the monthly maxima mm over 
+%42 months we can also try to fit a
+% GEV distribution directly:
+mm = zeros(1,41);
+for i=1:41
+  mm(i)=max(Hs(((i-1)*14+1):i*14)); 
+end
+
+gev=fitgev(mm);
+hold on
+plotedf(mm)
+plot(x,cdfgev(x,gev),'--')
+hold off
+wafostamp([],'(ER)')
+disp('Block = 12, Last block'),pause(pstate)
+

pywafo/src/wafo/transform/models.py

@@ -22,7 +22,7 @@ from numpy import (sqrt, atleast_1d, abs, imag, sign, where, cos, arccos, ceil,
     expm1, log1p, pi) #@UnresolvedImport
 import numpy as np
 import warnings
-from core import TrCommon
+from core import TrCommon, TrData
 __all__ = ['TrHermite', 'TrLinear', 'TrOchi']
 
 _example = '''
@@ -41,10 +41,40 @@ _example = '''
     >>> g2 = <generic>(mean=me, var=va, skew=sk, kurt=ku, ysigma=std)
     >>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world 
     '''
+class TrCommon2(TrCommon):
+    __doc__ = TrCommon.__doc__
+    def trdata(self,x=None, xnmin=-5, xnmax=5, n=513):
+        """
+        Return a discretized transformation model.
+
+        Parameters
+        ----------
+        x : vector  (default sigma*linspace(xnmin,xnmax,n)+mean)
+        xnmin : real, scalar
+            minimum on normalized scale
+        xnmax : real, scalar
+            maximum on normalized scale
+        n : integer, scalar
+            number of evaluation points
+
+        Returns
+        -------
+        t0 : real, scalar
+            a measure of departure from the Gaussian model calculated as
+            trapz(xn,(xn-g(x))**2.) where int. limits is given by X.
+        """
+        if x is None:
+            xn = np.linspace(xnmin, xnmax, n)
+            x = self.sigma*xn+self.mean
+        else:
+            xn = (x-self.mean)/self.sigma
 
+        yn = (self._dat2gauss(x)-self.ymean)/self.ysigma
+        
+        return TrData(yn, x, mean=self.mean, sigma=self.sigma)
 
-class TrHermite(TrCommon):
+class TrHermite(TrCommon2):
-    __doc__ = TrCommon.__doc__.replace('<generic>', 'Hermite') + """
+    __doc__ = TrCommon2.__doc__.replace('<generic>', 'Hermite') + """
     pardef : scalar, integer
         1  Winterstein et. al. (1994) parametrization [1]_ (default)
         2  Winterstein (1988) parametrization [2]_
@@ -294,8 +324,8 @@ class TrHermite(TrCommon):
 
 
 
-class TrLinear(TrCommon):
+class TrLinear(TrCommon2):
-    __doc__ = TrCommon.__doc__.replace('<generic>', 'Linear') + """
+    __doc__ = TrCommon2.__doc__.replace('<generic>', 'Linear') + """
     Description
     -----------
     The linear transformation model is monotonic linear polynomial, calibrated
@@ -333,8 +363,8 @@ class TrLinear(TrCommon):
         return x
     
    
-class TrOchi(TrCommon):
+class TrOchi(TrCommon2):
-    __doc__ = TrCommon.__doc__.replace('<generic>', 'Ochi') + """
+    __doc__ = TrCommon2.__doc__.replace('<generic>', 'Ochi') + """
 
     Description
     -----------