clear,clc,close
alpha_e1=15.11*pi/180;
alpha_e2=alpha_e1;
alpha_i1=11.11*pi/180;
alpha_i2=alpha_e1;
alpha_f1=78.35*pi/180;
alpha_f2=alpha_e1;
dc=25;
R2=dc/2;
dm_1=85;
dm_2=85;
Z1=19;
Z2=19;
Lwe_1=14.59;
Lwe_2=14.59;
Dw1=12.425;
Dw2=13.507;
Dwe_1=0.5*(Dw1+Dw2);
Dwe_2=Dwe_1;
Fa=5000;
Fr=15000;
F0=1000;
Mk=300000;
ff=@(X)PrimaryFun(X,dc,dm_1,Z1,dm_2,Z2,Lwe_1,Lwe_2,alpha_e1,alpha_e2,alpha_i1,alpha_i2,alpha_f1,alpha_f2,F0,Fa,Fr,Mk);
X0=[0.1,0.1,0.1];
x=fsolve(ff,X0);
delta_a=x(1);
delta_r=x(2);
theta=x(3);K_ne1=K_ne(alpha_e1,alpha_i1,alpha_f1,Lwe_1);
delta_0=(F0/(Z1*K_ne1*sin(alpha_e1)^2.11))^0.9;
P_phi=zeros(1,Z1);
phi=zeros(1,Z1);
for i=1:Z1%phi_1i=2*pi*(i-1)/Z1;phi(i)=phi_1i;delta_r_theta_L_i=R2*cos(phi_1i)*theta;detla_a_theta_L_i=0.5*dm_1*cos(phi_1i)*theta;delta_r_L_i=delta_r*cos(phi_1i)+delta_r_theta_L_i;delta_a_L_i=delta_a+delta_0+detla_a_theta_L_i;delta_n_L_i=delta_r_L_i*cos(alpha_e1)+delta_a_L_i*sin(alpha_e1);if delta_n_L_i >=0Q_e_L_i=K_ne1*delta_n_L_i^1.11;P_phi(i)=Q_e_L_i;end
end
P_max=max(P_phi);
plot(phi,P_phi);function y=PrimaryFun(X,dc,dm_1,Z1,dm_2,Z2,Lwe_1,Lwe_2,alpha_e1,alpha_e2,alpha_i1,alpha_i2,alpha_f1,alpha_f2,F0,Fa,Fr,Mk)
%Writen by ZSY
delta_a=X(1);
delta_r=X(2);
theta=X(3);
R2=dc/2;
K_ne1=8.06e4*Lwe_1^0.89*(1+(sin(alpha_e1+alpha_f1)/sin(alpha_i1+alpha_f1))^0.9*cos(alpha_e1-alpha_i1))^(-1.11);
K_ne2=8.06e4*Lwe_2^0.89*(1+(sin(alpha_e2+alpha_f2)/sin(alpha_i2+alpha_f2))^0.9*cos(alpha_e2-alpha_i2))^(-1.11);delta_0=(F0/(Z1*K_ne1*sin(alpha_e1)^2.11))^0.9;
F1_L=0;
F2_L=0;
F3_L=0;
for i=1:Z1phi_1i=2*pi*(i-1)/Z1;delta_r_theta_L_i=R2*cos(phi_1i)*theta;detla_a_theta_L_i=0.5*dm_1*cos(phi_1i)*theta;delta_r_L_i=delta_r*cos(phi_1i)+delta_r_theta_L_i;delta_a_L_i=delta_a+delta_0+detla_a_theta_L_i;delta_n_L_i=delta_r_L_i*cos(alpha_e1)+delta_a_L_i*sin(alpha_e1);if delta_n_L_i >=0Q_e_L_i=K_ne1*delta_n_L_i^1.11;F1_L=F1_L+Q_e_L_i*cos(alpha_e1)*cos(phi_1i);F2_L=F2_L+Q_e_L_i*sin(alpha_e1);F3_L=F3_L+0.5*dm_1*cos(phi_1i)*Q_e_L_i*sin(alpha_e1)+0.5*dc*cos(phi_1i)*Q_e_L_i*cos(alpha_e1);end
end
F1_R=0;
F2_R=0;
F3_R=0;
for i=1:Z2phi_2i=2*pi*(i-1)/Z2;delta_r_theta_R_i=R2*cos(phi_2i)*theta;detla_a_theta_R_i=0.5*dm_2*cos(phi_2i)*theta;delta_r_R_i=delta_r*cos(phi_2i)-delta_r_theta_R_i;delta_a_R_i=-delta_a-delta_0-detla_a_theta_R_i;delta_n_R_i=delta_r_R_i*cos(alpha_e2)+delta_a_R_i*sin(alpha_e2);if delta_n_R_i >=0Q_e_R_i=K_ne2*delta_n_R_i^1.11;F1_R=F1_R+Q_e_R_i*cos(alpha_e2)*cos(phi_2i);F2_R=F2_R+Q_e_R_i*sin(alpha_e2);F3_R=F3_R+0.5*dm_2*cos(phi_2i)*Q_e_R_i*sin(alpha_e2)+0.5*dc*cos(phi_2i)*Q_e_R_i*cos(alpha_e2);end
end
F1=-F1_L-F1_R+Fr;
F2=-F2_L+F2_R+Fa;
F3=-F3_L+F3_R+Mk;
y=[F1;F2;F3];
end
