# Solution 1 with the geometry package
restart;

# Calculation

with(geometry):
local I:
point(A,0,0): point(B,1,0): point(C,x0,y0):
assume(y0>0,-y0*(-1+x0-((1-x0)^2+y0^2)^(1/2))+y0*((x0^2+y0^2)^(1/2)+x0) <> 0):
triangle(t,[A,B,C]):
incircle(ic,t, 'centername'=I):
Cn:=coordinates(I):
centroid(G,t):
CG:=coordinates(G):
a:=-expand(tan(2*arctan(1/3))):
solve({Cn[2]=CG[2],y0/(x0-1)=a}, explicit);
point(C,eval([x0,y0],%)):
answer=FindAngle(line(AB,[A,B]),line(AC,[A,C]));

# Visualization (G is the point of medians intersection)

triangle(t,[A,B,C]):
incircle(ic,t, 'centername'=I):
centroid(G,t):
segment(s,[I,G]):
median(mB,B,t): median(mC,C,t):
draw([A(symbol=solidcircle,color=black),B(symbol=solidcircle,color=black),C(symbol=solidcircle,color=black),I(symbol=solidcircle,color=green),G(symbol=solidcircle,color=blue),t(color=black),s(color=red,thickness=2),ic(color=green),mB(color=blue,thickness=0),mC(color=blue,thickness=0)], axes=none, size=[800,500], printtext=true,font=[times,20]);

Warning, The imaginary unit, I, has been renamed _I

Warning, solve may be ignoring assumptions on the input variables.

# Solution 2 by a direct calculation

# Calculation

restart;
local I;
sinB:=y0/sqrt(x0^2+y0^2):
cosB:=x0/sqrt(x0^2+y0^2):
Sol1:=eval([x,y],solve({y=-(x-1)/3,y=(sinB/(1+cosB))*x}, {x,y})):
tanB:=expand(tan(2*arctan(1/3))):
Sol2:=solve({y0/3=Sol1[2],y0=-tanB*(x0-1)},explicit);
A:=[0,0]: B:=[1,0]: C:=eval([x0,y0],Sol2[2]):
AB:=<(B-A)[]>: AC:=<(C-A)[]>:
answer=arccos(AB.AC/sqrt(AB.AB)/sqrt(AC.AC));

# Visualization

with(plottools): with(plots):
ABC:=curve([A,B,C,A]):
I:=simplify(eval(Sol1,Sol2[2]));
c:=circle(I,eval(Sol1[2],Sol2[2]),color=green):
G:=(A+B+C)/~3;
IG:=line(I,G,color=red,thickness=2):
P:=pointplot([A,B,C,I,G], color=[black$3,green,blue], symbol=solidcircle):
T:=textplot([[A[],"A"],[B[],"B"],[C[],"C"],[I[],"I"],[G[],"G"]], font=[times,20], align=[left,below]):
M:=plot([[(C+t*~((A+B)/2-C))[],t=0..1],[(B+t*~((A+C)/2-B))[],t=0..1]], color=blue, thickness=0):
display(ABC,c,IG,P,T,M, scaling=constrained, axes=none,size=[800,500]);

Warning, The imaginary unit, I, has been renamed _I