function [ output_args ] = lab3( input_args )
    clc, clf, clear
    format long g

    hold on
    axis equal

	plot([-100, 100], [0, 0])
	plot([0, 0], [-100, 100])

	% ORIGINAL DOTS START
		figure(1) % figure showing the cross cut
		testpunkterx =  [0 8 12 6.5 2];%[0 1 3 4 2]
		testpunktery =  [0 -3 2 9.5 5.5];%[0 7 8 5 2]
		% plotting original dots
		plot(testpunkterx, testpunktery, 'x')

		% guessing the b and c values
		b = [8 27.5]; c = [22.5 -14];
	% ORIGINAL DOTS END

	% TEST SHIT
% 		testpunkterx =  [0 1 3 4 2]
% 		testpunktery =  [0 7 8 5 2]
%
% % 		b = [3 13]; c = [12 8];
%
% 		b = [13 7]; c = [12 -4];
% 		tempx = cos(-pi/3) * testpunkterx - sin(-pi/3) * testpunktery;
% 		tempy = sin(-pi/3) * testpunkterx + cos(-pi/3) * testpunktery;
% 		testpunkterx = tempx;
% 		testpunktery = tempy;

%
% 		testpunkterx = [0 1 1 0 0];
% 		testpunktery = [0 0 0 1 1];
% 		b = [1 0]; c = [0 1];
% 		tempx = cos(-pi/3) * testpunkterx - sin(-pi/3) * testpunktery;
% 		tempy = sin(-pi/3) * testpunkterx + cos(-pi/3) * testpunktery;
% 		testpunkterx = tempx;
% 		testpunktery = tempy;
	% TEST SHIT


	% plotting original dots
	plot(testpunkterx, testpunktery, 'x')
    t = 0 : 0.01 : 1;
	p1 = [0 0]; p2 = [0 0];
    crosscut = cubicbezier(p1, p2, b, c, t');
    % plotting my best guess
    plot(crosscut(:, 1), crosscut(:, 2), 'red')

    % guessing the t values
    t = [0.1 0.3 0.5 0.9];
    dxnorm = 1;
    iter = 0;
    x = [b c t]';

	% symbolisk lösning ser ut såhär:
%     syms bx by cx cy t1 t2 t3 t4
%     Fsymb = [3*t1*bx*(1 - t1)^2 + 3*cx*(1 - t1)*t1^2 - testpunkterx(2)
%              3*t1*by*(1 - t1)^2 + 3*cy*(1 - t1)*t1^2 - testpunktery(2)
%              3*t2*bx*(1 - t2)^2 + 3*cx*(1 - t2)*t2^2 - testpunkterx(3)
%              3*t2*by*(1 - t2)^2 + 3*cy*(1 - t2)*t2^2 - testpunktery(3)
%              3*t3*bx*(1 - t3)^2 + 3*cx*(1 - t3)*t3^2 - testpunkterx(4)
%              3*t3*by*(1 - t3)^2 + 3*cy*(1 - t3)*t3^2 - testpunktery(4)
%              3*t4*bx*(1 - t4)^2 + 3*cx*(1 - t4)*t4^2 - testpunkterx(5)
%              3*t4*by*(1 - t4)^2 + 3*cy*(1 - t4)*t4^2 - testpunktery(5)];
%     symbolicJ = jacobian(Fsymb, [bx by cx cy t1 t2 t3 t4])
%	och detta nedan ska göras i varje iteration för att uppdatera f och J
%		f = subs(Fsymb, [bx by cx cy t1 t2 t3 t4], x)
% 		J = subs(symbolicJ, [bx by cx cy t1 t2 t3 t4], x)

    while dxnorm > 1e-14 && iter < 15
        f = [	3*x(5)*x(1)*(1 - x(5)).^2 + 3*x(3)*(1 - x(5))*x(5)^2 - testpunkterx(2)
				3*x(5)*x(2)*(1 - x(5)).^2 + 3*x(4)*(1 - x(5))*x(5)^2 - testpunktery(2)
				3*x(6)*x(1)*(1 - x(6)).^2 + 3*x(3)*(1 - x(6))*x(6)^2 - testpunkterx(3)
				3*x(6)*x(2)*(1 - x(6)).^2 + 3*x(4)*(1 - x(6))*x(6)^2 - testpunktery(3)
				3*x(7)*x(1)*(1 - x(7)).^2 + 3*x(3)*(1 - x(7))*x(7)^2 - testpunkterx(4)
				3*x(7)*x(2)*(1 - x(7)).^2 + 3*x(4)*(1 - x(7))*x(7)^2 - testpunktery(4)
				3*x(8)*x(1)*(1 - x(8)).^2 + 3*x(3)*(1 - x(8))*x(8)^2 - testpunkterx(5)
				3*x(8)*x(2)*(1 - x(8)).^2 + 3*x(4)*(1 - x(8))*x(8)^2 - testpunktery(5)	];

        J = [	3*x(5)*(1 - x(5))^2, 0, 3*(1 - x(5))*x(5)^2, 0, 3*(x(1)*((1 - x(5))^2 - 2*x(5)*(1-x(5))) + x(3)*(2*x(5) - 3*x(5)^2)), 0, 0, 0
				0, 3*x(5)*(1 - x(5))^2, 0, 3*(1 - x(5))*x(5)^2, 3*(x(2)*((1 - x(5))^2 - 2*x(5)*(1-x(5))) + x(4)*(2*x(5) - 3*x(5)^2)), 0, 0, 0
				3*x(6)*(1 - x(6))^2, 0, 3*(1 - x(6))*x(6)^2, 0, 0, 3*(x(1)*((1 - x(6))^2 - 2*x(6)*(1-x(6))) + x(3)*(2*x(6) - 3*x(6)^2)), 0, 0
				0, 3*x(6)*(1 - x(6))^2, 0, 3*(1 - x(6))*x(6)^2, 0, 3*(x(2)*((1 - x(6))^2 - 2*x(6)*(1-x(6))) + x(4)*(2*x(6) - 3*x(6)^2)), 0, 0
				3*x(7)*(1 - x(7))^2, 0, 3*(1 - x(7))*x(7)^2, 0, 0, 0, 3*(x(1)*((1 - x(7))^2 - 2*x(7)*(1-x(7))) + x(3)*(2*x(7) - 3*x(7)^2)), 0
				0, 3*x(7)*(1 - x(7))^2, 0, 3*(1 - x(7))*x(7)^2, 0, 0, 3*(x(2)*((1 - x(7))^2 - 2*x(7)*(1-x(7))) + x(4)*(2*x(7) - 3*x(7)^2)), 0
				3*x(8)*(1 - x(8))^2, 0, 3*(1 - x(8))*x(8)^2, 0, 0, 0, 0, 3*(x(1)*((1 - x(8))^2 - 2*x(8)*(1-x(8))) + x(3)*(2*x(8) - 3*x(8)^2))
				0, 3*x(8)*(1 - x(8))^2, 0, 3*(1 - x(8))*x(8)^2, 0, 0, 0, 3*(x(2)*((1 - x(8))^2 - 2*x(8)*(1-x(8))) + x(4)*(2*x(8) - 3*x(8)^2))	];

        dx = -J\f

        plot(x(1), x(2), '*')
        plot(x(3), x(4), '.')

        % uppdaterar x med nya fina värden
        x = x + dx;

        % den här ska bli liten och fin
        dxnorm = norm(dx, inf)
        iter = iter + 1;
	end, x

	% plotting the result of the maximization
    t = 0 : 0.01 : 1;
    crosscut = cubicbezier(p1, p2, [x(1) x(2)], [x(3) x(4)], t');
    plot(crosscut(:, 1)+12, crosscut(:, 2), 'green')

	figure(2) % rotation of the bezier around z axis


	X = [], Y = [], Z = []
	x = crosscut(:, 1)+12;
	y = crosscut(:, 2);
	for angle = 0 : pi/64 : 2*pi;
		X = [X, x * cos(angle)];
		Y = [Y, x * sin(angle)];
		Z = [Z crosscut(:, 2)];
	end

	colormap('copper')
	surf(X, Y, Z, 'EdgeColor', 'none')
	axis equal

	hold on

	plotGuideCircles()
end

function out = cubicbezier(P1, P2, b, c, t)
    F3 = [(1-t).^3 3*t.*(1-t).^2 3*(1-t).*t.^2 t.^3];

    out = F3 * [P1; b; c; P2];
end

function out = plotGuideCircles()

	circs = [40, -3; 48, 2; 37, 9.5; 28, 5.5];
	for i = 1 : length(circs)
		[I, J, K] = circle(pi/20, circs(i, 1)/2, circs(i, 2));
		h = plot3(I, J, K, 'black')
		set(h, 'linewidth', 1.2)
	end
end

function [X, Y, Z] = circle(accuracy, x, y)
	X = [], Y = [], Z = []
	for angle = 0 : accuracy : 2*pi;
		X = [X, x * cos(angle)];
		Y = [Y, x * sin(angle)];
		Z = [Z y];
	end
end

function out = funcBez(p1, p2, b, c, t)

end

function out = redbez(b, c, t)
    F = [3*t.*(1 - t).^2 3*(1 - t).*t.^2];

    out = F * [b; c];
end