clear all
clc
x=1e7;
chou=1e5;
C=0;
for i=0:100000
R_o = normrnd(199.9588,0.1675,[1 1e7]);
L_o = normrnd(399.9615,0.1699,[1 1e7]);
Ca_o = normrnd(0.043,0.0129,[1 1e7]);
Cb_o = normrnd(0.043,0.0129,[1 1e7]);
Cc_o = normrnd(0.043,0.0129,[1 1e7]);
k_R = unidrnd(x,[1,chou]);
k_L = unidrnd(x,[1,chou]);
k_Ca = unidrnd(x,[1,chou]);
k_Cb = unidrnd(x,[1,chou]);
k_Cc = unidrnd(x,[1,chou]);
R_k = R_o(k_R);
L_k = L_o(k_L);
Ca_k = Ca_o(k_Ca);
Cb_k = Cb_o(k_Cb);
Cc_k = Cc_o(k_Cc);
R_mean = mean(R_k);
L_mean = mean(L_k);
Ca_mean = mean(Ca_k);
Cb_mean = mean(Cb_k);
Cc_mean = mean(Cc_k);
%%
y_mean(i+1)=sqrt((L_mean+R_mean-Ca_mean-Cb_mean-Cc_mean-3*C).^2-10000)
y=sqrt((L_k+R_k-Ca_k-Cb_k-Cc_k-3*C).^2-10000);
k=length(find(y<590));
Pf=1-k/chou
C=C+0.05
if Pf<0.95
    break;
end
end
plot(linspace(0.05,C,length(y_mean)),y_mean)
xlabel('¼ä϶´óС')
ylabel('»¬¿é×î´óÊä³öλÒÆƽ¾ùÖµ/mm')
title('¿¼ÂÇ·ÇÀÛ»ýÐÔÒòËØʱ»¬¿é×î´óÊä³öλÒÆ')