I want to contour plot the data as shown in attached figure from my equation. I am not able to get the gaussian type behavior. What's wrong with my code?

%%%Simulation of Ion beam profile
M=40*1.67e-27;%%Mass of Argon in kg
Ib=40/1000;%%Beam Current in A
vb=sqrt((2*1000*1.6e-19)/(M));%%Velocity of Argon ion in m/s
R=100/1000;%%Initial Beam width in m
Rho=(M*Ib)/(2*pi*vb*R^2);  %%Space Charge density
h_t=linspace(100,500,100);
r=linspace(-100,100,500);

[H,R]=meshgrid(h_t,r);
Rho_m=(Rho./(H).^2).*exp(-(R.^2./(2.*R.^2.*(H).^2)));

figure(55)
contourf(R,H,Rho_m,50,'edgecolor','none')
shading interp
c=colorbar;
c.FontSize=12;
c.Label.String = 'J(\muA/cm^{2})';
% caxis([0 1200])
colormap jet
set (gca,'fontsize',10, 'fontweight', 'b')
xlabel('Radial distance (mm)');
ylabel('Axial distance (mm)')
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [8 8]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0.5 0.5 7 6]);
set(gcf, 'renderer', 'painters');
axis equal
print(gcf, '-dmeta', '1000V_51mA.emf')

Answer 1

[VBBV]

h_t=linspace(100,500,100);
r=linspace(-100,100,500);
[H,R]=meshgrid(h_t,r);
Rho_m=Ib.*(Rho./(H)).^2.*exp(-(2*R.^2./(H))); % Use  this for Gaussian dist

Check with this for Gaussian

[VBBV]

You may have to change the units to improve or get the exact Gaussian profile

[VBBV]

M=40*1.67e-27;%%Mass of Argon in kg
Ib=40/1000;%%Beam Current in A
vb=sqrt((2*1000*1.6e-19)/(M));%%Velocity of Argon ion in m/s
R=100/1000;%%Initial Beam width in m
Rho=(M*Ib)/(2*pi*vb*R^2);  %%Space Charge density
h_t=linspace(100,500,100);
r=linspace(-100,100,500);
[H,R]=meshgrid(h_t,r);
Rho_m=Ib.*(Rho./(H)).^2.*exp(-(R.^2./(H))); % Use  this for Gaussian dist
figure(55)
contourf(R,H,Rho_m,50,'edgecolor','none')
shading interp

[BHARAT SINGH RAWAT]

Thank you it was all what I needed!