% This code was written by Mayur Birla for Mechanics of Microsystem (ME237)
% course assignment 4, Problem 2
clear all;
clc;

E=155e9;    % Young's modulus
w=20e-6;    % width of cantilever beam
t=5e-6;     % thickness of cantilever beam
l=200e-6;   % length of cantilever beam
I=w*t^3/12; % Moment of inertia about axis perpendicular to t (z-aixs)

k=3*E*I/l^3;% stiffness of cantilever

q0=1e-15;   % charge on substrate which just below ther cantilever tip

g0=2e-6;    % gap between q0 and cantilever tip
e0=8.854e-12; % permissivity of free space

c2= -2*g0;  % co-efficient of x^2
c3= g0^2;   % co-efficient of x

q=0:1e-14:2e-12; % charge at cantilever tip

a=size(q,2);
x0=zeros(3,a);
x=zeros(3,a);
y=zeros(3,a);

col=[1 0 0; 0.5 0.5 0; 0 0 1]; % red, green, blue

qcr_th= 16*pi*e0*k*g0^3/(27*q0); % theoritcal Q critical

for i=1:a
    c4= -q(i)*q0/(4*pi*e0*k);
    x0(:,i)= sort(roots([1  c2 c3 c4]));
end

x=real(x0);
y=imag(x0);

xmin=min(min(x));
xmax=max(max(x));
ymin=min(min(y));
ymax=max(max(y));

for i=1:a
    clf
    hold on
    plot(x(1,1:i),y(1,1:i),'-','color',col(1,:),'linewidth',1.5);
    plot(x(1,i),y(1,i),'d','color',col(1,:),'markersize',5);
    
    plot(x(2,1:i),y(2,1:i),'-','color',col(2,:),'linewidth',1.5);
    plot(x(2,i),y(2,i),'d','color',col(2,:),'markersize',5);
    
    plot(x(3,1:i),y(3,1:i),'-','color',col(3,:),'linewidth',1.5);
    plot(x(3,i),y(3,i),'d','color',col(3,:),'markersize',5);
    text(xmax*0.7,ymax*0.7,['Q = ', num2str(q(i)), ' C'])
    
    xlabel('Real value of x (\mum)');
    ylabel('Imaginary value of x (\mum)');
    title('Root locus plot of x wrt Q');
    
    axis([xmin xmax ymin ymax]);
    grid on
    pause(0.001)
    hold off
    
    if (max(i==[1,2,160,161]))
        pause
    end
end


