% This program solves the beam differential equation: 
%
%      EIy'''' = -eV^2w/2y^2 (1 + 0.65y/w)
%
% for a cantilever beam over a ground plane using the finite difference
% method.  Note that the fringing term is included in the above equation.
% There are two Modes to this program.  The Normal Mode is the familiar
% single voltage/single deflection ouput mode.  The new Pull-In Mode
% will automatically iteratively calculate the exact Pull-In Voltage for
% the input'ed case to within 2 decimal points precison without the user
% having to manually iterate in the Normal Mode as was the case in the past.
% Peter Osterberg, MIT
% 4/1/94
%
% This script is now updated to model stepped up boundaries found in
% fabrication processes with conformal deposition, like the MUMPs 
% process of MCNC.  The step is modeled by a vertical height, specified
% by the user, and has the same width and thickness as the cantilever.  
% By fixing  the step height to zero an ideal fixed end can still be 
% modeled.
%
% Also included in this script is the ability to model a curled cantilever
% subjected to a stress gradient before release.  By convention, a
% positive stress gradient deforms a beam upward, away from the ground
% electrode.  The model accounts for this deformation by varying the
% gap as a function of the axial dimension, and keeping the initial
% profile of the cantilever flat.  Applying a moment to the end of
% the cantilever to give the beam an initial curvature also works fine
% except for voltages near pull-in.
%
% Raj Gupta, MIT
% 1st October 1995

trial=input('New case ? (1=yes or 0=no): ');

if trial==1

clear;
d = input('Enter discretization number (50 to 100 range is suggested): ');
%d = 50;
N = d+4;
Z = input('Enter newton iteration number (5 to 10 range is suggested): ');
%Z = 5;
L = 1e-6*input('Enter length in um: ');
%L = 100e-6;
w = 1e-6*input('Enter width in um: ');
%w = 20e-6;
t = 1e-6*input('Enter thickness in um: ');
%t = 2e-6;
g = 1e-6*input('Enter gap in um: ');
%g = 1e-6;
E = 1e9*input('Enter youngs modulus in GPa: ');
%E = 160e9;
step = 1e-6*input('Enter step height in um: ');
%step = 0;
dsdy = 1e12*input('Enter transverse stress gradient (MPa/um): ');
%dsdy = 0;
end

tip = d+2;
e = 8.8541878e-12;
h = L/(d-1);
if dsdy~=0
    % Following assumption is only true for a beam with a
    % rectangular cross-section
    Radius = E/dsdy;
    mxphi = g/L + L/(2*Radius);
else
    mxphi = g/L;
end
%  Two widths allow for trapezoidal geometry,
%  w1 = top width, w = bottom width
w1 = w;
I = t^3/36 * (w^2+4*w*w1+w1^2)/(w+w1);
k = 2*E*t^3/(3*L^4);
V_pull_in_estimate = (8*k*g^3/(27*e))^0.5
g1 = 0.28;
c1 = 0.42;
%  Following needs to be changed to reflect
%  trapezoidal cross sectional geometry:
V_pull_in_HH = real((g1*E*t^3*g^3/(e*L^4*(1+c1*g/w)))^0.5)

y = g*ones(N,1);
y_old = y;
phi_old = 0;
r = zeros(N,1);
A = zeros(N,N);
C = zeros(N,1);
y4p = zeros(N,1);
con = zeros(N,1);

if dsdy ~= 0
  for i = 1:N
    con(i) = Radius - sign(Radius)*(Radius^2-((i-3)*h)^2)^0.5;
  end
else
    con = zeros(N,1);
end

A(1,3) =  1;
A(2,2) = -1;
A(2,4) =  1;
A(N-1,N-3) =  1;
A(N-1,N-2) = -2;
A(N-1,N-1) =  1;
A(N,N-4) = -1;
A(N,N-3) =  2;
A(N,N-1) = -2;
A(N,N) =  1;

mode=input('Normal Mode or Pull-In Mode ? (1=Normal or 0=Pull-In): ');

if mode == 1

% BEGIN MODE 1 = NORMAL

V = input('Enter voltage: ');
b = -e*V^2*w/(2*E*I);
c = 0.65*b/w;
phi = 0;
count = 0;
r(tip) = 1;
while abs(r(tip)/y(tip)) > 1e-5
    for i=3:N-2
	y4p(i) = (y(i-2)-4*y(i-1)+6*y(i)-4*y(i+1)+y(i+2))/h^4;
	A(i,i-2) =  1;
	A(i,i-1) = -4;
	A(i,i) = 6+2*h^4/(y(i)+con(i))*(y4p(i)-c/(2*(y(i)+con(i))));
	A(i,i+1) = -4;
	A(i,i+2) =  1;
	C(i) = h^4*(b/(y(i)+con(i))^2-y4p(i)+c/(y(i)+con(i)));
    end
    r=A\C;
    if count > Z
      y=y+(r/(count)^0.5);
    else
      y=y+r;
    end
    y(N-2)
    count = count + 1;
    if step ~= 0
        Mo = (y(4)-2*y(3)+y(2))/(h^2)*E*I;
	phi = -Mo*(step)/(E*I);
	y1p3 = -(y(4) -  y(2))/(2*h);
	C(2) = -2*h*(phi - y1p3);
    end
end

% END OF NEWTON ITERATIONS FOR MODE 1

y=(y+con);
cap = 0;
for i=4:N-2
  f = w/y(i);
  cap = cap+(e*w*h/y(i))*4*((4/(pi*f))^(3/20)+pi/4*f+(1+4/(pi*f))^(-1/2)*atanh((1+4/(pi*f))^(-1/2)))/(pi*f);
end
Capacitance = cap

plot(1e6*y);
axis([3 tip 0 (2*g)*1e6]);
xlabel('Discretization in x');
ylabel('Gap (um)');
Tip_gap = y(tip)
Error = r(tip)/y(tip)

% END OF MODE 1 = NORMAL

else

% BEGIN MODE 2 = PULLIN

VMIN = 0;
VMAX = 500*V_pull_in_HH;
V = 0.01*VMAX;
phi = 0;

% BEGIN BINARY SEARCH FOR PULLIN VOLTAGE

while abs(VMAX-VMIN) > 0.001

b = -e*V^2*w/(2*E*I);
c = 0.65*b/w;

% BEGIN NEWTON ITERATIONS FOR VOLTAGE = V

r(tip) = 0.5;
err = 1;
count = 1;
lim = 2*floor((log(VMAX/(VMAX-VMIN)))) + Z;
while (((abs(r(tip)) < err & count <= lim)|(count < Z))&(phi>=0)&(phi<=mxphi))
    if (count < 5)
        err = abs(r(tip));
    end
    for i=3:N-2
	y4p(i) = (y(i-2)-4*y(i-1)+6*y(i)-4*y(i+1)+y(i+2))/h^4;
	A(i,i-2) =  1;
	A(i,i-1) = -4;
	A(i,i) = 6+2*h^4/(y(i)+con(i))*(y4p(i)-c/(2*(y(i)+con(i))));
	A(i,i+1) = -4;
	A(i,i+2) =  1;
	C(i) =  h^4*(b/(y(i)+con(i))^2-y4p(i)+c/(y(i)+con(i)));
    end
    r=A\C;
    if count > 10
      y=y+(r/(count)^0.5);
    else
      y=y+r;
    end
    count = count + 1;
    if step ~= 0
        Mo = (y(4)-2*y(3)+y(2))/(h^2)*E*I;
        phi = -Mo*(step)/(E*I);
	y1p3 = -(y(4) -  y(2))/(2*h);
	C(2) = -2*h*(phi - y1p3);
    end
end

% END NEWTON ITERATIONS FOR VOLTAGE = V

V_pull_in = V
Tip_gap = y(tip)+con(tip)
if step ~= 0 
  rotation = phi
end
Error = r(tip)/y(tip)
iterations = count-1

if ((abs(Error)>1e-4)|(y(tip)<-con(tip))|(y(tip)>g)|(phi<0)|(phi>mxphi))
  VMAX = V;
  V = (VMIN+V)/2;
  y = y_old;
  phi = phi_old;
  C(2) = 0;
else
  VMIN = V;
  V = (VMAX+V)/2;
  y_old = y;
  phi_old = phi;
end

end

% END BINARY SEARCH FOR PULLIN VOLTAGE

% y=(y+con);
cap = 0;
for i=4:N-2
  f = w/y(i);
  cap = cap+(e*w*h/y(i))*4*((4/(pi*f))^(3/20)+pi/4*f+(1+4/(pi*f))^(-1/2)*atanh((1+4/(pi*f))^(-1/2)))/(pi*f);
end
Capacitance = cap

plot(y*1e6);
axis([3 tip 0 (2*g)*1e6]);
xlabel('Discretization in x');
ylabel('Gap (um)');

end