function tau=solve_tau(Pb,Pu,radius_k)

A=(norm(Pb))^2+(norm(Pu))^2-2*Pb.'*Pu;

B=6*Pb.'*Pu-4*(norm(Pu))^2-2*(norm(Pb))^2;

C=4*(norm(Pu))^2+(norm(Pb))^2-4*Pu.'*Pb-radius_k^2;

temp=roots([A,B,C]);

if (imag(temp(1,1))~=0 || imag(temp(2,1))~=0)
    temp(1,1)=real(temp(1,1))+imag(temp(1,1));
    temp(2,1)=real(temp(2,1))+imag(temp(2,1));
end

if 0<=temp(1,1)<=2
    tau=temp(1,1);
else if 0<=temp(2,1)<=2
        tau=temp(2,1);
    end
end









