001 function f_rfc = mc2rfc(f_xy,paramv,paramu)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020 if nargin<2
021 paramv=[-1, 1, length(f_xy)];
022 paramu=paramv;
023 end
024
025 if nargin<3
026 paramu=paramv;
027 end
028 dd=diag(rot90(f_xy));
029 N=length(f_xy);
030 Splus=sum(f_xy');
031 Sminus=fliplr(sum(f_xy));
032 Max_rfc=zeros(N,1);
033 Min_rfc=zeros(N,1);
034 norm=zeros(N,1);
035 for i=1:N
036 Spm=Sminus(i)+Splus(i)-dd(i);
037 if Spm>0.
038 Max_rfc(i)=(Splus(i)-dd(i))*(Splus(i)-dd(i))/(1-dd(i)/Spm)/Spm;
039 Min_rfc(i)=(Sminus(i)-dd(i))*(Sminus(i)-dd(i))/(1-dd(i)/Spm)/Spm;
040 norm(i)=Spm;
041 end
042 end
043
044
045
046
047
048
049 f_rfc=zeros(N,N);
050 f_rfc(N-1,1)=Max_rfc(N-1);
051 f_rfc(1,N-1)=Min_rfc(2);
052
053 for k=3:N-1
054 for i=2:k-1,
055
056
057
058 AA = f_xy(N-k+1:N-k+i,k-i+1:k);
059 RAA=f_rfc(N-k+1:N-k+i,k-i+1:k);
060 nA=length(AA);
061 MA= Splus(N-k+1:N-k+i);
062 mA=Sminus(N-k+1:N-k+i);
063 normA=norm(N-k+1:N-k+i);
064 MA_rfc=Max_rfc(N-k+1:N-k+i);
065 mA_rfc=Min_rfc(k-i+1:k);
066 SA=sum(sum(AA));
067 SRA=sum(sum(RAA));
068 SMA_rfc=sum(MA_rfc);
069 SMA=sum(MA);
070 DRFC=SA-SMA-SRA+SMA_rfc;
071
072 NT=MA_rfc(1)-sum(RAA(1,:));
073
074
075
076
077
078
079 NT=max(NT,0);
080
081 if NT>1e-6*MA_rfc(1)
082
083 NN=MA-sum(AA');
084 e=(fliplr(mA)-sum(AA))';
085 e=flipud(e);
086 AA=AA+flipud(rot90(AA,-1));
087 AA=rot90(AA);
088 AA=AA-0.5*diag(diag(AA));
089
090
091 for j=1:nA,
092 if normA(j)~=0
093 AA(j,:)=AA(j,:)/normA(j);
094 e(j)=e(j)/normA(j);
095 end
096 end
097 fx=0.;
098
099 if max(abs(e))>1e-7 & max(abs(NN))>1e-7*MA_rfc(1)
100
101
102 I=eye(size(AA));
103
104 if nA==1
105 fx=NN/(1-AA)*e;
106 else
107 fx=NN*((I-AA)\e);
108 end
109 end
110
111 f_rfc(N-k+1,k-i+1)=DRFC+fx;
112 end
113 end
114 m0=max(0,Min_rfc(N)-sum(f_rfc(N-k+2:N,1)));
115 M0=max(0,Max_rfc(N-k+1)-sum(f_rfc(N-k+1,2:k)));
116 f_rfc(N-k+1,1)=min(m0,M0);
117
118 end
119
120 for k=2:N
121 M0=max(0,Max_rfc(1)-sum(f_rfc(1,N-k+2:N)));
122 m0=max(0,Min_rfc(k)-sum(f_rfc(2:k,N-k+1)));
123 f_rfc(1,N-k+1)=min(m0,M0);
124 end
125 f_rfc=f_rfc+rot90(diag(dd),-1);
126 clf
127 subplot(1,2,2)
128 pcolor(levels(paramv),levels(paramu),flipud(f_xy+flipud(rot90(f_xy,-1))))
129 axis([paramv(1), paramv(2), paramu(1), paramu(2)])
130 title('MC-kernel f(x,y)')
131 ylabel('y'), xlabel('x')
132 axis('square')
133
134 subplot(1,2,1)
135 pcolor(levels(paramv),levels(paramu),flipud(f_rfc))
136 axis([paramv(1), paramv(2), paramu(1), paramu(2)])
137 title('Rainflow matrix')
138 ylabel('max'), xlabel('rfc-min')
139 axis('square')
140
141
142
143
144
145
146
147
148