Home > wafo > wsim > private > rfc2load_fat.m

## PURPOSE

Recontructs a load process given the frequency matrix (and residual).

## DESCRIPTION

``` RFC2LOAD_FAT  Recontructs a load process given the frequency matrix (and residual).

where

f        = the frequency matrix for the rainflow count,
residual = the residual (optional input argument),
if left out the program will generate a
stationary residual using the frequency
matrix.
num_cc   = the expected number of cycles,
if num_cc=[] the program will continue
until the residual is empty,```

## CROSS-REFERENCE INFORMATION

This function calls:
 comb2pro Probability vector of placing a cycle on a specific place given all combinations dtp2rfm Calculates rainflow matrix from discrete turning points. fr2comb Combination matrix for placing out cycles given the frequency matrix fr2res Generates a stationary residual from the frequency matrix. tp2rfc Finds the rainflow cycles from the sequence of turning points. error Display message and abort function. triu Extract upper triangular part. warning Display warning message; disable or enable warning messages.
This function is called by:
 rfm2dtp Reconstructs a sequence of turning points from a rainflow matrix.

## SOURCE CODE

```001 function [X,res,comb,f]=rfc2load_fat(f,res,num_cc)
002 %RFC2LOAD_FAT  Recontructs a load process given the frequency matrix (and residual).
003 %
004 % CALL:  X = rfc2load_fat(f,num_cc,residual);
005 %
006 %  where
007 %
008 %        X        = the reconstructed load,
009 %        f        = the frequency matrix for the rainflow count,
010 %        residual = the residual (optional input argument),
011 %                   if left out the program will generate a
012 %                   stationary residual using the frequency
013 %                   matrix.
014 %        num_cc   = the expected number of cycles,
015 %                   if num_cc=[] the program will continue
016 %                   until the residual is empty,
017
018 %  Copyright 1993, Mats Frendahl & Igor Rychlik,
019 %  Dept. of Math. Stat., University of Lund.
020
021 % Copyright (c) 2004 by Pär Johannesson
022
023 % Tested  on Matlab  6.5
024 %
025 % History:
026 % Adapted to WAFO by PJ (Pär Johannesson) 16-Feb-2004
027 %   The function 'rfc2load' originally from FAT (Fatigue Analysis Toolbox)
028 %   FAT is a predecessor of WAFO
029 % Changed by PJ 19-Feb-2004
030
031 %%%%
032 % Check input arguments
033
034 ni = nargin;
035 no = nargout;
036 error(nargchk(1,3,ni));
037
038 if ni<2, res=[]; end
039 if ni<3, num_cc=[]; end
040
041 % Check if freq. matrix is interger matrix.
042 if ~(round(f)==f)
043     warning('The frequency matrix is not interger matrix.   The matrix will be scaled.')
044     f=floor(100*f);
045 end
046
047 N=length(f);
048
049 % If empty residual, then calculate a stat. residual
050 % using the frequency matrix.
051 if isempty(res)
052     res0=fr2res(f);
053     % Don't give a proper residual, modify it!
054     % Correction by PJ
055     [RFC,RFC1,res] = tp2rfc(res0(:));  % Ger proper residual
056     resWAFO = N-res+1;              % Convert to WAFO-def
057     RFMres = dtp2rfm(resWAFO,N,'CS'); % Rainflow matrix of residual
058     RFMres = fliplr(RFMres)';       % Convert to FAT-def
059     f = f-RFMres;                   % Remove cycles in res from rainflow matrix
060 end
061
062 % Find the largest amplitude. Look up how many there are,
063 % generate this number of large amplitudes, reset the
064 % frequency matrix to 0 for this kind of amplitude.
065 largest_res=sort([min(min(res)) max(max(res))]);
066 index_i=largest_res(1); index_j=N-largest_res(2)+1;
067 n_large_res=f(index_i,index_j); f(index_i,index_j)=0;
068 xtra_cc=[]; for i=1:n_large_res, xtra_cc=[xtra_cc largest_res]; end
069
070 r=res(:)';
071
072 % Find the larges amplitude in the residual and put the
073 % extra large amplitudes here.
074 undone=1;
075 for i=1:length(r)-1
076     if undone
077         if ( sum(largest_res==[r(i) r(i+1)])==2 )
078             r=[r(1:i-1) xtra_cc r(i:length(r))];
079             undone=0;
080         elseif ( sum(fliplr(largest_res)==[r(i) r(i+1)])==2 )
081             r=[r(1:i-1) fliplr(xtra_cc) r(i:length(r))];
082             undone=0;
083         end
084     end
085 end
086
087
088 % If the number of requested cycles are not empty then find
089 % the number of extected runs to give this number of cycles.
090 if num_cc~=[]
091     % Norm frequency matrix to probability matrix.
092     f=f/sum(sum(f));
093
094     % Calculate the number of extected runs to get the requested number
095     % of cycles.
096     num_cc_limit=2*num_cc+4;
097
098     % Scale the frequency matrix to whole numbers and make it large to
099     % be stationary.
100     f=floor(max([1000*N 1e6])*f);
101 end
102
103 % Calculated the combination matrix
104 comb=fr2comb(f,r);
105
106 % Initiate the residual.
107 X=r(1);
108
109 again=1; num_runs=0;
110
111 while again
112
113     num_runs=num_runs+1;
114
115     h=0;
116
117     if r(1)>r(2)
118         % disp('upslope')
119
120         if r(1)>(r(2)+1)
121             m=r(1); M=r(2); j=N-m+1; i=m-1:-1:M+1;
122             combinations=comb(i,j);
123             freq=f(i,j);
124             prob=comb2pro(combinations,freq);
125             h=min(find(rand<cumsum(prob)))-1;
126         end
127
128         if h>0
129             C=m-h;
130             i=C:m;
131             if f(C,j)>0
132                 f(C,j)=f(C,j)-1;
133                 i=(C+1):r(1);
134                 j=N+1-r(1);
135                 comb(i,j)=comb(i,j)-1;
136                 r=[m-1 C r];
137             else
138                 'Warning 1', break
139                 r=[r(1)-1 r(2:length(r))];
140             end
141         else
142             j=N+1-r(1);
143             if r(1)>(r(2)+1)
144                 i=(r(2)+1):(r(1)-1);
145                 comb(i,j)=comb(i,j)-1;
146             end
147             r(1)=r(1)-1;
148
149         end
150
151     else
152
153         % disp('downslope')
154         if r(1)<r(2)-1
155             m=r(2); M=r(1); i=M; j=M+1:m-1; j=N-j+1;
156             combinations=comb(i,j)';
157             freq=f(i,j)';
158             prob=comb2pro(combinations,freq);
159             h=min(find(rand<cumsum(prob)))-1;
160         end
161
162         if h>0
163             j=M:M+h; j=N-j+1;
164             c0=M+h;
165             c=N+1-c0;
166             if f(i,c)>0
167                 f(i,c)=f(i,c)-1;
168                 i=r(1);
169                 j=N+1-((i+h-1:-1:r(1)));
170                 comb(i,j)=comb(i,j)-1;
171                 r=[M+1 c0 r];
172             else
173                 'Warning 2', break
174                 r=[r(1)+1 r(2:length(r))];
175             end
176         else
177             i=r(1);
178             j=N+1-(r(1):(r(2)-1));
179             comb(i,j)=comb(i,j)-1;
180             r(1)=r(1)+1;
181
182         end
183     end
184
185     if r(1)==r(2), r=r(2:length(r)); end
186
187     if (length(r)>=2)
188         again=1;
189     else
190         disp('   The residual is empty. Program will terminate.')
191         again=0;
192     end
193
194     X=[X r(1)];
195
196     n_tmp=length(X);
197
198     % This is an equivalent of shave, but faster
199
200     if n_tmp>2
201         if (X(n_tmp-2)<X(n_tmp-1))&(X(n_tmp-1)<X(n_tmp))
202             X=[X(1:n_tmp-2) X(n_tmp)];
203         elseif (X(n_tmp-2)>X(n_tmp-1))&(X(n_tmp-1)>X(n_tmp))
204             X=[X(1:n_tmp-2) X(n_tmp)];
205         end
206     end
207
208     comb=fliplr(triu(fliplr(comb),1));
209
210     % If the number of requested cycles are not empty then check
211     % if the number of step (length of load) exceeds the expected
212     % number.
213     if num_cc~=[]
214         if length(X)>num_cc_limit
215             disp('   The number of requested cycles has been reached.')
216             disp('   Program will terminate.')
217             again=0;
218         end
219     end
220
221 end
222
223 X = X';
224```

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