% --- Parámetros del problema ---
S0 = 50;      % kg/m2
alpha = 3;
L = 500;      % m
R = 1500;     % radio del embalse (m) -> círculo completo

% --- Mallado cartesiano (cuadrado) y máscara circular ---
n = 500;                                      % resolución (ajusta si necesitas más/menos)
x = linspace(-R, R, n);
y = linspace(-R, R, n);
[X, Y] = meshgrid(x, y);
inside = (X.^2 + Y.^2) <= R^2;                % máscara: verdadero dentro del círculo

% --- Campo S(x,y) y enmascarar el exterior ---
S = S0 * (1 + alpha * exp(-(X.^2 + Y.^2) / L^2));
S(~inside) = NaN;                             % fuera del círculo no se dibuja

% --- Gradiente analítico (dS/dx, dS/dy) y enmascarado ---
dSdx = S0 * alpha * exp(-(X.^2 + Y.^2)/L^2) .* (-2 .* X / L^2);
dSdy = S0 * alpha * exp(-(X.^2 + Y.^2)/L^2) .* (-2 .* Y / L^2);
dSdx(~inside) = NaN;
dSdy(~inside) = NaN;

% --- Figura con 3 subplots: 1 y 2 arriba, 3 ocupando toda la fila inferior ---
figure('Name','Sedimentación - mapas y campos','Units','normalized','Position',[0.05 0.05 0.9 0.8]);

% 1) Mapa de color (círculo completo)
subplot(2,2,1);
pcolor(X, Y, S);
shading interp;
axis equal tight;
xlabel('x (m)');
ylabel('y (m)');
title('sedimentación en el fondo','FontWeight','bold','FontSize',14);   % número arriba
colormap(gca, jet);
cb1 = colorbar;
cb1.Label.String = 'S (kg/m^2)';
set(gca,'XLim',[-R R],'YLim',[-R R]);


