% *************************************************************
% *************************************************************
%
% SCRIPT que calcula pelo modelo do Vanio os elementos da 
% geometria diferencial na superficie gerada pela funcao Shoe 
%
% *************************************************************
% *************************************************************

% =============================================================
% Definicao da funcao Shoe
% =============================================================
syms u v;

f = u;

g = v;

h = ((1/3) * (u^3)) - ((1/2) * (v^2));

% =============================================================
% Geracao dos pontos da superficie
% =============================================================
piso_u = -2;
teto_u =  2;
piso_v = -2;
teto_v =  2;
delta_u = 0.2;
delta_v = 0.2;

coluna = 1;

for u = piso_u : delta_u : teto_u
    linha = 1;
    for v = piso_v : delta_v : teto_v
        X(linha,coluna) = eval(f);
        Y(linha,coluna) = eval(g);
        Z(linha,coluna) = eval(h);
        linha = linha + 1;
    end
    coluna = coluna+1;
end

% =============================================================
% Calculo dos vetores de base (BUx, BUy, BUz, BVx, BVy, BVz)
% =============================================================
if (piso_u > 0)
  m = (piso_u + teto_u) / delta_u + 1;
else
  m = (teto_u - piso_u) / delta_u + 1;
end

if (piso_v > 0)
  n = (piso_v + teto_v) / delta_v + 1;
else
  n = (teto_v - piso_v) / delta_v + 1;
end

BUx = zeros(m,n);
BUy = zeros(m,n);
BUz = zeros(m,n);

BVx = zeros(m,n);
BVy = zeros(m,n);
BVz = zeros(m,n);

base1 = zeros(1,3);
base2 = zeros(1,3);

for i = 1 : 1 : m
    for j = 1 : 1 : n

        % Direcao para frente
        if ( i+1 < m )

            base1(1) = ( X(i+1,j) - X(i,j) ) / delta_u;
            base1(2) = ( Y(i+1,j) - Y(i,j) ) / delta_u;
            base1(3) = ( Z(i+1,j) - Z(i,j) ) / delta_u;

            BUx(i,j) = base1(1);
            BUy(i,j) = base1(2);
            BUz(i,j) = base1(3);
        end


        % Direcao e para baixo
        if ( j+1 < n )

            base1(1) = ( X(i,j+1) - X(i,j) ) / delta_v;
            base1(2) = ( Y(i,j+1) - Y(i,j) ) / delta_v;
            base1(3) = ( Z(i,j+1) - Z(i,j) ) / delta_v;

            BVx(i,j) = base1(1);
            BVy(i,j) = base1(2);
            BVz(i,j) = base1(3);
        end




    end
end

% =============================================================
% Calculo dos vetores normais (N1x, N1y, N1z, N2x, N2y, N2z)
% =============================================================
N1x = zeros(m,n);
N1y = zeros(m,n);
N1z = zeros(m,n);

VNB = zeros(m*n,3);

