Chapter 8 M-Files

M-File Home

Script Files

ShowMySqrt Plots relative error associated with MySqrt.
ShowBisect Illustrates the method of bisection.
ShowNewton Illustrates the classical Newton iteration.
FindConj Uses fzero to compute Mercury-Earth conjunctions.
FindTet Applies fmin to three different objective functions.
ShowGolden Illustrates Golden section search.
ShowSD Steepest descent test environment.
ShowN Newton test environment.
ShowFDN Finite difference Newton test environment.
ShowMixed Globalized Newton test environment.
ShowGN Gauss-Newton test environment.
ShowFmin Illustrates fmin.
ShowFmins Illustrates fmins.

Function Files

MySqrt Canonical square root finder.
Bisection The method of bisection.
StepIsIn Checks next Newton step.
GlobalNewton Newton method with bisection globalizer.
GraphSearch Interactive graphical search environment.
Golden Golden section search.
SDStep Steepest descent step.
NStep Newton step.
FDNStep Finite difference Newton step.
GNStep Gauss-Newton step.
SineMercEarth The sine of the Mercury-Sun-Earth angle.
DistMercEarth Mercury-Earth distance.
TA Surface area of tetrahedron.
TV Volume of tetrahedron.
TE Total edge length of tetrahedron.
TEDiff Difference between tetrahedron edges.
Orbit Generates and plots orbit points.
Sep Separation between points on two orbits.
gSep Gradient of Sep.
gHSep Gradient and Hessian of sep.
SepV Vector from a point on one orbit to a point on another orbit.
JSepV Jacobian of SepV.
rho Residual of orbit fit.
Jrho Jacobian of rho.

 


ShowMySqrt

% Script File: ShowMySqrt
% Plots the error in the function MySqrt.

close all
L = input('Enter left endpoint of test interval:');
R = input('Enter right endpoint of test interval:');
Avals = linspace(L,R,300);
s = zeros(1,300);
for k=1:300
   s(k) = MySqrt(Avals(k));
end
error = abs(s-sqrt(Avals))./sqrt(Avals);
figure
plot(Avals,error+eps*.01)
Title('Relative Error in the Function MySqrt(A)')
xlabel('A');

ShowBisect

% Script File: ShowBisect
% Illustrates 6 steps of the method of bisection.

close all
fName   = input('Enter the name of the function (in quotes):');
a0      = input('Enter the left  endpoint of interval of interest:');
b0      = input('Enter the right endpoint of interval of interest:');

% Get initial interval.
x = linspace(a0,b0);
y = feval(fName,x);
figure
plot(x,y,x,0*x);
title('Click in left endpoint a.');
[a,y] = ginput(1);
title('Click in right endpoint b.');
[b,y] = ginput(1);
% Display the initial situation.
x  = linspace(a-.1*(b-a),b+.1*(b-a));
y  = feval(fName,x);
fa = feval(fName,a); 
fb = feval(fName,b);  
plot(x,y,x,0*x,[a b], [fa fb],'*')
xlabel(sprintf('At the start:    a = %10.6f   b = %10.6f',a,b));
for k=1:6
   mid = (a+b)/2;
   fmid = feval(fName,mid);  
   title('Strike any key to continue');  
   pause
   if fa*fmid<=0
      % There is a root in [a,mid].
      b  = mid; 
      fb = fmid;
   else
      % There is a root in [mid,b].
      a  = mid; 
      fa = fmid;
   end   
   x = linspace(a-.1*(b-a),b+.1*(b-a));
   y = feval(fName,x);
   figure
   plot(x,y,x,0*x,[a b], [fa fb],'*')
   xlabel(sprintf('After %1.0f steps:     a = %10.6f      b = %10.6f',k,a,b));
end 

ShowNewton

% Script File: ShowNewton
% A pure Newton method test environment.

close all
fName  = input('Enter the name of the function (with quotes);');
fpName = input('Enter the name of the derivative function (with quotes);');

a = input('Enter left endpoint of interval of interest:');
b = input('Enter right endpoint of interval of interest:');
xvals = linspace(a,b,100);
fvals = feval(fName,xvals);
figure
plot(xvals,fvals,xvals,0*xvals)
v = axis;
title('Click in Starting Value')
[xc,y] = ginput(1);
fc    = feval(fName,xc);
fpc   = feval(fpName,xc);
Lvals = fc + fpc*(xvals-xc);
hold on
plot(xvals,Lvals,xc,fc,'*')
axis(v)
xlabel(sprintf(' xc = %10.6f      fc = %10.3e',xc,fc))
hold off
for k=1:3
   step = -fc/fpc;
   xc = xc+step;
   a = xc-abs(step);
   b = xc+abs(step);
   xvals = linspace(a,b,100);
   fvals = feval(fName,xvals);
   fc    = feval(fName,xc);
   fpc   = feval(fpName,xc);
   Lvals = fc + fpc*(xvals-xc);
   figure
   plot(xvals,fvals,xvals,Lvals,xvals,0*xvals,xc,fc,'*')
   xlabel(sprintf(' xc = %10.6f      fc = %10.3e',xc,fc))
end 

 


FindConj

% Script File: FindConj
% Estimates spacing between Mercury-Earth Conjunctions

clc
GapEst = input('Spacing Estimate = ');
disp(sprintf('Next ten conjunctions, Spacing Estimate = %8.3f',GapEst))
disp(' ')
t = zeros(11,1);
disp('Conjunction    Time       Spacing')
disp('-----------------------------------')
for k=1:10;
   t(k+1) = fzero('SineMercEarth',k*GapEst);
   disp(sprintf('    %2.0f       %8.3f    %8.3f',k,t(k+1),t(k+1)-t(k)))
end

FindTet

% Script File: FindTet
% Applies fmin to three different objective functions.

tol = 10e-16;

reps = 20;
disp(' ')
disp('Objective Function     Time           Theta')
disp('----------------------------------------------------')

