%% Nonlinear Solutions Example - Bracketing
% Interval Halving, False Position

toler=0.0001;
alpha=40;

% Trial and Error
a=30;           % Lower Guess
da=0.01;        % Increment    
b=35;           % Upper Guess
k=0;
phi=a:da:b;    % Plot range

while (abs(f(alpha,a))>toler) & (a<b);
    k=k+1;
    a=a+da;
    % Plot
    figure(1); clf;
    dphi=(b-a)/200;
    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--');    
    plot([b b], get(gca, 'ylim'),'r--');   
    text(a+10*dphi,f(alpha,a)+0.03,['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(0.001);
end;

% Interval Halving
a=30;           % Lower Guess
b=40;           % Upper Guess
c=(a+b)/2;      % Midpoint
k=0;
while (abs(f(alpha,c))>toler);
    k=k+1;
    if (f(alpha,a)*f(alpha,c)<0);
        a=a; 
        b=c; 
    else
        a=c; 
        b=b; 
    end;
    c=(a+b)/2;
    % Plot
    figure(2); clf;
    dphi=(b-a)/200;
    phi=a-dphi*10:dphi:b+dphi*10;    % Plot range
    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--');    
    plot([b b], get(gca, 'ylim'),'r--');   
    plot([c c], get(gca, 'ylim'),'k--');   
    text(c+10*dphi,f(alpha,c)+0.03,['c=',num2str(c),char(10),'f(c)=',num2str(f(alpha,c)),char(10),'k=',num2str(k)],'FontSize',14);
    xlabel('\phi (deg.)');
    ylabel('f(\phi)');
    pause;
end;

% False Position (linear interpolation)
a=20;           % Lower Guess
b=30;           % Upper Guess
c=0;
k=0;
while (abs(f(alpha,c))>toler);
    k=k+1;
    g=(f(alpha,b)-f(alpha,a))/(b-a);
    c=b-f(alpha,b)/g;
    %c=b-(f(alpha,b)*(b-a))/(f(alpha,b)-f(alpha,a));
    if (f(alpha,a)*f(alpha,c)<0);
        a=a; 
        b=c; 
    else
        a=c; 
        b=b; 
    end;
    % Plot
    figure(3); clf;
    dphi=(b-a)/200;
    phi=a-dphi*10:dphi:b+dphi*10;    % Plot range
    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--');    
    plot([b b], get(gca, 'ylim'),'r--');   
    plot([c c], get(gca, 'ylim'),'k--');   
    text(c+10*dphi,f(alpha,c)+0.03,['c=',num2str(c),char(10),'f(c)=',num2str(f(alpha,c)),char(10),'k=',num2str(k)],'FontSize',14);
    xlabel('\phi (deg.)');
    ylabel('f(\phi)');
    pause;
end;