for i = 1 : 1 : m
    for j = 1 : 1 : n

        normal_total = zeros(1,3);
        
        % Direcao d = 0, pra frente
        if ( (i+1 < m) & (j+1 < n) )
            base1(1) = X(i+1,j) - X(i,j);
            base1(2) = Y(i+1,j) - Y(i,j);
            base1(3) = Z(i+1,j) - Z(i,j);

            base2(1) = X(i,j+1) - X(i,j);
            base2(2) = Y(i,j+1) - Y(i,j);
            base2(3) = Z(i,j+1) - Z(i,j);
        
            normal = cross(base1,base2);
            normal = normal./norm(normal);
            normal_total = normal_total + normal;
        end

        % Direcao d = 1, para baixo        
        if ( (j+1 < n) & (i-1 > 0) )
            base1(1) = X(i,j+1) - X(i,j);
            base1(2) = Y(i,j+1) - Y(i,j);
            base1(3) = Z(i,j+1) - Z(i,j);

            base2(1) = X(i-1,j) - X(i,j);
            base2(2) = Y(i-1,j) - Y(i,j);
            base2(3) = Z(i-1,j) - Z(i,j);
        
            normal = cross(base1,base2);
            normal = normal./norm(normal);
            normal_total = normal_total + normal;
        end
        
        % Direcao d = 2, para esquerda        
        if ( (i-1 > 0) & (j-1 > 0) )
            base1(1) = X(i-1,j) - X(i,j);
            base1(2) = Y(i-1,j) - Y(i,j);
            base1(3) = Z(i-1,j) - Z(i,j);

            base2(1) = X(i,j-1) - X(i,j);
            base2(2) = Y(i,j-1) - Y(i,j);
            base2(3) = Z(i,j-1) - Z(i,j);
        
            normal = cross(base1,base2);
            normal = normal./norm(normal);
            normal_total = normal_total + normal;            
        end
            
        % Direcao d = 3, para cima        
        if ( (j-1 > 0) & (i+1 < m) )
            base1(1) = X(i,j-1) - X(i,j);
            base1(2) = Y(i,j-1) - Y(i,j);
            base1(3) = Z(i,j-1) - Z(i,j);

            base2(1) = X(i+1,j) - X(i,j);
            base2(2) = Y(i+1,j) - Y(i,j);
            base2(3) = Z(i+1,j) - Z(i,j);
        
            normal = cross(base1,base2);
            normal = normal./norm(normal);
            normal_total = normal_total + normal;            
        end
       
        normal_total = normal_total./norm(normal_total);
        N1x(i,j) = normal_total(1);
        N1y(i,j) = normal_total(2);
        N1z(i,j) = normal_total(3);

        
        VNB( ( (i-1) * n + (j-1) + 1), 1) = N1x(i,j);
        VNB( ( (i-1) * n + (j-1) + 1), 2) = N1y(i,j);
        VNB( ( (i-1) * n + (j-1) + 1), 3) = N1z(i,j);
        
    end
    
end





% =============================================================
% Calculo dos tensores de metrica (TM00, TM01, TM10, TM11)
% =============================================================

TM00 = zeros(m,n);
TM01 = zeros(m,n);
Tm10 = zeros(m,n);
Tm11 = zeros(m,n);

for i = 1 : 1 : m
    for j = 1 : 1 : n

        base1 = [ BUx(i,j) BUy(i,j) BUz(i,j) ];
        base2 = [ BVx(i,j) BVy(i,j) BVz(i,j) ];

        TM00(i,j) = dot(base1,base1);
        TM01(i,j) = dot(base1,base2);
        TM10(i,j) = dot(base2,base1);
        TM11(i,j) = dot(base2,base2);
    end
end


% Calculo das derivadas segundas (DD00, DD01, DD10, DD11)

for i = 1 : 1 : m
    for j = 1 : 1 : n

        % Direcao para frente
        if ( (i-1 > 1) & (i+1 < m) )

            base1(1) = ( X(i+1,j) - X(i,j) ) / delta_u;
            base1(2) = ( Y(i+1,j) - Y(i,j) ) / delta_u;
            base1(3) = ( Z(i+1,j) - Z(i,j) ) / delta_u;

            BUx(i,j) = base1(1);
            BUy(i,j) = base1(2);
            BUz(i,j) = base1(3);
        end


        % Direcao e para baixo
        if ( j+1 < n )

            base1(1) = ( X(i,j+1) - X(i,j) ) / delta_v;
            base1(2) = ( Y(i,j+1) - Y(i,j) ) / delta_v;
            base1(3) = ( Z(i,j+1) - Z(i,j) ) / delta_v;

            BVx(i,j) = base1(1);
            BVy(i,j) = base1(2);
            BVz(i,j) = base1(3);
        end




    end
end







% Calculo dos tensores de curvatura (TC00, TC01, TC10, TC11)

TC00 = zeros(m,n);
TC01 = zeros(m,n);
TC10 = zeros(m,n);
TC11 = zeros(m,n);

for i = 1 : 1 : m
    for j = 1 : 1 : n

        base1 = [ BUx(i,j) BUy(i,j) BUz(i,j) ];
        base2 = [ BVx(i,j) BVy(i,j) BVz(i,j) ];

        TM00(i,j) = dot(base1,base1);
        TM01(i,j) = dot(base1,base2);
        TM10(i,j) = dot(base2,base1);
        TM11(i,j) = dot(base2,base2);
    end
end




% Desenho da superficie
% surf(X,Y,Z)
