bh.m

Three dimensional hydrogen atom spherical harmonics visualization.

clear all;
clc;

N = 6;
L = 5;
M = 3;

if N <= L
    error('L i M nie moga byc wieksze badz rowne N');
end
if N < 0 
    error('N nie moze byc mniejsze zero');
end
if M < -L || M > L
    error('M nie moze byc spoza zakresu (-L .. +L)');
end

%nt = 240;
%np = 360;
nt = 60;
np = 90;

t1 = 0;
t2 = pi;
dt = (t2 - t1) / (nt - 1);
theta = [t1 : dt : t2];

p1 = 0;
p2 = 2 * pi;
dp = (p2 - p1) / (np - 1);
phi = [p1 : dp : p2];

%for N = 1 : 10
%N
%for r = 0.5 : 0.05 : 3
%    r
r = 10;   % r = wielokrotnosc ao

%rr = [0 : 50];
%RR = radial(N, L, rr);
%plot(rr, RR);
%figure;

sph = radial(N, L, r) * harmonics(L, M, cos(theta))' * cos(M * phi);
%sph = abs(radial(N, L, r) * hh(L, M, cos(theta))' * cos(M * phi));

x = sph .* (sin(theta)' * sin(phi));
y = sph .* (sin(theta)' * cos(phi));
z = sph .* (cos(theta)' * ones(size(phi)));

handle = surfl(x, y, z);
set(handle, 'FaceColor', [1 0 0]);
set(handle, 'EdgeColor', 'none');
set(handle, 'FaceAlpha', 0.6);
set(handle)
hold on;

sph2 = zeros(size(sph));
[w k] = size(sph);
for j = 1 : w
    for n = 1 : k
        if sph(j, n) <= 0
            sph2(j, n) = sph(j, n);
            sph(j, n) = 0;
        end
    end
end

sph2 = abs(sph2);

x = sph2 .* (sin(theta)' * sin(phi));
y = sph2 .* (sin(theta)' * cos(phi));
z = sph2 .* (cos(theta)' * ones(size(phi)));

handle = surfl(x, y, z);
set(handle, 'FaceColor', [0 1 0]);
set(handle, 'EdgeColor', 'none');
set(handle, 'FaceAlpha', 0.6);

light;
lighting phong;

%SHADING FLAT;

%plot(y)
%figure

%handle = surf(phi, theta, sph);
%set(handle, 'EdgeColor', 'none');
%set(handle, 'FaceAlpha', 0.6);
%axis('off')
%xlabel('phi')
%ylabel('theta')
%zlabel('Y')

%figure

%disp('x');
%size(x)
%disp('y');
%size(y)
%disp('z');
%size(z)

%handle = surf(x, y, z);
%title(['Funkcja falowa atomu H dla stanu [', num2str(N), ' ', num2str(L), ' ', num2str(M), '] w odleglosi r = ao * ', num2str(r)])
%axis('off');
%set(handle, 'EdgeColor', 'none');
%set(handle, 'FaceAlpha', 0.6);
%SHADING INTERP;
%view([1, 1, 0]);    

%xlim([-0.1 0.1])
%ylim([-0.1 0.1])
%zlim([-0.1 0.1])
 
%pause(0.05)
%end
%end