Empty spherical plot - strange error

I find it strange that I get an empty plot with the give command, and get the given error:

Error using matlab.graphics.chart.primitive.Surface
Value must be a vector or 2D array of numeric type.

Error in surf (line 145)
hh = matlab.graphics.chart.primitive.Surface(allargs{:});

Error in polar_coord_soln_Manz (line 59)
surf(X, Y, Z, Psi);

% Constants
hbar = 1.0545718e-34;  
m = 9.10938356e-31;    
E_ion = 5.139 * 1.60218e-19;   
k_f = 2 * m * E_ion / hbar^2;   
% Define alpha (renamed to avoid conflict with MATLAB function)
alpha_val = sqrt(k_f);
% Radial wave function
function R = radial_wavefunction(r, n, l, alpha)
L = laguerreL(n-l-1, 2*l+1, alpha * r.^2);
R = sqrt((2 * alpha)^(l+1) / factorial(n-l-1)) .* exp(-alpha * r.^2) .* (alpha * r).^l .* L;
end
% Spherical harmonic (assuming it's defined elsewhere)
function Y = spherical_harmonic(theta, phi, l, m)
Y = legendre(l, cos(theta)) .* exp(1i * m * phi);
end
% Total wave function in spherical coordinates
function psi = spherical_wavefunction(r, theta, phi, n, l, m, alpha)
R = radial_wavefunction(r, n, l, alpha);
Y = spherical_harmonic(theta, phi, l, m);
psi = R .* Y;
end
% Define grid
r = linspace(0, 10, 50);      % Radial coordinate r
theta = linspace(0, pi, 50);  % Polar angle theta
phi = linspace(0, 2*pi, 50);  % Azimuthal angle phi
% Create grid for 3D plotting
[R, Theta, Phi] = meshgrid(r, theta, phi);
 
n = 1;   
l = 0;   
m = 0;  
 
Psi = zeros(size(R));
 
for i = 1:numel(R)
    Psi(i) = abs(spherical_wavefunction(R(i), Theta(i), Phi(i), n, l, m, alpha_val))^2;  % Taking absolute value and squaring
end
% Reshape Psi to be 2D
Psi = reshape(Psi, size(R));
% Spherical to Cartesian Conversion
X = R .* sin(Theta) .* cos(Phi);
Y = R .* sin(Theta) .* sin(Phi);
Z = R .* cos(Theta);
% Plotting 3D surface
figure;
surf(X, Y, Z, Psi);
xlabel('x');
ylabel('y');
zlabel('z');
title(['|\psi_{', num2str(n), ',', num2str(l), ',', num2str(m), '}(r, \theta, \phi)|^2 for Sodium']);
colorbar;
axis equal;

Answer 1

[Manikanta Aditya]

Hi @Sergio,

It looks like you're encountering an error because surf expects X, Y, and Z to be 2D arrays, but you are passing 3D arrays. Additionally, the way you are computing Psi seems to be incorrect as you are using a 1D loop over a 3D grid.

Check this code that should fix the issue you encountered:

% Constants
hbar = 1.0545718e-34;  
m = 9.10938356e-31;    
E_ion = 5.139 * 1.60218e-19;   
k_f = 2 * m * E_ion / hbar^2;   

% Define alpha (renamed to avoid conflict with MATLAB function)
alpha_val = sqrt(k_f);

% Radial wave function
function R = radial_wavefunction(r, n, l, alpha)
    L = laguerreL(n-l-1, 2*l+1, alpha * r.^2);
    R = sqrt((2 * alpha)^(l+1) / factorial(n-l-1)) .* exp(-alpha * r.^2) .* (alpha * r).^l .* L;
end

% Spherical harmonic (assuming it's defined elsewhere)
function Y = spherical_harmonic(theta, phi, l, m)
    Y = legendre(l, cos(theta)) .* exp(1i * m * phi);
end

% Total wave function in spherical coordinates
function psi = spherical_wavefunction(r, theta, phi, n, l, m, alpha)
    R = radial_wavefunction(r, n, l, alpha);
    Y = spherical_harmonic(theta, phi, l, m);
    psi = R .* Y;
end

% Define grid
r = linspace(0, 10, 50);      % Radial coordinate r
theta = linspace(0, pi, 50);  % Polar angle theta
phi = linspace(0, 2*pi, 50);  % Azimuthal angle phi

% Create grid for 3D plotting
[R, Theta, Phi] = ndgrid(r, theta, phi);

n = 1;   
l = 0;   
m = 0;  

% Compute the wave function on the grid
Psi = abs(spherical_wavefunction(R, Theta, Phi, n, l, m, alpha_val)).^2;

% Spherical to Cartesian Conversion
X = R .* sin(Theta) .* cos(Phi);
Y = R .* sin(Theta) .* sin(Phi);
Z = R .* cos(Theta);

% Plotting 3D surface
figure;
surf(X(:,:,1), Y(:,:,1), Z(:,:,1), Psi(:,:,1)); % Plotting a slice of the 3D data
xlabel('x');
ylabel('y');
zlabel('z');
title(['|\psi_{', num2str(n), ',', num2str(l), ',', num2str(m), '}(r, \theta, \phi)|^2 for Sodium']);
colorbar;
axis equal;

Changed meshgrid to ndgrid to properly handle 3D grids. Used ndgrid to generate a 3D grid and compute the wave function on this grid. Since surf expects 2D inputs, I plotted a slice of the 3D data (Psi(:,:,1)).

I hope this helps!

[Sergio]

Thanks!! Is it possible to plot the 3D plot as some surface?

[Manikanta Aditya]

Hi, Yes we can do it.

To visualize the wave function as a 3D surface, we first use meshgrid to create a 2D grid of r and theta values while fixing phi to a constant value (e.g., pi/4). We then use surf to plot the 3D surface, setting 'EdgeColor', 'none' for a smoother appearance. The colormap('jet') function is applied to enhance the visibility of variations in the plot, and view(45, 30) is set to provide an optimal initial viewing angle.

Check this code:

% Constants
hbar = 1.0545718e-34;
m = 9.10938356e-31;
E_ion = 5.139 * 1.60218e-19;
k_f = 2 * m * E_ion / hbar^2;
alpha_val = sqrt(k_f);

% Radial wave function
function R = radial_wavefunction(r, n, l, alpha)
    L = laguerreL(n-l-1, 2*l+1, alpha * r.^2);
    R = sqrt((2 * alpha)^(l+1) / factorial(n-l-1)) .* exp(-alpha * r.^2) .* (alpha * r).^l .* L;
end

% Spherical harmonic (assuming it's defined elsewhere)
function Y = spherical_harmonic(theta, phi, l, m)
    Y = legendre(l, cos(theta)) .* exp(1i * m * phi);
end

% Total wave function in spherical coordinates
function psi = spherical_wavefunction(r, theta, phi, n, l, m, alpha)
    R = radial_wavefunction(r, n, l, alpha);
    Y = spherical_harmonic(theta, phi, l, m);
    psi = R .* Y;
end

% Define grid
[r, theta] = meshgrid(linspace(0, 10, 50), linspace(0, pi, 50));
phi = pi/4; % Fixed value for phi

% Set quantum numbers
n = 1;
l = 0;
m = 0;

% Compute the wave function on the grid
Psi = abs(spherical_wavefunction(r, theta, phi, n, l, m, alpha_val)).^2;

% Spherical to Cartesian Conversion
X = r .* sin(theta) .* cos(phi);
Y = r .* sin(theta) .* sin(phi);
Z = r .* cos(theta);

% Plotting 3D surface
figure;
surf(X, Y, Z, Psi, 'EdgeColor', 'none');
colormap('jet');
xlabel('x');
ylabel('y');
zlabel('z');
title(['|\psi_{', num2str(n), ',', num2str(l), ',', num2str(m), '}(r, \theta, \phi)|^2 for Sodium']);
colorbar;
axis equal;
view(45, 30);