%Script prova generazione pointcloud con algoritmo quadtree unbalanced
%Profilo NACA-2412 su dominio circolare
%Scritto da Pavan Andrea - 12/02/2022
clear;
clc;

%nodi profilo
airfoilNpoints = 100;
airfoil = naca4_generator('2412',airfoilNpoints);

%nodi contorno
% 'O' mesh
boundaryRadius = 3;     %raggio bordo esterno
boundaryNpoints = 50;      %numero punti bordo esterno
boundary(:,1) = 0.25 + boundaryRadius*cos(0:(2*pi/boundaryNpoints):2*pi);       %coordinate punti bordo esterno
boundary(:,2) = boundaryRadius*sin(0:(2*pi/boundaryNpoints):2*pi);
boundary(end,:) = [];

nodes = [airfoil(1:end-1,:); boundary];      %vettore coordinate nodi


%quadtree ricorsivo
nnodes = 2;     %numero punti ammessi per cella
lmax = 5;       %lunghezza massima lato cella
p = [0];        %vettore parent: p(i)=j indica che cella i è figlia della cella j
x = [0.5*(min(nodes(:,1)+max(nodes(:,1))))];        %coordinata x centro cella i
y = [0.5*(min(nodes(:,2)+max(nodes(:,2))))];        %coordinata y centro cella i
l = [max(nodes(:,1))-min(nodes(:,1))];       %lunghezza lato cella i
n = [length(nodes(:,1))];       %numero punti dentro cella i
while max(n)>nnodes
    nold = n;
    for j=1:length(nold)
        if nold(j)>nnodes || l(j)>lmax
            %decompongo nodo j
            p = [p; j; j; j; j];
            l = [l; l(j)/2; l(j)/2; l(j)/2; l(j)/2];
            x = [x; x(j)+l(j)/4; x(j)-l(j)/4; x(j)-l(j)/4; x(j)+l(j)/4];
            y = [y; y(j)+l(j)/4; y(j)+l(j)/4; y(j)-l(j)/4; y(j)-l(j)/4];
            n = [n; 
                sum((nodes(:,1)>=x(j)).*(nodes(:,2)>=y(j)).*(nodes(:,1)<=x(j)+l(j)/2).*(nodes(:,2)<=y(j)+l(j)/2));
                sum((nodes(:,1)<x(j)).*(nodes(:,2)>=y(j)).*(nodes(:,1)>x(j)-l(j)/2).*(nodes(:,2)<=y(j)+l(j)/2));
                sum((nodes(:,1)<x(j)).*(nodes(:,2)<y(j)).*(nodes(:,1)>x(j)-l(j)/2).*(nodes(:,2)>y(j)-l(j)/2));
                sum((nodes(:,1)>=x(j)).*(nodes(:,2)<y(j)).*(nodes(:,1)<=x(j)+l(j)/2).*(nodes(:,2)>y(j)-l(j)/2))];
            l(j) = -1;
            n(j) = -1;
        end
    end
end
nold = [];
x = x(n>=0);
y = y(n>=0);
l = l(n>=0);
n = n(n>=0);


%generazione pointcloud
P = [];
for j=1:length(n)
    %verifico quali dei 5 punti della cella j sono interni al dominio

    %punto centrale
    Ptrial = [x(j) y(j)];
    if check_internal_points(Ptrial,airfoil)==false && check_internal_points(Ptrial,boundary)==true
        P = [P; Ptrial];
    end

    %punto NE
    Ptrial = [x(j)+l(j)/2 y(j)+l(j)/2];
    if check_internal_points(Ptrial,airfoil)==false && check_internal_points(Ptrial,boundary)==true
        P = [P; Ptrial];
    end

    %punto NW
    Ptrial = [x(j)-l(j)/2 y(j)+l(j)/2];
    if check_internal_points(Ptrial,airfoil)==false && check_internal_points(Ptrial,boundary)==true
        P = [P; Ptrial];
    end

    %punto SW
    Ptrial = [x(j)-l(j)/2 y(j)-l(j)/2];
    if check_internal_points(Ptrial,airfoil)==false && check_internal_points(Ptrial,boundary)==true
        P = [P; Ptrial];
    end

    %punto SE
    Ptrial = [x(j)+l(j)/2 y(j)-l(j)/2];
    if check_internal_points(Ptrial,airfoil)==false && check_internal_points(Ptrial,boundary)==true
        P = [P; Ptrial];
    end
end


%grafico pointcloud
figure(1);
plot(airfoil(:,1),airfoil(:,2),'r.');
hold on;
plot(boundary(:,1),boundary(:,2),'b.');
for j=1:length(P(:,1))
    plot(P(j,1),P(j,2),'k.');
end
title('Grafico pointcloud');
xlabel('x/c');
ylabel('y/c');
axis equal;
axis square;
hold off;


%grafico quadtree
figure(2);
plot(airfoil(:,1),airfoil(:,2),'r.');
hold on;
plot(boundary(:,1),boundary(:,2),'b.');
for j=1:length(l)
    plot([x(j)-l(j)/2 x(j)+l(j)/2], [y(j)+l(j)/2 y(j)+l(j)/2],'b-');
    plot([x(j)-l(j)/2 x(j)-l(j)/2], [y(j)+l(j)/2 y(j)-l(j)/2],'b-');
    plot([x(j)-l(j)/2 x(j)+l(j)/2], [y(j)-l(j)/2 y(j)-l(j)/2],'b-');
    plot([x(j)+l(j)/2 x(j)+l(j)/2], [y(j)+l(j)/2 y(j)-l(j)/2],'b-');
end
title('Grafico suddivisione quadtree');
xlabel('x/c');
ylabel('y/c');
axis equal;
axis square;
hold off;