[X, Y] = meshgrid(linspace(-1, 1, gridSize), linspace(-1, 1, gridSize));

Ex_r = X ./ sqrt(X.^2 + Y.^2);

Ey_r = Y ./ sqrt(X.^2 + Y.^2);

Ex_a = -Y ./ sqrt(X.^2 + Y.^2);

Ey_a = X ./ sqrt(X.^2 + Y.^2);

intensityImage = getIntensityImage(X, Y, gridSize, 0.5);

imagesc(unique(X), unique(Y), intensityImage, 'interp', 'bilinear');

quiver(X, Y, Ex_r, Ey_r, 0.5, "filled", "Color", "black");

title('Radial Polarization');

imagesc(unique(X), unique(Y), intensityImage, 'interp', 'bilinear');

quiver(X, Y, Ex_a, Ey_a, 0.5, "filled", "Color", "black");

title('Angular Polarization');

function image = getIntensityImage(X, Y, gridSize, sigma)

I = exp(-((X.^2 + Y.^2) / (2 * sigma^2)));

I_normalized = (I - min(I(:))) / (max(I(:)) - min(I(:)));

cmap = colormap(parula(256));

I_mapped = round(I_normalized * 255) + 1;

image = zeros(gridSize, gridSize, 3);

image(i, j, :) = cmap(I_mapped(i, j), :);