#_______________________________________________________________________________
#
# DESCOMPOSICIÓN POLAR EN 2D
#
#   Este programa realiza la descomposición polar de un tensor gradiente de
#     movimientos en 2D arbitrario.
#
#   Los efectos de las transformaciones geométricas se visualizan representando
#     cómo se deforma una circunferencia.
#_______________________________________________________________________________

F = input("Introduce el tensor F = [f11 f12; f21 f22] = ");

[nf, nc] = size(F);

if ((nf != 2) || (nc != 2))
  printf("Error: Solo se admiten matrices 2x2\n")
  return
end

printf("\n")

f = det(F)

if (f <= 0)
  printf("Error: el determinante de F = I + J tiene que ser positivo\n")
  return
end

# Geometría inicial (circunferencia de radio 1)
np = 400;
for ip = 1 : np + 1
  t = 2 * pi * ((ip - 1) / np);
  dx0(ip) = cos(t);
  dy0(ip) = sin(t);
endfor

# Geometría final
for ip = 1 : np +1
  dr0 = [dx0(ip) dy0(ip)]';            % vector material inicial

  dr  = F * dr0;                       % vector material final
  dx(ip) = dr(1);
  dy(ip) = dr(2);
endfor

# Escala de la ventana de dibujo
wmin = 1.1* min([ min(dx0), min(dx), min(dy0), min(dy)]);
wmax = 1.1* max([ max(dx0), max(dx), max(dy0), max(dy)]);

# Comparación de la geometría inicial y la geometría final
figure(1);

subplot(1,2,1);                        % geometría inicial
plot(dx0, dy0, "color", "b", "linewidth", 2);
axis("square", [wmin wmax wmin wmax]);
line([0 dx0(1)     ], [0, dy0(1)     ], "color", "g", "linewidth", 3);
line([0,dx0(np/4+1)], [0, dy0(np/4+1)], "color", "r", "linewidth", 3);
title("Geometría Inicial");

subplot(1,2,2);                        % geometría final
plot(dx, dy, "color", "b", "linewidth", 2);
axis("square", [wmin wmax wmin wmax]);
line([0 dx(1)     ], [0, dy(1)     ], "color", "g", "linewidth", 3);
line([0,dx(np/4+1)], [0, dy(np/4+1)], "color", "r", "linewidth", 3);
title("Geometría Final");


# Descomposicion Polar
C = F' * F                             % C = Tensor derecho de Cauchy-Green
[U, Lambda2] = eig(C)                  % Autovectores y autovalores de C
Lambda = sqrt(Lambda2)                 % Raíz cuadrada de los autovalores de C
I = eye(2)                             % Matrix identidad
E = U * (Lambda - I) * U'              % Tensor de deformación de Biot
R = F * U * inv(Lambda) * U'           % Tensor de rotación
Ed = f^(-1/2)* (E + (1 - f^(1/2)) * I) % Tensor de distorsión
In = f^(1/2) * I                       % Tensor de inflación

# Comprobación
F

R * (I + E)

R * (I + Ed) * In

# Cambios geométricos provocados por la Inflación, la distorsión y la rotación
for ip = 1 : np +1
  dr0 = [dx0(ip) dy0(ip)]';            % vector material inicial
  
  dr1 = In * dr0;                      % vector material tras la inflación
  dx1(ip) = dr1(1);
  dy1(ip) = dr1(2);
  
  dr2 = (I + Ed) * dr1;                % vector material tras la distorsión
  dx2(ip) = dr2(1);
  dy2(ip) = dr2(2);
  
  dr3 = R * dr2;                       % vector material tras la rotación
  dx3(ip) = dr3(1);
  dy3(ip) = dr3(2);
endfor



# Evolución de los cambios geométricos
figure(2);

subplot(1,4,1);                        % geometría inicial
plot(dx0, dy0, "color", "b", "linewidth", 2);
axis("square", [wmin wmax wmin wmax]);
line([0 dx0(1)     ], [0, dy0(1)     ], "color", "g", "linewidth", 3);
line([0,dx0(np/4+1)], [0, dy0(np/4+1)], "color", "r", "linewidth", 3);
title("Geometría Inicial");

subplot(1,4,2);                        % geometría tras la inflación
plot(dx1, dy1, "color", "b", "linewidth", 2);
axis("square", [wmin wmax wmin wmax]);
line([0 dx1(1)     ], [0, dy1(1)     ], "color", "g", "linewidth", 3);
line([0,dx1(np/4+1)], [0, dy1(np/4+1)], "color", "r", "linewidth", 3);
title("Geometría tras la INFLACIÓN");

subplot(1,4,3);                        % geometría tras la distorsión
plot(dx2, dy2, "color", "b", "linewidth", 2);
axis("square", [wmin wmax wmin wmax]);
line([0 dx2(1)     ], [0, dy2(1)     ], "color", "g", "linewidth", 3);
line([0,dx2(np/4+1)], [0, dy2(np/4+1)], "color", "r", "linewidth", 3);
title("Geometría tras la DISTORSIÓN");

subplot(1,4,4);                        % geometría tras la rotación
plot(dx3, dy3, "color", "b", "linewidth", 2);
axis("square", [wmin wmax wmin wmax]);
line([0 dx3(1)     ], [0, dy3(1)     ], "color", "g", "linewidth", 3);
line([0,dx3(np/4+1)], [0, dy3(np/4+1)], "color", "r", "linewidth", 3);
title("Geometría tras la ROTACIÓN");


