Matlab equivalent of Maple densityplot

Question

I would like to create a density plot of this function:

In Maple, one could use the densityplot function to achieve this (code at the end), which gives:

However, I am not sure what to use for plotting a similar figure in MATLAB.

Here is my current MATLAB code:

x = [0:10:100];
y = [-50:10:50];
s = [10, 0];
i = [50,25];
for ii = 1 : length(x)
    sir(ii) = -10 * 9.8 * log10((power((x(ii) - s(1)),2) + power((y(ii) - s(2)),2)) / (power((x(ii) - i(1)),2) + power((y(ii) - i(2)),2)));  
end

Could someone suggest an equivalent in MATLAB?

For the density plot in Maple, I used

densityplot(sir(x,y), x=0..100, y=-50..50, axes=boxed, style=patchnogrid, scaletorange=-5..50, colorscheme = [black, "green", "white"])

Please see my edit, I've shown how to smooth out the colour banding and make it look more like your original (if you're interested!) — Wolfie
This seems to be more like the original. However I'm still missing the white spot with the light green color at S's position. Do you happen to know how I could achieve this? scaletorange=-5..50, colorscheme = [black, "green", "white"] did the trick in maple. — smyslov
What value of g are you using? If I use g=9.8 I get a bigger white spot (bigger patch where density > 50). If you answer that I'll edit my question with a fix to get the exact plot :) — Wolfie

Wolfie Wolfie · Accepted Answer · 2017-08-18T08:30:47

You can use surf (a 3D surface plot) to achieve this, but you will need a finer grid than steps of 10 for it to look good!

Also you will need meshgrid to get all combinations of the x and y coordinates.

Please see the comments for further details.

% Set up grid points
x = 0:0.1:100;
y = -50:0.1:50;
[x,y] = meshgrid(x,y);
% Set up parameters i, s and g
i = [50 25]; s = [10 0]; g = 9.8;
% Work out density  
% - no need for loop if we use element-wise operations ./ and .^
% - power(z,2) replaced by z.^2 (same function, more concise)
% - You forgot the sqare roots in your question's code, included using .^(1/2)
% - line continuation with "...", could remove and have on one line
sir = -10*g*log10( ((x-s(1)).^2 + (y-s(2)).^2).^(1/2) ./ ...
                   ((x-i(1)).^2 + (y-i(2)).^2).^(1/2)        ); 
% Plot, and set to a view from above
surf(x,y,sir,'edgecolor','none','facecolor','interp');
view(2);
% Change the colour scheme
colormap('bone')

Result:

Matching your example

You used the Maple command scaletorange=-5..50. This limits the scale between -5 and 50 (docs), so since sir is our scale variable, we should limit it the same. In MATLAB:

% Restrict sir to the range [-5,50]
sir = min(max(sir,-5),50);
% Of course we now have to replot
surf(x,y,sir,'edgecolor','none','facecolor','interp');
view(2);

Now, if you wanted the black/green colours, you can use a custom colormap, this would also smooth out the banding caused by the 'bone' colormap only having 64 colours.

% Define the three colours to interpolate between, and n interpolation points
black = [0 0 0]; green = [0 1 0]; white = [1 1 1]; 
n = 1000; 
% Do colour interpolation, equivalent to Maple's 'colorscheme = [black, "green", "white"]'
% We need an nx3 matrix of colours (columns R,G,B), which we get using interp1
colormap(interp1(1:3, [black; green; white], linspace(1,3,n)));

With g=3.5 (not sure what you used), we get an almost identical plot

Matlab equivalent of Maple densityplot

1 Answers