%This program solves the differential equation:
%Dy''''+Ty''=-eV^2/2y^2+P 
%for a clamped circular plate (or diaphragm) over a ground plane under voltage 
%and/or pressure load using the finite difference method.
%It also solves for the capacitance of the plate to ground.
%Note that the radius of both the top plate and Gnd Plate are both user
%specified.
%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
%4/1/94

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

if case==1

clear;
%d=input('Enter discretization number (50 to 100 range is suggested): ');
d=128;
N=d+3;
%Z=input('Enter newton iteration number (5 to 10 range is suggested): ');
Z=51;
%R=1e-6*input('Enter movable conductor radius in um: ');
R=100e-6;
%RF=1e-6*input('Enter fixed ground conductor radius in um: ');
%RF=250e-6;
%if RF >= R
RF=R;
%end
%t=1e-6*input('Enter movable conductor thickness in um: ');
t=1e-6;
%g=1e-6*input('Enter gap in um: ');
g=1e-6;
%E=1e9*input('Enter movable conductor youngs modulus in GPa: ');
E=100e9;
%nu=input('Enter movable conductor poisson ratio: ');
nu=0;
s=1e6*input('Enter movable conductor residual stress in MPa: ');
%s=43e6;
%P=input('Enter applied pressure in Pa:');
P=0;
h=R/(d-1);
D=E*t^3/(12*(1-nu^2));
T=s*t;
a=T/D;

end

e=8.8541878e-12;

g1 = 1.55;
g2 = 1.65;
g3 = 0.00;
S3 = s*t*g^3;
B3 = E/(1-nu^2)*t^3*g^3;
u3 = (12*S3/B3)^(1/2) * g2 * R/2;
D3 = 1 + (2*(1-cosh(u3))/u3/sinh(u3));
V_pull_in_HH=real(g1*S3/(e*R^2*D3))^(1/2)
g1a = 1.55;
g2a = 1.63978;
u3a = (12*S3/B3)^(1/2) * g2a * R/2;
D3a = 1 + (2*(1-cosh(u3a))/u3a/sinh(u3a));
V_pull_in_HH2=real(g1a*S3/(e*R^2*D3a))^(1/2)
Bckl_stress_in_MPa = -14.682*E*t^2/12/R^2/1e6

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

if mode==1

b=P/D;
A=zeros(N,N);
C=zeros(N,1);

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

for i=3:N-2
r=(i-2)*h;
A(i,i-2)=1-h/r;
A(i,i-1)=-4+2*h/r-h^2*(1/r^2+a)-.5*h^3*(1/r^3-a/r);
A(i,i)=6+2*h^2*(1/r^2+a);
A(i,i+1)=-4-2*h/r-h^2*(1/r^2+a)+.5*h^3*(1/r^3-a/r);
A(i,i+2)=1+h/r;
C(i)=b*h^4;
end

ig=A\C;

V=input('Enter applied voltage: ');

b=-e*V^2/(2*D);
c=P/D;

y=g*ones(N,1);
y=y+ig;
A=zeros(N,N);
C=zeros(N,1);
cap=zeros(N,1);
y1p=zeros(N,1);
y2p=zeros(N,1);
y3p=zeros(N,1);
y4p=zeros(N,1);

PN=round((RF/R)*d+1);

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

for j=1:Z

for i=3:PN
y1p(i)=(-y(i-1)+y(i+1))/(2*h);
y2p(i)=(y(i-1)-2*y(i)+y(i+1))/h^2;
y3p(i)=(-y(i-2)+2*y(i-1)-2*y(i+1)+y(i+2))/(2*h^3);
y4p(i)=(y(i-2)-4*y(i-1)+6*y(i)-4*y(i+1)+y(i+2))/h^4;
r=(i-2)*h;
A(i,i-2)=1-h/r;
A(i,i-1)=-4+2*h/r-h^2*(1/r^2+a)-h^3*(1/r^3-a/r)/2;
q=2*h^4/y(i)*(y4p(i)+2*y3p(i)/r-y2p(i)*(1/r^2+a)+y1p(i)*(1/r^3-a/r)-c);
A(i,i)=6+2*h^2*(1/r^2+a)+q;
A(i,i+1)=-4-2*h/r-h^2*(1/r^2+a)+h^3*(1/r^3-a/r)/2;
A(i,i+2)=1+h/r;
C(i)=h^4*(b/y(i)^2-y4p(i)-2*y3p(i)/r+y2p(i)*(1/r^2+a)-y1p(i)*(1/r^3-a/r)+c);
end

for i=PN+1:N-2
y1p(i)=(-y(i-1)+y(i+1))/(2*h);
y2p(i)=(y(i-1)-2*y(i)+y(i+1))/h^2;
y3p(i)=(-y(i-2)+2*y(i-1)-2*y(i+1)+y(i+2))/(2*h^3);
y4p(i)=(y(i-2)-4*y(i-1)+6*y(i)-4*y(i+1)+y(i+2))/h^4;
r=(i-2)*h;
A(i,i-2)=1-h/r;
A(i,i-1)=-4+2*h/r-h^2*(1/r^2+a)-h^3*(1/r^3-a/r)/2;
q=2*h^4/y(i)*(y4p(i)+2*y3p(i)/r-y2p(i)*(1/r^2+a)+y1p(i)*(1/r^3-a/r)-c);
A(i,i)=6+2*h^2*(1/r^2+a)+q;
A(i,i+1)=-4-2*h/r-h^2*(1/r^2+a)+h^3*(1/r^3-a/r)/2;
A(i,i+2)=1+h/r;
C(i)=h^4*(-y4p(i)-2*y3p(i)/r+y2p(i)*(1/r^2+a)-y1p(i)*(1/r^3-a/r)+c);
end

f=A\C;
y=y+f;
end

cap(3)=e*pi*h^2/y(2);
for i=4:PN
r=(i-2)*h;
cap(i)=cap(i-1)+e*pi*h^2*(2*i-5)/(.5*(y(i)+y(i-1)));
end

y=1e6*y;
cap=1e12*cap;
f=1e6*f;

plot(y);
axis([2 d+1 0 g*2e6]);
xlabel('Discretization in x');
ylabel('Gap (um)');
Center_gap=y(2)
%Center_gap_estimate=1e6*(g+P/k)
Error=f(2)
Cap_pF=cap(PN)
k1=16*E*t^3/(3*R^4*(1-nu^2));
k2=4*T/R^2;
k=k1+k2;
V_pull_in_estimate=((8*k*(g+P/k)^3)/(27*e))^.5

else

b=P/D;
A=zeros(N,N);
C=zeros(N,1);

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

for i=3:N-2
r=(i-2)*h;
A(i,i-2)=1-h/r;
A(i,i-1)=-4+2*h/r-h^2*(1/r^2+a)-.5*h^3*(1/r^3-a/r);
A(i,i)=6+2*h^2*(1/r^2+a);
A(i,i+1)=-4-2*h/r-h^2*(1/r^2+a)+.5*h^3*(1/r^3-a/r);
A(i,i+2)=1+h/r;
C(i)=b*h^4;
end

ig=A\C;

e=8.854e-12;
k1=16*E*t^3/(3*R^4*(1-nu^2));
k2=4*T/R^2;
k=k1+k2;
VMIN=((8*k*(g+P/k)^3)/(27*e))^.5;
VMAX=2*VMIN;
V=0;
c=P/D;

while abs(VMAX-VMIN) > 0.0001

b=-e*V^2/(2*D);
y=g*ones(N,1);
y=y+ig;
A=zeros(N,N);
C=zeros(N,1);
cap=zeros(N,1);
y1p=zeros(N,1);
y2p=zeros(N,1);
y3p=zeros(N,1);
y4p=zeros(N,1);

PN=round((RF/R)*d+1);

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

for j=1:Z

for i=3:PN
y1p(i)=(-y(i-1)+y(i+1))/(2*h);
y2p(i)=(y(i-1)-2*y(i)+y(i+1))/h^2;
y3p(i)=(-y(i-2)+2*y(i-1)-2*y(i+1)+y(i+2))/(2*h^3);
y4p(i)=(y(i-2)-4*y(i-1)+6*y(i)-4*y(i+1)+y(i+2))/h^4;
r=(i-2)*h;
A(i,i-2)=1-h/r;
A(i,i-1)=-4+2*h/r-h^2*(1/r^2+a)-h^3*(1/r^3-a/r)/2;
q=2*h^4/y(i)*(y4p(i)+2*y3p(i)/r-y2p(i)*(1/r^2+a)+y1p(i)*(1/r^3-a/r)-c);
A(i,i)=6+2*h^2*(1/r^2+a)+q;
A(i,i+1)=-4-2*h/r-h^2*(1/r^2+a)+h^3*(1/r^3-a/r)/2;
A(i,i+2)=1+h/r;
C(i)=h^4*(b/y(i)^2-y4p(i)-2*y3p(i)/r+y2p(i)*(1/r^2+a)-y1p(i)*(1/r^3-a/r)+c);
end

for i=PN+1:N-2
y1p(i)=(-y(i-1)+y(i+1))/(2*h);
y2p(i)=(y(i-1)-2*y(i)+y(i+1))/h^2;
y3p(i)=(-y(i-2)+2*y(i-1)-2*y(i+1)+y(i+2))/(2*h^3);
y4p(i)=(y(i-2)-4*y(i-1)+6*y(i)-4*y(i+1)+y(i+2))/h^4;
r=(i-2)*h;
A(i,i-2)=1-h/r;
A(i,i-1)=-4+2*h/r-h^2*(1/r^2+a)-h^3*(1/r^3-a/r)/2;
q=2*h^4/y(i)*(y4p(i)+2*y3p(i)/r-y2p(i)*(1/r^2+a)+y1p(i)*(1/r^3-a/r)-c);
A(i,i)=6+2*h^2*(1/r^2+a)+q;
A(i,i+1)=-4-2*h/r-h^2*(1/r^2+a)+h^3*(1/r^3-a/r)/2;
A(i,i+2)=1+h/r;
C(i)=h^4*(-y4p(i)-2*y3p(i)/r+y2p(i)*(1/r^2+a)-y1p(i)*(1/r^3-a/r)+c);
end

f=A\C;
y=y+f/5;
end

cap(3)=e*pi*h^2/y(2);
for i=4:PN
r=(i-2)*h;
cap(i)=cap(i-1)+e*pi*h^2*(2*i-5)/(.5*(y(i)+y(i-1)));
end

y=1e6*y;
cap=1e12*cap;
f=1e6*f;

Center_gap=y(2)
%Center_gap_estimate=1e6*(g+P/k)
Error=f(2)
V_max=VMAX
V_min=VMIN
V_pull_in=V
k1=16*E*t^3/(3*R^4*(1-nu^2));
k2=4*T/R^2;
k=k1+k2;
V_pull_in_estimate=((8*k*(g+P/k)^3)/(27*e))^.5

if abs(f(2)) < 5e-4
VMIN=V;
V=(VMAX+V)/2;
end

if abs(f(2)) >= 5e-4
VMAX=V;
V=(VMIN+V)/2;
end

end

plot(y);
axis([2 d+1 0 g*2e6]);
xlabel('Discretization in r');
ylabel('Gap (um)');

end