diagram1

clear
close all
R = 200;
theta = linspace(140,200,501)*pi/180;
z = R*cos(theta) + R;
y = R*sin(theta);
plot(z,y,'LineWidth',2);
axis equal
hold on
plot([-100 R+10],[0 0],'g','LineWidth',2);
k = [2 3];
k = k/norm(k);
z = [0 200*k(2)]-100;
y = [0 200*k(1)]+20;
plot(z,y,'k','LineWidth',2);
hold off

% (y-yr)^2+(z-zr-R)^2 = R^2
% (yo+ky*q)^2 +(-R + zo + kz*q)^2 = R^2
%  yo^2+2ky*q*yo+ky^2*q^2  +  zo^2 -2zo*R + 2zo*kz*q  -2kz*q*R + R^2 + kz^2 = R^2
%  yo^2+zo^2 + (2ky*yo+2kz*zo - 2Kz*R)q + q^2 = 0
%  B = ky*yo+kz*zo-kz*R;
%  C = xo^2+yo^2;

B = y(1)*k(1) + z(1)*k(2) - k(2)*R;
C = y(1)^2 + z(1)^2 - 2*R*z(1);
D = sqrt(B^2-C);
if B<0
    q = -B + D;
else
    q = -B - D;
end
q = -B - sqrt(B^2-C);
yp = y(1) + k(1)*q;
zp = z(1) + k(2)*q;
hold on
plot(zp,yp,'ko');
hold off

hold on
plot([R zp],[0 yp],'m','LineWidth',2);
hold off

kn = [-yp R-zp];
kn = kn/norm(kn);
kperp = [kn(2) -kn(1)];
A = [kn(1) -kperp(1); kn(2) -kperp(2)];
B = [0; -R];
q = A\B;
yi = q(2)*kperp(1);
zi = q(2)*kperp(2);
hold on
plot([0 zi],[0 yi],'k','LineWidth',2);
s = 10;
zs = s*[-kperp(2) -kperp(2)+kn(2) +kn(2)];
ys = s*[-kperp(1) -kperp(1)+kn(1) +kn(1)];
plot(zs+zi,ys+yi,'k');
hold off

zs = [s s 0];
ys = [0 s s];
hold on
plot([zp zp],[yp 0],'k','LineWidth',2);
plot(zp+zs,ys,'k');
hold off

A = [1 -kperp(1); 0 -kperp(2)];
B = -[ zp; zp];
q = A\B;
yo = q(2)*kperp(1);
zo = q(2)*kperp(2);
hold on
plot(zo,yo,'ko');
hold off

G1 = sqrt(zo^2+yo^2);
G2 = yp - yo;
fprintf('G1 %g G2 %g\n',G1,G2);

G1 58.3931 G2 58.3931

hold on
plot([-50 100]*k(2),[-50 100]*k(1),'k--');
plot(-[-100 100]*k(1)+zo,[-100 100]*k(2)+yo,'k--');
hold off

% find Q
q = -(y(1)*k(1) + z(1)*k(2));
yq = y(1)+q*k(1);
zq = z(1)+q*k(2);
Q = sqrt(yq^2+zq^2);
fprintf('Q %g\n',Q);
hold on
%plot([0 zq],[0 yq],'y');
hold off

Q 72.111

% find U, beta, and I - then find G
U = atan(k(1)/k(2));
beta = asin(yp/R);
I = U + beta;
G = Q/(cos(U)+cos(I));
fprintf('G %g\n',G);

G 58.3931