d=152.4;
a=50.8;
b=177.8;
c=228.6;
Rpa=152.4;
w2=10;
alpha2=0;
delta3=30;
theta2=30:360:5;
t2=theta2/180*pi;

Ax=a*cos(theta2);
Ay=a*sin(theta2);
By1=(-Q+sqrt(Q^2-4*P*R))/2*P;
By2=(-Q-sqrt(Q^2-4*P*R))/2*P;
Bx1=S-(2*Ay*By1)/2(Ax-d);
Bx2=S-(2*Ay*By2)/2(Ax-d);
j=sqrt(-1)

S=(a^2-b^2+c^2-d^2)/2*(Ax-d);
P=(Ay*Ay)/((Ax-d)*(Ax-d))+1;
Q=2*Ay(d-S)/(Ax-d);
R=(d-S)*(d-S)-c*c;

theta31=atan((By1-Ay)/(Bx1-Ax));
theta32=atan((By2-Ay)/(Bx2-Ax));
theta41=atan(By1/(Bx1-d));
theta42=atan(By2/(Bx2-d));

w31=a*w2*sin(theta41-theta2)/b/sin(theta31-theta41);
w32=a*w2*sin(theta42-theta2)/b/sin(theta32-theta42);
w41=a*w2*sin(theta2-theta31)/c/sin(theta41-theta32);
w42=a*w2*sin(theta2-theta32)/c/sin(theta42-theta32);

theta31=atan((By1-Ay)/(Bx1-Ax));
theta32=atan((By2-Ay)/(Bx2-Ax));
theta41=atan(By1/(Bx1-d));
theta42=atan(By2/(Bx2-d));

aa1=c*sin(theta41);
aa2=c*sin(theta42);
bb1=b*sin(theta31);
bb2=b*sin(theta32);
cc1=a*alpha2*sin(theta2)+w2*w2*cos(theta2)+b*w31*w31*cos(theta31)-c*w41*w41*cos(theta41);
cc2=a*alpha2*sin(theta2)+w2*w2*cos(theta2)+b*w32*w32*cos(theta32)-c*w42*w42*cos(theta42);
dd1=c*cos(theta41);
dd2=c*cos(theta42);
ee1=b*cos(theta31);
ee2=b*cos(theta32);
ff1=a*alpha2*cos(theta2)-a*w2*w2*sin(theta2)-b*w31*w31*sin(theta31)+c*w41*w41*sin(theta41);
ff2=a*alpha2*cos(theta2)-a*w2*w2*sin(theta2)-b*w32*w32*sin(theta32)+c*w42*w42*sin(theta42);

alpha31=(cc1*dd1-aa1*ff1)/(aa1*ee1-bb1*dd1);
alpha41=(cc1*ee1-bb1*ff1)/(aa1*ee1-bb1*dd1);
alpha32=(cc2*dd2-aa2*ff2)/(aa2*ee2-bb2*dd2);
alpha42=(cc2*ee2-bb2*ff2)/(aa2*ee2-bb2*dd2);
Rpaa1=Rpa*alpha31(-sin(theta31-delta3)+j*cos(theta31+delta3))-Rpa*w31(cos(theta31+delta3)+j*sin(theta31+delta3));
Rpaa2=Rpa*alpha32(-sin(theta32-delta3)+j*cos(theta32+delta3))-Rpa*w32(cos(theta32+delta3)+j*sin(theta32+delta3));

여기 까지 했는데 매트릭스를 사용해서 계산해야되는데 매트릭스를 못쓰겠어... ㅠ