From 5dbc448368144c145cb13f8f12320c5a9dbd45c8 Mon Sep 17 00:00:00 2001 From: "per.andreas.brodtkorb" Date: Tue, 18 Jan 2011 06:54:23 +0000 Subject: [PATCH] Added tutorial_scripts --- pywafo/src/wafo/__init__.py | 3 +- pywafo/src/wafo/doc/__init__.py | 0 .../src/wafo/doc/tutorial_scripts/chapter1.py | 188 ++++++ .../src/wafo/doc/tutorial_scripts/chapter2.py | 356 ++++++++++ .../src/wafo/doc/tutorial_scripts/chapter3.py | 611 ++++++++++++++++++ .../src/wafo/doc/tutorial_scripts/chapter4.py | 388 +++++++++++ .../src/wafo/doc/tutorial_scripts/chapter5.py | 195 ++++++ pywafo/src/wafo/transform/models.py | 44 +- 8 files changed, 1777 insertions(+), 8 deletions(-) create mode 100644 pywafo/src/wafo/doc/__init__.py create mode 100644 pywafo/src/wafo/doc/tutorial_scripts/chapter1.py create mode 100644 pywafo/src/wafo/doc/tutorial_scripts/chapter2.py create mode 100644 pywafo/src/wafo/doc/tutorial_scripts/chapter3.py create mode 100644 pywafo/src/wafo/doc/tutorial_scripts/chapter4.py create mode 100644 pywafo/src/wafo/doc/tutorial_scripts/chapter5.py diff --git a/pywafo/src/wafo/__init__.py b/pywafo/src/wafo/__init__.py index e54ef8c..af89cb9 100644 --- a/pywafo/src/wafo/__init__.py +++ b/pywafo/src/wafo/__init__.py @@ -1,7 +1,8 @@ from info import __doc__ import misc -import data +import data +import kdetools import objects import spectrum import transform diff --git a/pywafo/src/wafo/doc/__init__.py b/pywafo/src/wafo/doc/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter1.py new file mode 100644 index 0000000..afd2d19 --- /dev/null +++ b/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') diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py new file mode 100644 index 0000000..b93be3a --- /dev/null +++ b/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) + diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py new file mode 100644 index 0000000..34dd4ff --- /dev/null +++ b/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.4Hs/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) 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') \ No newline at end of file diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter5.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter5.py new file mode 100644 index 0000000..2583ac9 --- /dev/null +++ b/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) + diff --git a/pywafo/src/wafo/transform/models.py b/pywafo/src/wafo/transform/models.py index 7384abd..15d0b7d 100644 --- a/pywafo/src/wafo/transform/models.py +++ b/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 = (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): - __doc__ = TrCommon.__doc__.replace('', 'Hermite') + """ +class TrHermite(TrCommon2): + __doc__ = TrCommon2.__doc__.replace('', '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): - __doc__ = TrCommon.__doc__.replace('', 'Linear') + """ +class TrLinear(TrCommon2): + __doc__ = TrCommon2.__doc__.replace('', 'Linear') + """ Description ----------- The linear transformation model is monotonic linear polynomial, calibrated @@ -333,8 +363,8 @@ class TrLinear(TrCommon): return x -class TrOchi(TrCommon): - __doc__ = TrCommon.__doc__.replace('', 'Ochi') + """ +class TrOchi(TrCommon2): + __doc__ = TrCommon2.__doc__.replace('', 'Ochi') + """ Description -----------