%% Nonlinear Solutions Example - Open Domain
% Newton's Method and Secant Method

toler=0.000001;
alpha=40;
maxiter=100;

% Newton's Method Iteration
a=25;           % Guess
xl=20;          % Upper Bound for Plotting
xh=40;          % Lower Bound for Plotting
k=0;            % Zeroth Iteration
dphi=(xh-xl)/200;
phi=xl:dphi:xh;     % Plot range

while (abs(f(alpha,a))>toler) & (k<maxiter);
    k=k+1;
    a=a-(360/(2*pi))*f(alpha,a)/fprime(alpha,a);
    % Plot
    figure(1); clf;
    plot(phi,f(alpha,phi),'k-');
    axis([min(phi) max(phi) -0.2 0.2]);
    hold on;
    plot(get(gca, 'xlim'),[0 0],'k:');   
    plot([a a], get(gca, 'ylim'),'b--');    
    text(a+10*dphi,f(alpha,a)-0.05*dphi,['a=',num2str(a),char(10),'f(a)=',num2str(f(alpha,a)),char(10),'k=',num2str(k)],'FontSize',14);
    xlabel('\phi (deg.)');
    ylabel('f(\phi)');
    pause;
end;

% Newton's Method Iteration, Numerical Derivative
a=50;           % Guess
xl=25;          % Upper Bound for Plotting
xh=40;          % Lower Bound for Plotting
k=0;            % Zeroth Iteration
dphi=(xh-xl)/200;
phi=xl:dphi:xh;     % Plot range

while (abs(f(alpha,a))>toler) & (k<maxiter);
    k=k+1;
    da=toler;
    derivative=(f(alpha,a+da)-f(alpha,a))/da;
    a=a-f(alpha,a)/derivative;
    % Plot
    figure(1); clf;
    plot(phi,f(alpha,phi),'k-');
    axis([min(phi) max(phi) -0.2 0.2]);
    hold on;
    plot(get(gca, 'xlim'),[0 0],'k:');   
    plot([a a], get(gca, 'ylim'),'b--');    
    text(a+10*dphi,f(alpha,a)-0.05*dphi,['a=',num2str(a),char(10),'f(a)=',num2str(f(alpha,a)),char(10),'k=',num2str(k)],'FontSize',14);
    xlabel('\phi (deg.)');
    ylabel('f(\phi)');
    pause;
end;

% Secant Method Iteration
a0=30;          % Guess 0
ak=20;          % Guess 1
xl=25;          % Upper Bound for Plotting
xh=40;          % Lower Bound for Plotting
k=0;            % Zeroth Iteration
dphi=(xh-xl)/200;
phi=xl:dphi:xh;     % Plot range

while (abs(f(alpha,ak))>toler) & (k<maxiter);
    k=k+1;
    slope=(f(alpha,ak)-f(alpha,a0))/(ak-a0);
    a0=ak;
    ak=ak-f(alpha,ak)/slope;
    
    % Plot
    figure(3); clf;
    plot(phi,f(alpha,phi),'k-');
    axis([min(phi) max(phi) -0.2 0.2]);
    hold on;
    plot(get(gca, 'xlim'),[0 0],'k:');   
    plot([ak ak], get(gca, 'ylim'),'b--');    
    text(ak+10*dphi,f(alpha,ak)-0.05*dphi,['ak=',num2str(ak),char(10),'f(a)=',num2str(f(alpha,ak)),char(10),'k=',num2str(k)],'FontSize',14);
    xlabel('\phi (deg.)');
    ylabel('f(\phi)');
    pause;
end;
% 
