disp '                    ****** 轉軸的可靠性設計 *******'
M=input('            輸入轉軸危險截面上的彎矩(Nmm)     M = ');
T=input('            輸入轉軸危險截面上的扭矩(Nmm)     T = ');
Kmsa=32*M/pi;
fprintf (1,'        對稱循環彎曲應力幅系數             Kmsa = %3.3f \n',Kmsa) 
Kcsa=0.08*Kmsa;
fprintf (1,'        彎曲應力幅標準離差系數             Kcsa = %3.3f \n',Kcsa) 
Kmsm=16*sqrt(3)*T/pi;
fprintf (1,'        穩定扭轉平均應力系數               Kmsm = %3.3f \n',Kmsm) 
Kcsm=0.08*Kmsm;
fprintf (1,'        扭轉平均應力標準離差系數           Kcsm = %3.3f \n',Kcsm) 
rb=Kmsa/Kmsm;
fprintf (1,'        應力幅與平均應力的比值               rb = %3.3f \n',rb) 
Kmrb=sqrt(1+1/rb^2);
fprintf (1,'        應力比均值系數                     Kmrb = %3.3f \n',Kmrb) 
Kmsf=Kmsa*sqrt(1+1/rb^2);
fprintf (1,'        復合疲勞平均應力系數               Kmsf = %3.3f \n',Kmsf) 
Kcsf=Kcsa*sqrt(1+1/rb^2);
fprintf (1,'        復合疲勞平均應力標準離差系數       Kcsf = %3.3f \n',Kcsf) 
Cb=input('            輸入轉軸材料的彎曲強度極限(MPa)  Cb = ');
Csjdc=0.43*Cb;   % 袖珍機械設計師手冊(第2版),P17,表1-18,結構鋼
fprintf (1,'        試件的對稱循環彎曲疲勞極限(MPa)   Csjbc = %3.3f \n',Csjdc) 
B=input('            輸入轉軸的表面質量系數            B = ');
Ec=input('            輸入轉軸的彎曲絕對尺寸系數       Ec = ');
Et=input('            輸入轉軸的扭轉絕對尺寸系數       Et = ');
E=(Ec+Et)/2;
fprintf (1,'        轉軸的彎曲絕對尺寸系數                E = %3.3f \n',E) 
Kc=input('            輸入轉軸的彎曲疲勞應力集中系數   Kc = ');
Kt=input('            輸入轉軸的扭轉疲勞應力集中系數   Kt = ');
Q=input('            輸入轉軸的敏感系數                Q = ');
Kf=1+Q*(Kc*Kt-1);
fprintf (1,'        轉軸的復合疲勞應力集中系數           Kf = %3.3f \n',Kf) 
Cdc=Csjdc*B*E/Kf;
fprintf (1,'        轉軸的對稱循環彎曲疲勞極限(MPa)     Cdc = %3.3f \n',Cdc) 
% 復合疲勞應力下強度的均值Sj,按照應力線與最佳擬合均值線的交點求出
Sj=sqrt(Cdc^2*Cb^2*(1+rb^2)/(Cb^2*rb^2+Cdc^2));
fprintf (1,'        轉軸強度的均值(MPa)                  Sj = %3.3f \n',Sj) 
Cs=0.08*Sj;
fprintf (1,'        轉軸強度的標準離差(MPa)              Cs = %3.6f \n',Cs) 
R=input('            輸入可靠度                        R = ');
% 根據失效概率F求聯結系數z時,用累積分布反函數  z=norminv(F,mu,sigma)
% 根據聯結系數z求失效概率F時,用累積分布函數    F=normcdf(z,mu,sigma)
% 根據聯結系數z求失效頻數f時,用概率密度函數    f=normpdf(z,mu,sigma)
% 正態分布N(mu,sigma),x=z時,mu=0,sigma=1,為標準正態分布,mu和sigma可以省略
z=norminv(1-R);  
F=normcdf(z);    
fprintf (1,'        與可靠度R對應的失效概率               F = %3.6f \n',F)
fprintf (1,'        聯結系數(可靠性系數,安全指數)         z = %3.6f \n',z)
f=normpdf(z);
fprintf (1,'        與聯結系數z對應的失效頻數             f = %3.6f \n',f)
a6=z^2*Cs^2-Sj^2;
fprintf (1,'        聯結方程多項式中6次方項的系數        a6 = %3.4f \n',a6) 
a3=2*Sj*Kmsf;
fprintf (1,'        聯結方程多項式中3次方項的系數        a3 = %3.4f \n',a3) 
a0=z^2*Kcsf^2-Kmsf^2;
fprintf (1,'        聯結方程多項式中的常數項             a0 = %3.4f \n',a0) 
p=[a6 0 0 a3 0 0 a0]; 
% 求解多項式的根
d=roots(p);
disp '        轉軸危險截面直徑(聯結方程多項式的根)'
fprintf (1,'                          d1 = %3.3f mm \n',d(1)) 
fprintf (1,'                          d2 = %3.3f mm \n',d(2)) 
fprintf (1,'                          d3 = %3.3f mm \n',d(3)) 
fprintf (1,'                          d4 = %3.3f mm \n',d(4)) 
fprintf (1,'                          d5 = %3.3f mm \n',d(5)) 
fprintf (1,'                          d6 = %3.3f mm \n',d(6)) 
disp '          ****** 轉軸可靠性設計的有關參數分析 *******'
RR=[0.5 0.6 0.7 0.8 0.9 0.99 0.999 0.9999];
zz=norminv(1-RR);  
figure(1);
plot(zz,RR)
title('\bf   聯結系數z與可靠度R的關系曲線 \rm R=\int e^{-z^{2}/2}dz ')
xlabel('聯結系數 z')
ylabel('可靠度 R')
grid;
aa6=zz.^2*Cs^2-Sj^2;
aa3=2*Sj*Kmsf;
aa0=zz.^2*Kcsf^2-Kmsf^2;
p1=[aa6(1) 0 0 aa3 0 0 aa0(1)]; 
dr1=roots(p1);
fprintf (1,'    可靠度 R=0.5000 時的聯結系數      z1 = %3.6f \n',zz(1))
fprintf (1,'                        轉軸直徑     dr1 = %3.3f mm \n',dr1(6))
p2=[aa6(2) 0 0 aa3 0 0 aa0(2)]; 
dr2=roots(p2);
fprintf (1,'    可靠度 R=0.6000 時的聯結系數      z2 = %3.6f \n',zz(2))
fprintf (1,'                        轉軸直徑     dr2 = %3.3f mm \n',dr2(6))
p3=[aa6(3) 0 0 aa3 0 0 aa0(3)]; 
dr3=roots(p3);
fprintf (1,'    可靠度 R=0.7000 時的聯結系數      z3 = %3.6f \n',zz(3))
fprintf (1,'                        轉軸直徑     dr3 = %3.3f mm \n',dr3(6))
p4=[aa6(4) 0 0 aa3 0 0 aa0(4)]; 
dr4=roots(p4);
fprintf (1,'    可靠度 R=0.8000 時的聯結系數      z4 = %3.6f \n',zz(4))
fprintf (1,'                        轉軸直徑     dr4 = %3.3f mm \n',dr4(6))
p5=[aa6(5) 0 0 aa3 0 0 aa0(5)]; 
dr5=roots(p5);
fprintf (1,'    可靠度 R=0.9000 時的聯結系數      z5 = %3.6f \n',zz(5))
fprintf (1,'                        轉軸直徑     dr5 = %3.3f mm \n',dr5(6))
p6=[aa6(6) 0 0 aa3 0 0 aa0(6)]; 
dr6=roots(p6);
fprintf (1,'    可靠度 R=0.9900 時的聯結系數      z6 = %3.6f \n',zz(6))
fprintf (1,'                        轉軸直徑     dr6 = %3.3f mm \n',dr6(6))
p7=[aa6(7) 0 0 aa3 0 0 aa0(7)]; 
dr7=roots(p7);
fprintf (1,'    可靠度 R=0.9990 時的聯結系數      z7 = %3.6f \n',zz(7))
fprintf (1,'                        轉軸直徑     dr7 = %3.3f mm \n',dr7(6))
p8=[aa6(8) 0 0 aa3 0 0 aa0(8)]; 
dr8=roots(p8);
fprintf (1,'    可靠度 R=0.9999 時的聯結系數      z8 = %3.6f \n',zz(8))
fprintf (1,'                        轉軸直徑     dr8 = %3.3f mm \n',dr8(6))
dz=[dr1(6) dr2(6) dr3(6) dr4(6) dr5(6) dr6(6) dr7(6) dr8(6)];
figure(2);
plot(dz,RR)
title('\bf   軸的直徑d與可靠度R的關系曲線 \it')
xlabel('軸的直徑 d/mm')
ylabel('可靠度 R')
grid;