Home > wafo > trgauss > spec2cmat.m

# spec2cmat

## PURPOSE

Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough)

## SYNOPSIS

[f, fmm] = spec2cmat(spec,utc,def,paramt,paramu,nit)

## DESCRIPTION

``` SPEC2CMAT Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough)

CALL:  f   = spec2cmat(S,u,def,paramt,paramu,nit);

f    = pdf (density structure) of crests (trough) heights
S    = spectral density structure
u    = reference level (default the most frequently crossed level).
def   = 'Mm'  : gives maximum and the following minimum height.
'rfc' : gives maximum and the rainflow minimum height.
'AcAt': gives (crest,trough) heights  (this option needs
more work).
paramt  = [0 tn Nt] defines discretization of half period: tn is
the longest period considered while Nt is the number of
points, i.e. (Nt-1)/tn is the sampling frequnecy.
paramt=[0 10 51] implies that the halfperiods are
considered at 51 linearly spaced points in the interval
[0,10], i.e. sampling frequency is 5 Hz.
paramu  = [u v N] defines discretization of maxima and minima ranges:
u is the lowest minimum considered, v the heighest
maximum and N is the number of levles (u,v) included.
nit  =  0,...,9. Dimension of numerical integration (only
positive nit are allowed). (default nit=1).
[]    = default values are used.

The model for loads is a stationary Gaussian transformed process X(t),
where  Y(t) = g(X(t)) is a zero-mean Gaussian with spectrum, S.

Note: algorithm uses Markov Chain approximation to the sequence of
turning points in Y.

Example: % The intensity matrix of rainflow cycles is computed by:
S  = jonswap;
L0 = spec2mom(S,1);
paramu = [sqrt(L0)*[-4 4] 41];
frfc   = spec2cmat(S,[],'rfc',[],paramu);

## CROSS-REFERENCE INFORMATION

This function calls:
 createpdf PDF class constructor dat2gaus Transforms x using the transformation g. gaus2dat Transforms xx using the inverse of transformation g. levels Calculates discrete levels given the parameter matrix. mctp2rfc Rainflow matrix given a Markov matrix of a Markov chain of turning points mctp2tc Frequencies of upcrossing troughs and crests using Markov chain of turning points. qlevels Calculates quantile levels which encloses P% of PDF tranproc Transforms process X and up to four derivatives wminmax Calculates joint density of minimum and following maximum wnormspec Normalize a spectral density such that m0=m2=1 clock Current date and time as date vector. error Display message and abort function. etime Elapsed time. isfield True if field is in structure array. lower Convert string to lowercase. num2str Convert number to string. (Fast version) warning Display warning message; disable or enable warning messages.
This function is called by:
 Chapter1 % CHAPTER1 demonstrates some applications of WAFO Chapter4 % CHAPTER4 contains the commands used in Chapter 4 of the tutorial itmkurs_lab3 Script to computer exercises 3

## SOURCE CODE

```001 function [f, fmm] = spec2cmat(spec,utc,def,paramt,paramu,nit)
002 %SPEC2CMAT Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough)
003 %
004 % CALL:  f   = spec2cmat(S,u,def,paramt,paramu,nit);
005 %
006 %        f    = pdf (density structure) of crests (trough) heights
007 %        S    = spectral density structure
008 %        u    = reference level (default the most frequently crossed level).
009 %       def   = 'Mm'  : gives maximum and the following minimum height.
010 %               'rfc' : gives maximum and the rainflow minimum height.
011 %               'AcAt': gives (crest,trough) heights  (this option needs
012 %                       more work).
013 %     paramt  = [0 tn Nt] defines discretization of half period: tn is
014 %               the longest period considered while Nt is the number of
015 %               points, i.e. (Nt-1)/tn is the sampling frequnecy.
016 %               paramt=[0 10 51] implies that the halfperiods are
017 %               considered at 51 linearly spaced points in the interval
018 %               [0,10], i.e. sampling frequency is 5 Hz.
019 %     paramu  = [u v N] defines discretization of maxima and minima ranges:
020 %               u is the lowest minimum considered, v the heighest
021 %               maximum and N is the number of levles (u,v) included.
022 %        nit  =  0,...,9. Dimension of numerical integration (only
023 %                positive nit are allowed). (default nit=1).
024 %       []    = default values are used.
025 %
026 %
027 %  The model for loads is a stationary Gaussian transformed process X(t),
028 %  where  Y(t) = g(X(t)) is a zero-mean Gaussian with spectrum, S.
029 %
030 %  Note: algorithm uses Markov Chain approximation to the sequence of
031 %  turning points in Y.
032 %
033 % Example: % The intensity matrix of rainflow cycles is computed by:
034 %       S  = jonswap;
035 %       L0 = spec2mom(S,1);
036 %       paramu = [sqrt(L0)*[-4 4] 41];
037 %       frfc   = spec2cmat(S,[],'rfc',[],paramu);
038 %
040
041
042 % Tested on : matlab 5.3
043 % History: by I. Rychlik 01.10.1998 with name minmax.m
044 % bounds by I.R. 02.01.2000.
045 % Revised by es 000322. Made call with directional spectrum possible.
046 % revised by ir 000612. Help and plots improved.
047 % revised by IR removing error in transformation 29 VI 2000
048 % revised by I.R. 01.20.2001 Change name minmax to wminmax
049 % revised by I.R. 6 II 2001 adapted for MATLAB 6
050 % revised pab 30nov2003
051
052 % TODO % AcAt option needs more work
053 startTime = clock;
054 [S, xl4, L0, L2, L4, L1]=wnormspec(spec);
055
056
057 A = sqrt(L0/L2);
058 SCIS=0;
059 if nargin<6|isempty(nit)
060   nit=1;
061 elseif nit<0
062   warning('Only postive nit allowed')
063  nit=1;
064 end
065
066 if isfield(spec,'tr')
067    g = spec.tr;
068 else
069    g = [];
070 end
071 if isempty(g)
072   g = [sqrt(L0)*(-5:0.02:5)', (-5:0.02:5)'];
073 end
074 S.tr = g;
075
076 if nargin<7|isempty(speed)
077    speed=5;
078 end
079 if nargin<2|isempty(utc)
080     utc_d=gaus2dat([0, 0],g); % most frequent crossed level
081     utc=utc_d(1,2);
082 end
083
084 % transform reference level into Gaussian level
085 uu = dat2gaus([0., utc],g);
086 u  = uu(2);
087 disp(['The level u for Gaussian process = ', num2str(u)])
088
089
090 if nargin<4|isempty(paramt)
091   distanceBetweenExtremes = 5*pi*sqrt(L2/L4); %(2.5 * mean distance between extremes)
092   paramt = [0 distanceBetweenExtremes,41];
093 end
094 t0     = paramt(1);
095 tn     = paramt(2);
096 Ntime  = paramt(3);
097 t      = levels([0, tn/A, Ntime]); % normalized times
098
099
100 if nargin<3|isempty(def)
101     defnr=0;
102 else
103  switch lower(def)
104  case  'mm',    defnr = 1;
105  case  'rfc',   defnr = 0;
106  case  'acat',  defnr =-1;
107  otherwise, error('Unknown def')
108  end
109 end
110
111
112 nr = 4;
113
114 if nargin<5|isempty(paramu)
115    paramu=[-4*sqrt(L0), 4*sqrt(L0), 41];
116 end
117  %Transform amplitudes to Gaussian levels:
118 h   = levels(paramu);
119 h   = reshape(h,length(h),1);
120 Nx  = length(h);
121 der = ones(Nx,1);
122 hg  = tranproc([h, der],g);
123 der = abs(hg(:,2));
124 hg  = hg(:,1); % Gaussian level
125
126 %if exist('h.in'), delete h.in, end
127 %fid=fopen('h.in','wt');
128 %fprintf(fid,'%12.10f\n',hg);
129 %fclose(fid);
130
131
132
133 %paru=paramu;
134 % paru(1:2)=paru(1:2)/sqrt(L0);
135  fmM = wminmax(S,nit,paramu,t);
136
137
138 Htxt=['Joint density of crest and trough'];
139 labx='crest [m]';
140 laby='trough [m]';
141 if (defnr==0)
142    Htxt = ['Joint density of maximum and rainflow minimum'];
143    labx='max [m]';
144    laby='rainflow min [m]';
145 end
146 f=createpdf;
147 f.title=Htxt;
148 f.nit=nit;
149 f.speed=speed;
150 f.SCIS=SCIS;
151 f.labx{1}=labx;
152 f.labx{2}=laby;
153 f.x{1}=h;
154 fmm=createpdf;
155 fmm.title = ['Joint density of maximum and minimum'];
156 fmm.labx{1}='max [m]';
157 fmm.labx{2}='min [m]';
158 fmm.nit=nit;
159 fmm.speed=speed;
160 fmm.SCIS=SCIS;
161 dh=h(2)-h(1);
162 f.x{1}=h;
163 f.x{2}=h;
164 fmm.x{1}=h;
165 fmm.x{2}=h;
166 ftmp=fmM;
167
168 for i=1:Nx
169   ftmp(:,i)=dh*dh*ftmp(:,i);%.*der*der(i);%* sqrt(-R(1,6)/R(1,4))/2/pi;
170 end
171 fmm.f=ftmp;
172
173 if (defnr==0)
174   f.f=fliplr(mctp2rfc(fliplr(ftmp)));%* sqrt(-R(1,6)/R(1,4))/2/pi;
175   fmm.f=ftmp;
176 end
177 if (defnr==1)
178   f.f=ftmp;
179 end
180 if (defnr==-1)
181   f.f    = fliplr(mctp2tc(fliplr(ftmp),utc,paramu));
182   index1 = find(f.x{1}>0);
183   index2 = find(f.x{2}<0);
184   f.f    = flipud(f.f(index2,index1));
185   f.x{1} = f.x{1}(index1);
186   f.x{2} = abs(flipud(f.x{2}(index2)));
187 end
188 [f.cl,f.pl] = qlevels(f.f,[10, 30, 50, 70, 90, 95, 99, 99.9],f.x{1},f.x{2});
189
190 f.elapsedTime = etime(clock,startTime);
191
192```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003