t0=clock;
for k=1:reps
   OptLat1 = fmin('TA',-pi/2,pi/2,[0 tol]);
end
t1 = etime(clock,t0);
disp(sprintf('      Area       %10.3f    %20.16f',1,OptLat1))
t0 = clock;
for k=1:reps
   OptLat2 = fmin('TV',-pi/2,pi/2,[0 tol]);
end
t2 = etime(clock,t0);
disp(sprintf('    Volume       %10.3f    %20.16f',t2/t1,OptLat2))
t0 = clock;
for k=1:reps
   OptLat3 = fmin('TE',-pi/2,pi/2,[0 tol]);
end
t3 = etime(clock,t0);
disp(sprintf('      Edge       %10.3f    %20.16f',t3/t1,OptLat3))
OptLat0 = fzero('TEdiff',-.3,eps);
disp(sprintf('\n                 fzero method: %20.16f',OptLat0))
disp(sprintf('                        Exact: %20.16f',asin(-1/3)))

 


ShowGolden

% Script File: ShowGolden
% Illustrates Golden Section Search

close all
tvals = linspace(900,950);
plot(tvals,DistMercEarth(tvals));
axis([900 950 40 150])
axis('off')
hold on
fname = 'DistMercEarth';
a = 900;
b = 950;
r = (3 - sqrt(5))/2;
c = a + r*(b-a);     fc = feval(fname,c);
d = a + (1-r)*(b-a); fd = feval(fname,d);
y = 78;
step = 0;
plot([900 950],[y y],[a c d b],[y y y y],'o')
text(952,y,num2str(step))
text(950.5,y+8,'Step')
title('Golden Section Search: (a c d b) Trace')
while step < 5
   if fc >= fd
      z = c + (1-r)*(b-c); 
      % [a c d b ] <--- [c d z b]
      a = c; 
      c = d; fc = fd;
      d = z; fd = feval(fname,z);
   else
      z = a + r*(d-a); 
      % [a c d b ] <--- [a z c d]
      b = d; 
      d = c; fd = fc;
      c = z; fc = feval(fname,z);
   end
   y = y - 7;
   step = step+1;
   plot([900 950],[y y],[a c d b],[y y y y],'o')
   text(952,y,num2str(step))
   
end
tmin = (c+d)/2;
hold off

ShowSD

% Script File: ShowSD
% Minimizes the sep function using steepest descent.

planet1 = struct('A',10,'P',2,'phi', pi/8); 
planet2 = struct('A', 4,'P',1,'phi',-pi/7);

close all
Lmax = 1;           % Max line step 
nSteps   = 25;      % Total number of steps.
manualSteps = 3;    % Number of manual steps with point-and-click line search.
con_size = 20;      % Size of the contour plot

% Plot the two orbits.
figure
t=linspace(0,2*pi);
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');
hold off

% Display contours of the sep function.
tRange = linspace(0,2*pi,con_size);
Z = zeros(con_size,con_size);
for i=1:con_size
   for j=1:con_size
      t = [tRange(i);tRange(j)];
      Z(i,j) = Sep(t,planet1,planet2);
   end
