
% kontur.m : contour curves for orbital of 3D osc, used in text and solution of exercise 7.
clear; format long; close all;

f1 = figure(1);
xlabel('x/b'); ylabel('y/b');
axis([0 4.5 -2.25 2.25]); axis('square'); 
xminstart = 0.16; xmaxstart = 9;                            % initial suggestion                                                             %  for FZERO
hold on                      % because we run a loop creating several curves
N = 100; fakt = sqrt(N);
nN = -N:1:N; arg = nN/fakt;
b = 0.09 : 0.1 : 1;                                         % psi/psi_max
for n = 1 : length(b) 
	  a = b(n)^2;                                       % (psi/psi_max)^2
[xmin2,eps] = fzero(@(x) x*exp(-x+1.0000001)-a,xminstart);  % min of x^2
% fprintf(' a = %10.5f, xmin2 = %10.5f\n ',a,xmin2);        
[xmax2,eps] = fzero(@(x) x*exp(-x+1.0000001)-a,xmaxstart); % max of x^2
% fprintf(' a = %10.5f, xmax2 = %10.5f\n ',a,xmax2);
xminstart = xmin2; xmaxstart = xmax2;  % new suggestions for FZERO
xmin = sqrt(xmin2); 
xmax = sqrt(xmax2); 
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(log(x.^2/a)+1.0001-x.^2); % upper oart of curve
yminus = -yplus;                       % lower part of curve
if n == 2 xxx = x; end          % chooses curve for psi/psi_max=b(n) 
                                % (for 2. diagram)
if n == 2 yyplus = yplus; end   %  - " -
 

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

end % of b-loop
hold off

print -dpng fig42155-7-2.png % creates png-file for exercise (1st diagram)
f2 = figure(2);
axis([-2.5 2.5 -2 2]); axis('square');
plot(xxx,yyplus,'r',xxx,-yyplus,'r',-xxx,yyplus,'b',-xxx,-yyplus,'b','linewidth',2);
xlabel('x/b'); ylabel('y/b');
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;
