%%
clear, clc, close all

syms t
x(t) = 2*t*sin(2*pi*t);
y(t) =   t*cos(2*pi*t);
z(t) = 3-t^2;
r(t) = 0.01+t^2/6;

%%
syms s
%  Quadrat des Eukldischen Abstands zwischen P(s) und P(t)
d_long = (x(t)-x(s))^2 +(y(t)-y(s))^2 +(z(t)-z(s))^2;
%  Ausmultiplizieren und vereinfachen (vielleicht hilft das)
d(s,t) = simplify(expand(d_long));

% Bedingung d(s,t) >= (r(s)+r(t))^2 
% <==> g(s,t) := d(s,t) - (r(s)+r(t))^2  >=0 
g(s,t) = d(s,t) - (r(s)+r(t))^2;

%% 
disp("Explizites auflösen nach t: g(s,t)=0")
tsol_0 = solve(g(s,t)==0,t)      % nicht explizit nach t

% naechster Versuch
t1=1.5;
disp("Explizites auflösen nach t: g(s="+t1+",t)=0")
tsol_1 = solve(g(t1,t)==0,t)

% und noch ein Versuch
interval_left  = t1*[1-0.1,1     ];
interval_right = t1*[1    ,1+0.1];
%tsol_2 = vpasolve(g(t1,t)==0, t, interval)
tsol_2 = vpasolve(g(t1,t)==0, t, interval_left);
tsol_2 = [tsol_2 vpasolve(g(t1,t)==0, t, interval_right)];

%
interval = [interval_left(1), interval_right(2)];
disp("Loesungen in ["+num2str(interval)+"] :")
disp(tsol_2)
%
fplot(g(t1,t), interval); 
hold on; 
plot(interval,[0,0],'r-', tsol_2,[0,0],'b*')


%% Teste die Funktion
disp("Teste die Funktion")
t2 = bsp_20_fkt(x,y,z,r,t1);

disp ("Solution: "+ t2)

