
% kontur.m : contour curves for 2px orbital of hydrogenic atom (used for fig42155-8-7.png)
clear; format long; close all;

%f1 = figure(1);
 
bbstart=0.94;
xminstart = 2*( 1-sqrt(2*(1-bbstart)) ); 
xmaxstart = 2*( 1+sqrt(2*(1-bbstart)) );                    % initial suggestion                                                                                          %  for FZERO
hold on                      % because we run a loop creating several curves
N = 200; fakt = sqrt(N);
nN = -N:1:N; arg = nN/fakt;
b = bbstart : -0.10 : 0.015;                                     % psi/psi_max
for n = 1 : length(b) 
	  bb=b(n);                                                   % (psi/psi_max)^2
[xmin,eps] = fzero(@(x) 0.5*x*exp(1-x/2)-bb ,xminstart);  % min of x
% fprintf('bb = %10.5f, xmin = %10.5f\n ',bb,xmin);        
[xmax,eps] = fzero(@(x)  0.5*x*exp(1-x/2)-bb,xmaxstart); % max of x
% fprintf(' bb = %10.5f, xmax = %10.5f\n ',bb,xmax);
xminstart = xmin; xmaxstart = xmax;  % new suggestions for FZERO
xm = 0.5*(xmin+xmax);  dx = 0.5*(xmax-xmin); 
xx = xm + tanh(arg).*dx;               % using tanh to get more dense spacing 
                                       % near ends
x =[xmin xx xmax];                     % end points included
yplus = sqrt((2.*log(x/(2*bb))+2.001).^2-x.^2); % upper part of curve (z really)
yminus = -yplus;                       % lower part of curve
if n == 10 xxx = x; end          % chooses curve for psi/psi_max=b(n) 
                                % (for 2. diagram)
if n == 10 yyplus = yplus; end   %  - " -
 

plot(x,yplus,'r',x,yminus,'r','linewidth',2);

end % of b-loop
hold off
xlabel('x/a'); ylabel('z/a'); ax=10;
axis([0 2*ax -ax ax]); axis('square');
%%%%%%%%%%%%%%%%%%%%print -dpng fig42155-8-7.png

%%%%print -dpng fig42155-7-2.png % creates png-file for exercise (1st diagram)
f2 = figure(2);
plot(xxx,yyplus,'r',xxx,-yyplus,'r',-xxx,yyplus,'b',-xxx,-yyplus,'b','linewidth',2);
xlabel('x/a'); ylabel('z/a');
axis([-ax ax -ax ax]); axis equal;
%print -dpng fig42155-7-3.png   % full orbital with psi/psi_max=b(2) (2nd diagram)

% forget about what follows:
% f3 = figure(3);clf;
% z = [0: .02:1]'; r = sin(3*pi*z)+2;
%		    cylinder(r),axis square
		    
% f4 = figure(4);
%		    view(-30,30);
% z =[-xxx xxx]';
% r = [yyplus yyplus];
%		    cylinder(r),axis square;
