% Plotting 3D harmonic oscillator orbital (or one of them), used in solution of exercise 7
close all;
% Quantum numbers =========================================================
   nx=1; ny=0; nz=0;
%==========================================================================
probability=0.19^2;
% hbar over m omega equal to 1
% Normalization
N=1;
psix=@(x) x.*exp(0.5-x.^2/2);
psiy=@(y) exp(-y.^2/2);
psiz=@(z) exp(-z.^2/2);
% Wave function
psi=@(x,y,z) psix(x).*psiy(y).*psiz(z);
% Setting the grid
border=3;      % must be increased for more excited orbitals
accuracy=100;
[x,y,z]=ndgrid(linspace(-border,border,accuracy),linspace(-border,border,accuracy),linspace(-border,border,accuracy));
% Conversion Cartesian to spherical coordinates
r=sqrt(x.^2+y.^2+z.^2);
theta=acos(z./r);
phi=atan2(y,x);
% Plot orbital,  - and + wave function phase
  colors=sign(psi(x,y,z));
isosurface(psi(x,y,z).^2,probability,colors);
title(['nx = ' num2str(nx) ', ny = ' num2str(ny) ', nz = ' num2str(nz)],'FontName','Times','FontSize',12);
set(gcf,'color',[1 1 1]);
daspect([1 1 1]);
grid;
axis off;
%view(3);
camlight('left');
camzoom(0.95);
lighting phong;
axis vis3d;
rotate3d on; 
view(90,0);      % change to view from different angles
brighten(1);
% print -dpng fig42155-7-4.png   % surface with psi/psi_max=sqrt(probability)