end
figure
C = contour(tRange,tRange,Z',10);
title('Enter Contour Labels (Optional). To quit, strike .')
xlabel('t1')
ylabel('t2')
Clabel(C,'manual') 
% Get ready for the iteration.
title('Enter Starting Point for Steepest Descent')
tc = zeros(2,1);
[tc(1),tc(2)] = ginput(1);
fc = Sep(tc,planet1,planet2);
gc = gSep(tc,planet1,planet2);
tvals = tc
fvals = fc;
gvals = norm(gc);
title('')
% Steepest Descent with line search:

clc
disp('Step       sep            t(1)         t(2)       norm(grad)')
disp('----------------------------------------------------------------')
disp(sprintf('%2.0f     %10.6f     %10.6f   %10.6f      %8.3e ',0,fc,tc,norm(gc)))
figure(3)

for step = 1:nSteps
   if step<=manualSteps  
      % Manual Line Search
      [tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,0);
   else
      % Automated Line Search
      [tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,1);
   end
   tvals = [tvals tnew];       tc = tnew;
   fvals = [fvals fnew];       fc = fnew;
   gvals = [gvals norm(gnew)]; gc = gnew;
   disp(sprintf('%2.0f     %10.6f     %10.6f   %10.6f      %8.3e ',step,fc,tc,norm(gc)))
end

% Show solution on orbit plot:

figure(1)
hold on
last = nSteps+1; 
pt1 = Orbit(tvals(1,last),planet1,'*');
pt2 = Orbit(tvals(2,last),planet2,'*');
plot([pt1.x pt2.x],[pt1.y pt2.y])
title(sprintf('min Sep = %8.4f',fvals(last)))
hold off

% Show the descent path on the contour plot:

figure(2)
hold on
plot(tvals(1,:),tvals(2,:),tvals(1,1),tvals(2,1),'o')
hold off
title(sprintf('tmin = (%8.3f,%8.3f)   norm(gmin)= %8.4e',tvals(1,last),tvals(2,last),gvals(last)))

% Plot the descent of the sep and its gradient:

figure(3)
subplot(2,1,1)
plot(fvals)
title('Value of Sep')
xlabel('Iteration')
subplot(2,1,2)
semilogy(gvals)
title('Value of norm(gSep).')
xlabel('Iteration')

ShowN

% Script File: ShowN
% Newton method test environment.

close all
planet1 = struct('A',15,'P',2,'phi', pi/10); 
planet2 = struct('A',20,'P',3,'phi',-pi/8);

itmax = 10;
tol = 1e-14;

% Plot Orbits

figure
axis equal off
t = linspace(0,2*pi);
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');

% Enter Starting Points
tc = rand(2,1)*2*pi;
Orbit(tc(1),planet1,'*');
Orbit(tc(2),planet2,'*');
title('Initial Guess')
hold off

% Initializations
Fc =  SepV(tc,planet1,planet2);   r = norm(Fc);
Jc = JSepV(tc,planet1,planet2);   c = cond(Jc);
clc
disp('Iteration        tc(1)                  tc(2)             norm(Fc)   cond(Jc)')
disp('------------------------------------------------------------------------------') 
disp(sprintf('    %2.0f   %20.16f   %20.16f     %8.3e   %8.3e',0,tc(1),tc(2),r,c)) 

% The Newton Iteration

k=0;
while (ktol)
   % Take a step and display
   [tc,Fc,Jc] = NStep(tc,Fc,Jc,planet1,planet2);
   
   k=k+1;
   figure
   Orbit(t,planet1,'-');
   hold on
   Orbit(t,planet2,'--');
   Orbit(tc(1),planet1,'*');
   Orbit(tc(2),planet2,'*');
   hold off 
   r = norm(Fc);
   c = cond(Jc);
   title(sprintf(' Iteration = %2.0f   norm(Fc) = %8.3e',k,r))
   disp(sprintf('    %2.0f   %20.16f   %20.16f     %8.3e   %8.3e',k,tc(1),tc(2),r,c)) 
end
sol = Orbit(tc(1),planet1);
disp(sprintf('\n Intersection = ( %17.16f , %17.16f )',sol.x,sol.y))

 


ShowFDN

% Script File: ShowFDN
% Finite Difference Newton method test environment.

close all
planet1 = struct('A',15,'P',2,'phi', pi/10); 
planet2 = struct('A',20,'P',3,'phi',-pi/8);

itmax = 10;
tol = 1e-14;

% Plot Orbits

figure
axis equal off
t = linspace(0,2*pi);
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');

% Enter Starting Points
tc = rand(2,1)*2*pi;
Orbit(tc(1),planet1,'*');
Orbit(tc(2),planet2,'*');
title('Initial Guess')
hold off


% Initializations

Fc =  SepV(tc,planet1,planet2);   r = norm(Fc);
clc
disp('Iteration        tc(1)                  tc(2)             norm(Fc)')
disp('--------------------------------------------------------------------') 
disp(sprintf('    %2.0f   %20.16f   %20.16f     %8.3e',0,tc(1),tc(2))) 

% The FDNewton Iteration


k=0;
while (ktol)
   % Take a step and display
   [tc,Fc] = FDNStep(tc,Fc,planet1,planet2);
   
   k=k+1;
   figure
   Orbit(t,planet1,'-');
   hold on
   Orbit(t,planet2,'--');
   Orbit(tc(1),planet1,'*');
   Orbit(tc(2),planet2,'*');
   hold off 
   r = norm(Fc);
   title(sprintf(' Iteration = %2.0f   norm(Fc) = %8.3e',k,r))
   disp(sprintf('    %2.0f   %20.16f   %20.16f     %8.3e   %8.3e',k,tc(1),tc(2),r)) 
end
sol = Orbit(tc(1),planet1);
disp(sprintf('\n Intersection = ( %17.16f , %17.16f )',sol.x,sol.y))

ShowMixed

% Script File: ShowMixed
% Minimizes the sep function by doing a few steepest descent steps
% (with line search) and then a few full Newton steps.


close all
SDSteps = 2;  % The number of steepest descent steps.
NSteps = 10;  % The number of Newton steps.

planet1 = struct('A',10,'P',2,'phi', pi/8); 
planet2 = struct('A', 4,'P',1,'phi',-pi/7);

Lmax = 1;           % Max line step 
con_size = 20;      % Size of the contour plot
% Plot the two orbits.
figure
t=linspace(0,2*pi);
Orbit(t,planet1,'-');
hold on
Orbit(t,planet2,'--');
hold off
% Display contours of the sep function.
tRange = linspace(0,2*pi,con_size);
Z = zeros(con_size,con_size);
for i=1:con_size
   for j=1:con_size
      t = [tRange(j);tRange(i)];
      Z(i,j) = Sep(t,planet1,planet2);
   end
end
figure
C = contour(tRange,tRange,Z,10);
title('Enter Contour Labels (Optional). To quit, strike .')
xlabel('t1')
ylabel('t2')
Clabel(C,'manual') 
% Get ready for the iteration.
title('Enter Starting Point for Steepest Descent')
tc = zeros(2,1);
[tc(1),tc(2)] = ginput(1);
fc = Sep(tc,planet1,planet2);
gc = gSep(tc,planet1,planet2);
tvals = tc;
fvals = fc;
gvals = norm(gc);
title('')
% Steepest Descent with line search:
clc
disp('Step     sep              t(1)                 t(2)          norm(grad)')
disp('------------------------------------------------------------------------')
disp(sprintf('%2.0f %10.6f     %18.15f   %18.15f    %5.1e ',0,fc,tc,norm(gc)))
figure(3)
for step = 1:SDSteps    
   [tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,0);      
   tvals = [tvals tnew];       tc = tnew;
   fvals = [fvals fnew];       fc = fnew;
   gvals = [gvals norm(gnew)]; gc = gnew;
   disp(sprintf('%2.0f %10.6f     %18.15f   %18.15f    %5.1e ',step,fc,tc,norm(gc)))
end
disp('------------------------------------------------------------------------')
[gc,Hc] = gHSep(tc,planet1,planet2);
for step=SDSteps+1:SDSteps+NSteps;
   tc = tc - Hc\gc;
   [gc,Hc] = gHSep(tc,planet1,planet2);
   fc = Sep(tc,planet1,planet2);	  
   tvals = [tvals tc];       
   fvals = [fvals fc];      
   gvals = [gvals norm(gc)]; 
   disp(sprintf('%2.0f %10.6f     %18.15f   %18.15f    %5.1e ',step,fc,tc,norm(gc)))
end
% Show solution on orbit plot:
figure(1)
hold on
last = length(gvals); 
pt1 = Orbit(tvals(1,last),planet1,'*');
pt2 = Orbit(tvals(2,last),planet2,'*');
plot([pt1.x pt2.x],[pt1.y pt2.y])
title(sprintf('min Sep = %8.4f',fvals(last)))
hold off
% Show the descent path on the contour plot:

figure(2)
hold on
plot(tvals(1,:),tvals(2,:),tvals(1,1),tvals(2,1),'o')
hold off
title(sprintf('tmin = (%8.3f,%8.3f)   norm(gmin)= %8.4e',tvals(1,last),tvals(2,last),gvals(last)))
% Plot the descent of the sep and its gradient:
figure(4)
subplot(2,1,1)
plot(fvals)
title('Value of Sep')
xlabel('Iteration')
subplot(2,1,2)
semilogy(gvals)
title('Value of norm(gSep).')
xlabel('Iteration')

ShowGN

% Script File: ShowGN 
% Illustrates Gauss-Newton with line search applied to the problem of fitting a
% two-parameter ellipse to some noisy, nearly-elliptical data.

A0 = 5; P0 = 1;  % Parameters of the "true" ellipse.
del = .1;        % Noise factor.
m = 50;          % Number of points.

con_size = 20;        % Size of contour plot
kmax     = 25 ;       % Maximum number of GN steps.
Lmax     = 8;         % Maximum step length factor
tol      = .000001;   % Size of gradient for termination.

% Generate a noisy orbit
theta = linspace(0,2*pi,m)';
c = cos(linspace(0,2*pi,m)');
r = ((2*A0*P0)./(P0*(1-c) + A0*(1+c))).*(1 + del*randn(m,1));
plist = [r c];

% Plot the true orbit and the radius vector samples
close all
figure
axis equal
planet = struct('A',A0,'P',P0,'phi',0);

%Orbit(linspace(0,2*pi,200),planet,'-');
hold on
x = r.*c;
y = r.*sin(theta);
plot(x,y,'o')
hold off  

% Generate Contour Plot
range = linspace(.9*min(r),1.5*max(r),con_size);
Z = zeros(con_size,con_size);
for i=1:con_size
   for j=1:con_size
      f = rho([range(j);range(i)],plist);
      Z(i,j) = (f'*f)/2;
   end
end
figure
C = contour(range,range,Z,20);
title('Enter Contour Labels (Optional). To quit, strike .')
xlabel('A'); ylabel('P'); Clabel(C,'manual')

% Initialize for iteration
title('Enter Starting Point for Gauss-Newton')
zc = zeros(2,1);
[zc(1),zc(2)] = ginput(1); 
Jc = Jrho(zc,plist);
Fc =  rho(zc,plist);
wc = (Fc'*Fc)/2;
k=0;
gnorm = norm(Jc'*Fc);

clc
disp(sprintf('"True A" = %5.3f,  "True P" = %5.3f',A0,P0))
disp(sprintf('Relative Noise = %5.3f,  Samples = %3.0f\n',del,m))
disp('Iteration     wc         A          P          norm(grad)')
disp('----------------------------------------------------------')
disp(sprintf('   %2.0f   %10.4f %10.6f %10.6f       %5.3e ',k,wc,zc(1),zc(2),gnorm))
figure(3)

% Gauss-Newton with linesearch
while (gnorm > tol ) & (k < kmax)
   if k<=2
      [zc,Fc,wc,Jc] = GNStep('rho','Jrho',zc,Fc,wc,Jc,plist,Lmax,0);
   else 
      [zc,Fc,wc,Jc] = GNStep('rho','Jrho',zc,Fc,wc,Jc,plist,Lmax,1);
   end
   k=k+1;
   gnorm = norm(Jc'*Fc);
   disp(sprintf('   %2.0f   %10.4f %10.6f %10.6f       %5.3e ',k,wc,zc(1),zc(2),gnorm))
end
figure(1)
planetModel = struct('A',zc(1),'P',zc(2),'phi',0);
hold on
Orbit(linspace(0,2*pi),planetModel,'-');
hold off

ShowFmin

% Script File: ShowFmin
% Illustrates the 1-d minimizer fmin.

clc
tol = .000001;
disp('Local minima of the Mercury-Earth separation function.')
disp(sprintf('\ntol = %8.3e\n\n',tol))
disp(' Initial Interval     tmin        f(tmin)    f evals')
disp('----------------------------------------------------')
options = zeros(18,1);
for k=1:8
   L = 100+(k-1)*112;
   R = L+112;
   options(2) = tol;
   [tmin options] = fmin('DistMercEarth',L,R,options);
   minfeval = options(8);
   nfevals = options(10);
   disp(sprintf('     [%3.0f,%3.0f]     %10.5f   %10.5f, %6.0f',L,R,tmin,minfeval,nfevals))
end

ShowFmins

% Script Files: ShowFmins
% Illustrates the multidimensional minimizer fmins.

close all
planet1 = struct('A',10,'P',2,'phi', pi/8); 
planet2 = struct('A', 4,'P',1,'phi',-pi/7);
% Plot the orbots
tVals = linspace(0,2*pi);
Orbit(tVals,planet1,'-');
hold on
Orbit(tVals,planet2,'--');

% Call fmins
trace = 0;
steptol = .000001;
ftol = .000001;
options = [trace steptol ftol];
t0=[5;2];
[t,options] = Fmins('Sep',t0,options,[],planet1,planet2);

% Show Solution
pt1 = Orbit(t(1),planet1,'o');
pt2 = Orbit(t(2),planet2,'o');
plot([pt1.x pt2.x],[pt1.y pt2.y])
hold off
its = options(10);
gval = norm(gSep(t,planet1,planet2));
title(sprintf('t_{*} = ( %8.6f , %8.6f ), Steps = %3.0f,  norm(gradient) = %5.3e',t,its,gval))
axis off equal

MySqrt

  function x = MySqrt(A)
% x = MySqrt(A)
% A is a non-negative real and x is its square root.

if A==0
   x = 0;
else
   TwoPower = 1;
   m = A;
   while m>=1
      m = m/4;
      TwoPower = 2*TwoPower;
   end
   while m < .25
      m = m*4;
      TwoPower = TwoPower/2;
   end	  
   % sqrt(A) = sqrt(m)*TwoPower	  
   x = (1+2*m)/3;
   for k=1:4
      x = (x + (m/x))/2;
   end
   x = x*TwoPower;
end

Bisection


  function root = Bisection(fname,a,b,delta)
% root = Bisection(fname,a,b,delta)
% Method of bisection.
%
% fname is a string that names a continuous function  f(x) of a single variable.
% a and b define an interval [a,b] and f(a)f(b) < 0
% delta   non-negative real.
%
% root is the midpoint of an interval [alpha,beta]
% with the property that f(alpha)f(beta)<=0 and
%            |
%             |beta-alpha| <= delta+eps*max(|alpha|,|beta|)

fa = feval(fname,a); 
fb = feval(fname,b);
if fa*fb > 0
   disp('Initial interval is not bracketing.')
   return
end
if nargin==3
   delta = 0;
end
while abs(a-b) > delta+eps*max(abs(a),abs(b))
   mid = (a+b)/2;
   fmid = feval(fname,mid);    
   if fa*fmid<=0
      % There is a root in [a,mid].
      b  = mid; 
      fb = fmid;
   else
      % There is a root in [mid,b].
      a  = mid; 
      fa = fmid;
   end  
end
root = (a+b)/2;

StepIsIn

  function ok = StepIsIn(x,fx,fpx,a,b)
% ok = StepIsIn(x,fx,fpx,a,b)
% Yields 1 if the next Newton iterate is in [a,b] and 0 otherwise.  
% x is the current iterate, fx is the value of f  at x, and fpx is 
% the value of f' at x.

if fpx > 0
   ok = ((a-x)*fpx <= -fx) & (-fx <= (b-x)*fpx);
elseif fpx < 0
   ok = ((a-x)*fpx >= -fx) & (-fx >= (b-x)*fpx);
else
   ok = 0;
end

GlobalNewton

  function [x,fx,nEvals,aF,bF] = GlobalNewton(fName,fpName,a,b,tolx,tolf,nEvalsMax)
% [ x,fx,nEvals,aF,bF] = GlobalNewton(fName,fpName,a,b,tolx,tolf,nEvalsMax)
% fName       string that names a function f(x).
% fpName      string that names the derivative function f'(x).
% a,b         A root of f(x) is sought in the interval [a,b]
%               and f(a)*f(b)<=0.
% tolx,tolf   Nonnegative termination criteria.
% nEvalsMax   Maximum number of derivative evaluations.
%
% x          An approximate zero of f.
% fx         The value of f at x.
% nEvals     The number of derivative evaluations required. 
% aF,bF      The final bracketing interval is [aF,bF].
%
% The iteration terminates as soon as x is within tolx of a true zero or
% if |f(x)|<= tolf or after nEvalMax f-evaluations.

fa  = feval(fName,a);
fb  = feval(fName,b);
if fa*fb>0
   disp('Initial interval not bracketing.')
   return
end
x   = a;
fx  = feval(fName,x);
fpx = feval(fpName,x);
disp(sprintf('%20.15f  %20.15f    %20.15f',a,x,b))

nEvals = 1;
while (abs(a-b) > tolx ) & (abs(fx) > tolf) & ((nEvals < nEvalsMax) | (nEvals==1))
   %[a,b] brackets a root and x = a or x = b.
   if StepIsIn(x,fx,fpx,a,b)
      %Take Newton Step
      disp('Newton')
      x   = x-fx/fpx;
   else
      %Take a Bisection Step:
      disp('Bisection')
      x = (a+b)/2;
   end
   fx  = feval(fName,x);
   fpx = feval(fpName,x); 
   nEvals = nEvals+1;  
   if fa*fx<=0
      % There is a root in [a,x]. Bring in right endpoint.
      b  = x; 
      fb = fx;
   else
      % There is a root in [x,b]. Bring in left endpoint.
      a  = x; 
      fa = fx;
   end 
   disp(sprintf('%20.15f  %20.15f    %20.15f',a,x,b))
end
aF = a;
bF = b;

GraphSearch

  function [L,R] = GraphSearch(fname,a,b,Save,nfevals)
% [L,R] = GraphSearch(fname,a,b,Save,nfevals)
%
% Graphical search. Produces sequence of plots of the function f(x). The user
% specifies the x-ranges by mouseclicks. 
%
% name is a string that names a function f(x) that is defined 
% on the interval [a,b]. 
% nfevals>=2 
%
% Save is used to determine how the plots are saved. If Save is
% nonzero, then each plot is saved in a separate figure window.
% If Save is zero or if GraphSearch is called with just three
% arguments, then only the final plot is saved. 
%
% [L,R] is the x-range of the final plot. The plots are based on nfevals
% function evaluations. If GraphSearch is called with less
% than five arguments, then nfevals is set to 100.


close all
if nargin==3
   Save=0;
end
if nargin&lt5
   nfevals=100;
end
AnotherPlot = 1;
L = a;
R = b;
while AnotherPlot
   x = linspace(L,R,nfevals);
   y = feval(fname,x);
   if Save
      figure
   end
   ymin = min(y);
   plot(x,y,[L R],[ymin ymin])
   title(['The Function ' fname])
   v = axis;
   v(1) = L;
   v(2) = R;
   v(3) = v(3)-(v(4)-v(3))/10;
   axis(v)
   xlabel('Enter New Left Endpoint. (Click off x-range to terminate.)')
   text(R,ymin,['  ' num2str(ymin)])
   [x1,y1] = ginput(1);
   if (x1 < L) | (R < x1) 	  
      AnotherPlot=0;
   else
      xlabel('Enter New Right Endpoint. (Click off x-range to terminate.)')
      [x2,y2] = ginput(1);
      if (x2 < L) | (R < x2) 
         AnotherPlot=0; 
      else
         L = x1;
         R = x2;
      end
   end
   xlabel(' ')
end

Golden

  function tmin = Golden(fname,a,b)
% tmin = Golden(fname,a,b)
% Golden Section Search
%
% fname   string that names function  f(t) of a single variable.
% a,b   define an interval [a,b] upon which f is unimodal.
%
% tmin   approximate global minimizer of f on [a,b].

r = (3 - sqrt(5))/2;
c = a + r*(b-a);     fc = feval(fname,c);
d = a + (1-r)*(b-a); fd = feval(fname,d);
while (d-c) > sqrt(eps)*max(abs(c),abs(d))
   if fc >= fd
      z = c + (1-r)*(b-c); 
      % [a c d b ] <--- [c d z b]
      a = c; 
      c = d; fc = fd;
      d = z; fd = feval(fname,z);
   else
      z = a + r*(d-a); 
      % [a c d b ] <--- [a z c d]
      b = d; 
      d = c; fd = fc;
      c = z; fc = feval(fname,z);
   end
end
tmin = (c+d)/2;

SDStep

  function [tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,auto)
% [tnew,fnew,gnew] = SDStep(tc,fc,gc,planet1,planet2,Lmax,auto)
% Generates a steepest descent step.
%
%  fName     the name of a function f(x,plist) that accepts column
%            n-vectors and returns a scalar.
%  gName     the name of a function g(x,plist) that returns the
%            gradient of f at x.
%     xc     a column-vector, an approximate minimizer
%     fc     f(xc,plist)
%     gc     g(xc,plist)
%  plist     orbit parameters for f and g
%   Lmax     max step length
%   auto     if auto is nonzero, then automated line search
%
%   xnew     a column n-vector, an improved minimizer
%   fnew     f(xnew,plist)
%   gnew     g(new,plist)
%

nplotvals = 20;

% Get the Steepest Descent Search Direction
sc = -gc;
% Line Search
% Try to get L<=Lmax so xc+L*sc is at least as good as xc.
L = Lmax;
Halvings = 0;
fL = Sep(tc+L*sc,planet1,planet2);
while (fL>=fc) & Halvings<=10
   L = L/2;
   Halvings = Halvings+1;
   fL = Sep(tc+L*sc,planet1,planet2);
end  
% Sample across [0.L]
lambdavals = linspace(0,L,nplotvals);
fvals = zeros(1,nplotvals);
for k=1:nplotvals
   fvals(k) = Sep(tc+lambdavals(k)*sc,planet1,planet2);
end
if auto==0
   % Manual line search.
   plot(lambdavals,fvals);
   xlabel('lambda');
   ylabel('Sep(tc+lambda*sc)');
   title('Click the Best lambda value.');
   [lambda,y] = ginput(1);
   tnew = tc+lambda*sc;
   fnew = Sep(tnew,planet1,planet2);
   gnew = gSep(tnew,planet1,planet2);
else
   % Automated line search.
   [fnew,i] = min(fvals(1:nplotvals));
   tnew = tc + lambdavals(i(1))*sc;
   gnew = gSep(tnew,planet1,planet2);
end 

NStep

  function [tnew,Fnew,Jnew] = NStep(tc,Fc,Jc,planet1,planet2)
% [tnew,Fnew,Jnew] = NStep(tc,Fc,Jc,planet1,planet2)
% Newton Step
%
% tc is a column 2-vector, Fc is the value of SepV(t,planet1,planet2) at t=tc
% and Jc is the value of JSepV(t,planet1,planet2) at t=tc. Does one Newton step
% applied to SepV rendering a new approximate root tnew. Fnew and Jnew are the 
% values of SepV(tnew,planet1,planet2) and jSepV(tnew,planet1,planet2) respectively.

sc = -(Jc\Fc);
tnew = tc + sc;
Fnew = SepV(tnew,planet1,planet2);   
Jnew = JSepV(tnew,planet1,planet2);

FDNStep

  function [tnew,Fnew] = FDNStep(tc,Fc,planet1,planet2)
% [tnew,Fnew] = FDNStep(tc,Fc,planet1,planet2)
% Finite difference Newton step.
% 
% tc is a column 2-vector, Fc is the value of SepV(t,planet1,planet2) at t=tc
% Does one finite difference Newton step applied to SepV rendering a new 
% approximate root tnew. Fnew is the values of SepV(tnew,planet1,planet2).

% Set up the FD Jacobian.
Jc = zeros(2,2);
delta = sqrt(eps);
for k=1:2
   tk = tc;
   tk(k) = tc(k) + delta;
   Jc(:,k) = (SepV(tk,planet1,planet2) - Fc)/delta;
end

% The Finite Difference Newton Step.
sc = -(Jc\Fc);
tnew = tc + sc;
Fnew = SepV(tnew,planet1,planet2);

GNStep

  function [xnew,Fnew,wnew,Jnew] = GNStep(FName,JName,xc,Fc,wc,Jc,plist,Lmax,auto)
% [xnew,Fnew,wnew,Jnew] = GNStep(FName,JName,xc,Fc,wc,Jc,plist,Lmax,auto)
% Generates a Gauss-Newton Step
% 
%  FName     the name of a function F(x,plist) that accepts  a column
%            n-vectors and column m-vector
%  JName     the name of a function JF(x,plist) that returns the
%            Jacobian of F at x.
%     xc     a column n-vector, an approximate minimizer.         
%     Fc     F(xc,plist)
%     wc     Fc'*Fc/2
%     Jc     JF(xc,plist)
%  plist     parameter list
%   auto     0 for manual line search and nonzero otherwise.
%
%   xnew     a column n-vector, an improved minimizer.
%   Fnew     F(xnew,plist)
%   wnew     Fnew'*Fnew/2
%   Jnew     JF(xnew,plist)

   nplotvals = 20; % Number of F evaluations per line search.

% Get the Gauss-Newton  Search Direction

sc = -Jc\Fc;

% Line Search

% Try to get L<=Lmax so xc+L*sc is at least as good as xc.

L = Lmax;
Halvings = 0;
FL = feval(FName,xc+L*sc,plist);
while ((FL'*FL/2)>=wc) & Halvings<=10
   L = L/2;
   Halvings = Halvings+1;
   FL = feval(FName,xc+L*sc,plist);
end  
lambdavals = linspace(0,L,nplotvals);

wvals = zeros(nplotvals,1);
for k=1:nplotvals
   F = feval(FName,xc+lambdavals(k)*sc,plist);
   wvals(k) = (F'*F)/2;
end
if auto==0
   % Manual line search
   plot(lambdavals,wvals);
   xlabel('lambda');
   ylabel('w(xc+lambda*sc)');
   title('Enter the Best lambda value.');
   [lambda,y] = ginput(1);
   xnew = xc+lambda*sc;  
else
   % Automated line search
   
   [fnew,i] = min(wvals(1:nplotvals));
   xnew = xc + lambdavals(i(1))*sc;
end
Fnew = feval(FName,xnew,plist); 
wnew = .5*(Fnew'*Fnew);
Jnew = feval(JName,xnew,plist);

SineMercEarth

  function s = SineMercEarth(t)
% s = SineMercEarth(t)
% The sine of the Mercury-Sun-Earth angle at time t

% Mercury location:
xm = -11.9084 +  57.9117*cos(2*pi*t/87.97);
ym =             56.6741*sin(2*pi*t/87.97);

% Earth location:
xe =  -2.4987 + 149.6041*cos(2*pi*t/365.25);
ye =            149.5832*sin(2*pi*t/365.25);

s = (xm.*ye - xe.*ym)./(sqrt(xm.^2 + ym.^2).*sqrt(xe.^2 + ye.^2));

DistMercEarth

  function s = DistMercEarth(t)
% s = DistMercEarth(t)
% The distance between Mercury and Earth at time t.

% Mercury location:
xm = -11.9084 +  57.9117*cos(2*pi*t/87.97);
ym =             56.6741*sin(2*pi*t/87.97);

% Earth location:
xe =  -2.4987 + 149.6041*cos(2*pi*t/365.25);
ye =            149.5832*sin(2*pi*t/365.25);

s = sqrt((xe-xm).^2 + (ye-ym).^2);

TA

  function A = TA(theta)
% A = TA(theta)
%
% theta is a real number in the interval [-pi/2,pi/2].
% A is a negative multiple of the surface area of a tetrahedron
% that is inscribed in the unit sphere. One vertex is at the north 
% pole (0,0,1) and the other three  form an equilateral triangle in the plane 
% z = sin(theta).  
 
c = cos(theta); 
s = sin(theta); 

A = -c.*(sqrt((3*s-5).*(s-1)) + c);

TV

  function V = TV(theta)
% V = TV(theta)
% 
% theta  real in interval [-pi/2,pi/2]
% V is a negative multiple of the volume of a tetrahedron
% that is inscribed in the unit sphere. One vertex is at the north 
% pole (0,0,1) and the other three form an equilateral triangle in the plane 
% z = sin(theta).  
  
s = sin(theta);     
V = (1-s.^2).*(s-1);

 


TE

  function E = TE(theta)
% E = TE(theta)
%
% theta is  a real number in the interval [-pi/2,pi/2]
% E is a  negative multiple of the total edge length of a tetrahedron
% that is inscribed in the unit sphere. One vertex is at the north 
% pole (0,0,1) and the other three form an equilateral triangle in the plane 
% z = sin(theta).  

E = -(sqrt((1-sin(theta))) + sqrt(3/2)*cos(theta));

TEDiff

  function E = TEDiff(theta)
% E = TEDiff(theta)
% 
% theta is a real number in the interval [-pi/2,pi/2]
%
% Let T be the tetrahedron with one vertex at the north pole (0.0,1) and the other 
% three forming an equilateral triangle in the plane  z = sin(theta). E is difference
% between a pole to base edge length and the length of a base edge. 
  
E = sqrt((1-sin(theta))) - sqrt(3/2)*cos(theta);

Orbit

  function pLoc = Orbit(t,planet,lineStyle)
% pLoc = Orbit(t,planet,lineStyle)
%
% t is a row vector and for k=1:length(t), pLoc.x(k) = x(t(k)) and 
% pLoc.y(k) = y(t(k)) where
% 
%  x(tau)    cos(phi) sin(phi)   (planet.A-planet.P)/2 + ((planet.A+planet.P)/2)cos(tau)
%         =                    *
%  y(tau)   -sin(phi) cos(phi)             sqrt(planet.A*planet.P)sin(tau)
%
% If nargin==3 then the points are plotted with line style defined by the
% string lineStyle.

c = cos(t); s = sin(t);
x0 = ((planet.P-planet.A)/2) + ((planet.P+planet.A)/2)*c;
y0 = sqrt(planet.A*planet.P)*s;
cphi = cos(planet.phi); sphi = sin(planet.phi);
pLoc = struct('x',cphi*x0 + sphi*y0,'y',-sphi*x0 + cphi*y0);
if nargin==3
   plot(pLoc.x,pLoc.y,lineStyle)
end

Sep

  function d = Sep(t,planet1,planet2)
% d = Sep(t,planet1,planet2)
%
% t is a 2-vector and planet1 and planet2 are orbit structures.
% t(1) defines a planet1 orbit point, t(2) defines a planet2 orbit point,
% and d is the distance between them.


pLoc1 = Orbit(t(1),planet1);
pLoc2 = Orbit(t(2),planet2);
d = ((pLoc1.x-pLoc2.x)^2 + (pLoc1.y-pLoc2.y)^2)/2;

 


gSep

  function g = gSep(t,planet1,planet2)
% g = gSep(t,planet1,planet2)
% The gradient of sep(t,planet1,planet2) with respect to 2-vector t.

A1 = planet1.A; P1 = planet1.P; phi1 = planet1.phi; 
A2 = planet2.A; P2 = planet2.P; phi2 = planet2.phi; 

alfa1  = (P1-A1)/2; beta1  = (P1+A1)/2; gamma1 = sqrt(P1*A1);
alfa2  = (P2-A2)/2; beta2  = (P2+A2)/2; gamma2 = sqrt(P2*A2);
s1     = sin(t(1)); c1     = cos(t(1));
s2     = sin(t(2)); c2     = cos(t(2));
cphi1  = cos(phi1); sphi1  = sin(phi1);
cphi2  = cos(phi2); sphi2  = sin(phi2);
Rot1   = [cphi1 sphi1; -sphi1 cphi1];
Rot2   = [cphi2 sphi2; -sphi2 cphi2];
P1     = Rot1*[alfa1+beta1*c1;gamma1*s1];
P2     = Rot2*[alfa2+beta2*c2;gamma2*s2];
dP1    = Rot1*[-beta1*s1;gamma1*c1];
dP2    = Rot2*[-beta2*s2;gamma2*c2]; 
g = [-dP1';dP2']*(P2-P1);

gHSep

  function [g,H] = gHSep(t,planet1,planet2)
% [g,H] = gHSep(t,planet1,planet2)
%
% t is a 2-vector and planet1 and planet2 are orbits.
% g is the gradient of Sep(t,planet1,planet2) and H is the Hessian.


A1 = planet1.A; P1 = planet1.P; phi1 = planet1.phi; 
A2 = planet2.A; P2 = planet2.P; phi2 = planet2.phi; 

alfa1  = (P1-A1)/2; beta1  = (P1+A1)/2; gamma1 = sqrt(P1*A1);
alfa2  = (P2-A2)/2; beta2  = (P2+A2)/2; gamma2 = sqrt(P2*A2);
s1     = sin(t(1)); c1     = cos(t(1));
s2     = sin(t(2)); c2     = cos(t(2));
cphi1  = cos(phi1); sphi1  = sin(phi1);
cphi2  = cos(phi2); sphi2  = sin(phi2);
Rot1   = [cphi1 sphi1; -sphi1 cphi1];
Rot2   = [cphi2 sphi2; -sphi2 cphi2];
P1     = Rot1*[alfa1+beta1*c1;gamma1*s1];
P2     = Rot2*[alfa2+beta2*c2;gamma2*s2];
dP1    = Rot1*[-beta1*s1;gamma1*c1];
dP2    = Rot2*[-beta2*s2;gamma2*c2];
g = [-dP1';dP2']*(P2-P1);
ddP1   = Rot1*[-beta1*c1;-gamma1*s1];
ddP2   = Rot2*[-beta2*c2;-gamma2*s2];
H = zeros(2,2);
H(1,1) = (P1(1)-P2(1))*ddP1(1) + (P1(2)-P2(2))*ddP1(2) + ...
   dP1(1)^2 + dP1(2)^2;
H(2,2) = (P2(1)-P1(1))*ddP2(1) + (P2(2)-P1(2))*ddP2(2) + ...
   dP2(1)^2 + dP2(2)^2;									    
H(1,2) = -dP1(1)*dP2(1) - dP1(2)*dP2(2);
H(2,1) = H(1,2);

 


SepV

  function z = SepV(t,planet1,planet2)
% z = SepV(t,planet1,planet2)
% The vector from the t(1) point on the planet1 orbit
% to the t(2) point on the planet2 orbit.

pt1 = Orbit(t(1),planet1);
pt2 = Orbit(t(2),planet2);
z = [pt2.x-pt1.x;pt2.y-pt1.y];

JSepV

  function J = JSepV(t,planet1,planet2) 
% J = JSepV(t,planet1,planet2)
% J is the Jacobian of sepV(t,planet1,planet2).

A1 = planet1.A; P1 = planet1.P; phi1 = planet1.phi; 
A2 = planet2.A; P2 = planet2.P; phi2 = planet2.phi; 

s1 = sin(t(1)); c1 = cos(t(1));
s2 = sin(t(2)); c2 = cos(t(2));
beta1 = (P1+A1)/2; gamma1 = sqrt(P1*A1);
beta2 = (P2+A2)/2; gamma2 = sqrt(P2*A2);
cphi1 = cos(phi1); sphi1  = sin(phi1);
cphi2 = cos(phi2); sphi2  = sin(phi2);
Rot1  = [cphi1 sphi1; - sphi1 cphi1];
Rot2  = [cphi2 sphi2; - sphi2 cphi2];

J = [ Rot1*[beta1*s1;-gamma1*c1] Rot2*[-beta2*s2;gamma2*c2]];

rho

  function f = rho(z,plist)
% f = rho(z,plist)
% z is a 2-vector whose components are the orbit parameters A and P
% respectively. plist is a 2-column matrix. plist(i,1) is the i-th
% observed radius vector length and plist(i,2) is the cosine of the
% corresponding observed polar angle.

A = z(1);
P = z(2);
r = plist(:,1);
c = plist(:,2);
f = r - (2*A*P)./(P*(1-c) + A*(1+c));

Jrho

  function J = Jrho(z,plist)
% J = Jrho(z,plist)  
% J is the Jacobian of the function rho at z.

A = z(1);
P = z(2);
c = plist(:,2);
denom = (P*(1-c) + A*(1+c)).^2;
J = -[ (2*P^2)*(1-c)./denom   (2*A^2)*(1+c)./denom];