%% Rotationskörper: Kreisel

clear; clc
syms t

f(t) = piecewise(t <= 2, t.^2 / 2, t > 2, (t - 7).^2 * 2 / 25);
x = linspace(0, 7, 50);
r = f(x);
r = eval(r);  % um von symbolischen Einträgen auf Gleitkommaeinträge umzuschalten

[spinner_x, spinner_y, spinner_z] = cylinder(r);
spinner_z = spinner_z * 5;
surf(spinner_x, spinner_y, spinner_z)

%% Rotationskörper: Maturaaufgabe Herbsttermin 2021/22, BHS-HTL2, Aufgabe 7B
% Quelle: https://www.matura.gv.at/index.php?eID=dumpFile&t=f&f=5601&token=c89c06caa876c1dfbe0eb705d5172cae4c980ee8

clear; clc
syms t

f2(t) = 5/4 * sqrt(-t^2 + 20.8*t - 50.4);
x_zeros = solve(f2(t) == 0, t);

x_max = max(x_zeros);
x_max = eval(x_max);

x_0 = 6.2;

f1(t) = f2(x_0) + (t - x_0)*subs(diff(f2, t), t, x_0);  % Tangentengleichung

f(t) = piecewise( ...
    t <= x_0, f1(t), ...  % wenn t <= x0, dann f1(t),
    t > x_0, f2(t) ...  % und sonst f2(t).
);

x = linspace(0, x_max, 50);
r = eval(f(x));

[balloon_x, balloon_y, balloon_z] = cylinder(r);
balloon_z = balloon_z * x_max;

surf(balloon_x, balloon_y, balloon_z)
axis equal

%% Extraktion von ebenen Kurven aus Surfaces

clear; clc

[cyl_x, cyl_y, cyl_z] = cylinder(5, 35);
curve_x = cyl_x(1, :);
curve_y = cyl_y(1, :);

clear cyl_x cyl_y cyl_z  % volle Matrizen werden nicht mehr gebraucht

plot(curve_x, curve_y)  % Grundfläche!
axis equal

%% Oktaeder, Runde 1: fill3

verts = [
    -1, -1, 0;
     1, -1, 0;
     1,  1, 0;
    -1,  1, 0;
     0,  0, sqrt(2);
     0,  0, -sqrt(2);
];

faces_up = [
    1 2 5;
    2 3 5;
    3 4 5;
    4 1 5;
];
faces_down = [
    1 2 6;
    2 3 6;
    3 4 6;
    4 1 6;
];

hold on
for face = faces_up'  % Iteration über Spalten von Transponierter = Iteration über Zeilen!
    fill3(verts(face, 1), verts(face, 2), verts(face, 3), 'y')
end
for face = faces_down'
    fill3(verts(face, 1), verts(face, 2), verts(face, 3), 'm')
end

axis equal
alpha 0.3
view([-36, 10])
hold off

%% Oktaeder, Runde 2: trisurf

verts = [
    -1, -1, 0;
     1, -1, 0;
     1,  1, 0;
    -1,  1, 0;
     0,  0, sqrt(2);
     0,  0, -sqrt(2);
];

faces_up = [
    1 2 5;
    2 3 5;
    3 4 5;
    4 1 5;
];
faces_down = [
    1 2 6;
    2 3 6;
    3 4 6;
    4 1 6;
];
faces = [faces_up; faces_down];

hold on

trisurf(faces, verts(:,1), verts(:,2), verts(:,3))
idx = [1 2 3 4];
fill3(verts(idx, 1), verts(idx, 2), verts(idx, 3), 'r')

alpha 0.3
axis equal
view([-36, 10])
hold off

%% Würfelexplosion / -zerlegung in Tetraeder

clear; clc

N = 500;
verts = [
    -N, -N, -N;
     N, -N, -N;
    -N,  N, -N;
     N,  N, -N;
    -N, -N,  N;
     N, -N,  N;
    -N,  N,  N;
     N,  N,  N
];
tetraeders = [0 1 2 4; 2 1 3 7; 1 4 5 7; 1 2 4 7; 4 2 6 7];  % C indexing
tetraeders = tetraeders + 1;  % Matlab indexing

num_tetraeders = size(tetraeders, 1);

figure(1)
tetramesh(tetraeders, verts, 'FaceAlpha', 0.3, 'visible', 'on')
axis equal
xlabel("x"); ylabel("y"); zlabel("z");
title(["Würfel in ", num2str(num_tetraeders), " Tetraeder zerschnitten"]);
view([-14,42])

figure(2)
hold on

cube_center = mean(verts);
for k = 1:num_tetraeders
    tetraeder_idx = tetraeders(k, :);
    T = verts(tetraeder_idx, :);  % Koordinatenmatrix von aktuellem Tetraeder
    T_center = mean(T);
    shift_vector = (T_center - cube_center) * 0.3;  % 0.3 kontrolliert Expansion
    T = T + repmat(shift_vector, 4, 1);  % repmat: wiederhole Matrix 4x in Spaltenrichtung, 1x in Zeilenrichtung
    
    local_connectivity = [1 2 3; 1 2 4; 2 3 4; 1 3 4];
    trisurf(local_connectivity, T(:, 1), T(:, 2), T(:, 3))
end

alpha 0.7
axis equal
view([-14, 42])
