Home > wafo > trgauss > mctp2tc.m

# mctp2tc

## PURPOSE

Frequencies of upcrossing troughs and crests using Markov chain of turning points.

## SYNOPSIS

F=mctp2tc(freqPVR,utc,param,freqPVL)

## DESCRIPTION

 MCTP2TC  Frequencies of upcrossing troughs and crests using Markov chain of turning points.

CALL: f_tc = mctp2tc(f_Mm,utc,param);

where

f_tc  = the matrix with frequences of upcrossing troughs and crests,
f_Mm  = the frequency matrix for the Max2min cycles,
utc   = the reference level,
param = a vector defining the discretization used to compute f_Mm,
note that f_mM has to be computed on the same grid as f_mM.

optional call: f_tc = mctp2tc(f_Mm,utc,param,f_mM);

f_mM  = the frequency matrix for the min2Max cycles.

## CROSS-REFERENCE INFORMATION

This function calls:
 cmat2nt Calculates a counting distribution from a cycle matrix. levels Calculates discrete levels given the parameter matrix. nt2cmat Calculates a cycle matrix from a counting distribution. display Display array. rot90 Rotate matrix 90 degrees.
This function is called by:
 spec2cmat Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) spec2mmtpdf Calculates joint density of Maximum, minimum and period. wafofig10 Intensity of trough-crest cycles computed from St

## SOURCE CODE

001 function F=mctp2tc(freqPVR,utc,param,freqPVL)
002 %MCTP2TC  Frequencies of upcrossing troughs and crests using Markov chain of turning points.
003 %
004 %  CALL: f_tc = mctp2tc(f_Mm,utc,param);
005 %
006 %  where
007 %
008 %        f_tc  = the matrix with frequences of upcrossing troughs and crests,
009 %        f_Mm  = the frequency matrix for the Max2min cycles,
010 %        utc   = the reference level,
011 %        param = a vector defining the discretization used to compute f_Mm,
012 %                note that f_mM has to be computed on the same grid as f_mM.
013 %
014 %  optional call: f_tc = mctp2tc(f_Mm,utc,param,f_mM);
015 %
016 %        f_mM  = the frequency matrix for the min2Max cycles.
017 %
018
019
020 if nargin<4
021   freqPVL=freqPVR;
022 end
023 if nargin<3
024    display('too few inputs parameters, stop')
025    break
026 end
027
028 u=levels(param);
029 udisc=fliplr(u);
030 ntc=sum(udisc>=utc);
031 %ntc1=sum(u<=utc);
032
033 %if udisc(ntc) > utc
034 %   ntc=ntc+1
035 %end
036
037 %if u(ntc1) >= utc
038 %   ntc1=ntc1-1
039 %end
040
041
042 if ntc < 2
043   display('the reference level out of range, stop')
044     break
045     end
046 if ntc > (param(3)-1)
047   display('the reference level out of range, stop')
048     break
049     end
050
051 % normalization of frequency matrices
052
053 nP=length(freqPVR);
054 for i=1:nP,
055     rowsum=sum(freqPVR(i,:));
056     if rowsum~=0
057        freqPVR(i,:)=freqPVR(i,:)/rowsum;
058     end
059 end
060 P=fliplr(freqPVR);
061
062 Ph=rot90(fliplr(freqPVL),-1);
063 for i=1:nP,
064     rowsum=sum(Ph(i,:));
065     if rowsum~=0
066        Ph(i,:)=Ph(i,:)/rowsum;
067     end
068 end
069 Ph=fliplr(Ph);
070
071
072 n=nP; F=zeros(n,n);
073
074 if ntc > n-1
075   display('index for mean-level out of range, stop')
076     break
077     end
078
079
080 %F(1:ntc,1:ntc1)=freqPVL(1:ntc,1:ntc1);
081 F(1:ntc-1,1:(n-ntc))=freqPVL(1:ntc-1,1:(n-ntc));
082
083
084 F=cmat2nt(F);
085
086 for i=2:ntc,
087     for j=ntc:n-1,
088
089     if i<ntc
090
091        Ap=P(i:ntc-1,i+1:ntc); Bp=Ph(i+1:ntc,i:ntc-1);
092        dim_p=ntc-i;
093        tempp=zeros(dim_p,1);
094        I=eye(size(Ap));
095        if i==2
096            e=Ph(i+1:ntc,1);
097          else
098            e=sum(Ph(i+1:ntc,1:i-1)')';
099        end
100        if max(abs(e))>1e-10
101           if dim_p==1
102          tempp(1,1)=(Ap/(1-Bp*Ap)*e);
103           else
104         tempp=Ap*((I-Bp*Ap)\e);
105           end
106         end
107      end
108
109      if j>ntc
110
111        Am=P(ntc:j-1,ntc+1:j); Bm=Ph(ntc+1:j,ntc:j-1);
112        dim_m=j-ntc;
113        tempm=zeros(dim_m,1);
114        Im=eye(size(Am));
115        if j==n-1
116          em=P(ntc:j-1,n);
117        else
118          em=sum(P(ntc:j-1,j+1:n)')';
119        end
120        if max(abs(em))>1e-10
121           if dim_m==1
122         tempm(1,1)=(Bm/(1-Am*Bm)*em);
123           else
124         tempm=Bm*((Im-Am*Bm)\em);
125           end
126        end
127     end
128
129      if (j>ntc)*(i<ntc)
130     F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-ntc:-1:n-j+1)*tempm;
131     F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-j:-1:1)*ones(n-j,1);
132     F(i,n-j+1)=F(i,n-j+1)+ones(1,i-1)*freqPVL(1:i-1,n-ntc:-1:n-j+1)*tempm;
133      end
134      if (j==ntc)*(i<ntc)
135     F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-j:-1:1)*ones(n-j,1);
136     for k=1:ntc
137         F(i,n-k+1)=F(i,n-ntc+1);
138     end
139      end
140      if (j>ntc)*(i==ntc)
141     F(i,n-j+1)=F(i,n-j+1)+ones(1,i-1)*freqPVL(1:i-1,n-ntc:-1:n-j+1)*tempm;
142     for k=ntc:n
143     F(k,n-j+1)=F(ntc,n-j+1);
144     end
145      end
146
147   end
148 end
149 F;
150 F=nt2cmat(F);
151
152
153 %fmax=max(max(F));
154
155 %  contour (u,u,flipud(F),...
156 %fmax*[0.005 0.01 0.02 0.05 0.1 0.2 0.4 0.6 0.8])
157 %  axis([param(1) param(2) param(1) param(2)])
158
159 %  title('Crest-trough density')
160 %  ylabel('crest'), xlabel('trough')
161 %  axis('square')
162 %if mlver>1, commers, end
163
164
165
166
167
168
169
170
171

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