function [P, T] = presekKrivulj(p, pdot, intp, q, qdot, intq, h)
%[P, T] = presekKrivulj(p, pdot, [a, b], q, qdot, [c, d], h) poisce 
%presek dveh parametriziranih krivulj, p: [a, b] -> R^2 in 
%q: [c, d] -> R^2. (Intervala [a, b] in [c, d] razdelimo s korakom h, 
%da dobimo delilne tocke za konstrukcijo lomljenk.)
%Vrne nabor krajevnih vektorjev presecisc P = [[x1; y1], ..., [xk; yk]]
%in nabor pripadajocih parametrov T = [[t1; u1], ..., [tk; uk]].

%pripravimo nabor parametrov t in u iz [a, b] in [c, d] s korakom h...
t = intp(1):h:intp(2);
u = intq(1):h:intq(2);
%... in poracunamo tocke na krivuljah pri tem naboru parametov.
A = p(t);
B = q(u);
%Nato poiscemo presecisca dobljenih lomljenk (skupaj s priblizki pripadajocih parametrov).
[P, T] = presecisca(A, B, h, intp(1), intq(1));
%Dobljene priblizke izboljsamo z Newtonovo metodo.
for k = 1:length(T)
	%Iscemo nicle F(t, u) = p(t) - q(u), ...
	F = @(T) (p(T(1)) - q(T(2)));
	%... ki ima Jacobijevo matriko JF(t, u) = [pdot(t), -qdot(u)].
	JF = @(T) [pdot(T(1)), -qdot(T(2))];
	%Pozenemo Newtonovo metodo, dobimo izboljsane vrednosti parametrov, ...
	T(:, k) = newton(F, JF, T(:, k), 1e-10, 20);
	%... ki jih nato vstavimo v originalno parametrizacijo, da dobimo presecisca.
	P(:, k) = p(T(1, k));
end