clc, clear, close

p21=2.416;
delta2=165.2*pi/180;
alpha2=43.3*pi/180;

P21=p21*[cos(delta2); sin(delta2)];


%% L

theta=71.6*pi/180;
beta2=38.4*pi/180;
phi=26.5*pi/180;

w1=[cos(theta); sin(theta)];
w2=[cos(theta+beta2); sin(theta+beta2)];

z1=[cos(phi); sin(phi)];
z2=[cos(phi+alpha2); sin(phi+alpha2)];

a1=w2'*w1-1;
b1=z2'*w1-z1'*w1;
c1=P21'*w1;

a2=w2'*z1-w1'*z1;
b2=z2'*z1-1;
c2=P21'*z1;

A=[a1 b1; a2 b2];
C=[c1; c2];

X_L=inv(A)*C;

W=X_L(1,1);
Z=X_L(2,1);

%% R

sigma=30*pi/180;
gama2=60*pi/180;
sai=150*pi/180;

u1=[cos(sigma); sin(sigma)];
u2=[cos(sigma+gama2); sin(sigma+gama2)];

s1=[cos(sai); sin(sai)];
s2=[cos(sai+alpha2); sin(sai+alpha2)];

a1=u2'*u1-1;
b1=s2'*u1-s1'*u1;
c1=P21'*u1;

a2=u2'*s1-u1'*s1;
b2=s2'*s1-1;
c2=P21'*s1;

A_R=[a1 b1; a2 b2];
C_R=[c1; c2];

X_R=inv(A_R)*C_R;

U=X_R(1,1);
S=X_R(2,1);

%%

W1=W*w1;
W2=W*w2;

Z1=Z*z1;
Z2=Z*z2;

U1=U*u1;
U2=U*u2;

S1=S*s1;
S2=S*s2;

V1=Z1-S1;
V2=Z2-S2;
V=norm(V1);

G=W1+V1-U1;
g=norm(G);

L1=g;
L2=W;
L3=V;
L4=U;

L=[L1, L2, L3, L4]

l=max(L);
s=min(L);
l_s=l+s;
p_q=sum(L)-l_s;

% Grashof=p_q-l_s
if p_q>l_s
    disp('Grashof')
else
    disp('Non-Grashof')
end
line([0, G(1,1)], [0, G(2,1)], 'color', 'k')

line([0, W1(1,1)], [0, W1(2,1)], 'color', 'b')
line([W1(1,1), W1(1,1)+Z1(1,1)], [W1(2,1), W1(2,1)+Z1(2,1)], 'color', 'b')
line([W1(1,1), W1(1,1)+V1(1,1)], [W1(2,1), W1(2,1)+V1(2,1)], 'color', 'b')
line([W1(1,1)+V1(1,1), W1(1,1)+Z1(1,1)], [W1(2,1)+V1(2,1), W1(2,1)+Z1(2,1)], 'color', 'b')
line([G(1,1), W1(1,1)+V1(1,1)], [G(2,1), W1(2,1)+V1(2,1)], 'color', 'b')

line([0, W2(1,1)], [0, W2(2,1)], 'color', 'r')
line([W2(1,1), W2(1,1)+Z2(1,1)], [W2(2,1), W2(2,1)+Z2(2,1)], 'color', 'r')
line([W2(1,1), W2(1,1)+V2(1,1)], [W2(2,1), W2(2,1)+V2(2,1)], 'color', 'r')
line([W2(1,1)+V2(1,1), W2(1,1)+Z2(1,1)], [W2(2,1)+V2(2,1), W2(2,1)+Z2(2,1)], 'color', 'r')
line([G(1,1), W2(1,1)+V2(1,1)], [G(2,1), W2(2,1)+V2(2,1)], 'color', 'r')