%
% Direct numerical finite central difference integration scheme to solve
% a second order ordinary differential equation encountered in 1 DOF
% dynamic system. Forcing function is dependent on time as well as the
% position in a nonlinear fashion. Source of force: electrostatic force
% in a parallel plate capacitor. -- G. K. Ananthasuresh
%
clear
hold off
clc
clg
echo off

% Data
%
% m = mass of the moving plate or approximation thereof
% k = spring constant or approximation thereof
% c = damping coefficient
% vdc = DC voltage
% vac = AC voltage
% g0 = initial capacitor gap
% eps = permittivity of the capacitor medium
% area = area of the capacitor plate
% N = number of time steps in integration
%    (N*0.1*sqrt(m/k) will be the final time)
% freq_fac*nat_freq = applied AC voltage frequency 
%
% Change also the initial conditions x(1) and v(1) if necessary.
%

%mi = input('Enter mi: ');
m = 1.1615e-10;
k = 164.0;
%c = input('Enter c: ');
c = 0;

% VdynPI = 27.545471
vdc = 2.0;
vac = 20.0;
freq_fac = input('Enter freq_fac: ');
%freq_fac = 0.1;
freq = sqrt(k/m)/2/pi*freq_fac;         % freq is in Hertz
zeta = c/(2*m*sqrt(k/m));
freqd = freq*sqrt(1-zeta^2);    % freqd is in Hertz

g0 = 1.6e-6;
eps = 8.854e-12;
area = 25e-9;

Vpi = sqrt( 8*k*g0^3 / (27*eps*area) );

% Derived data
epsa = eps*area;
omega = freq*2*pi;      % Omega is in rad/sec

% Time step = one tenth of the time period
h = 0.1*sqrt(m/k);
N = 500;

%
% Initial conditions
%
t(1) = 0;
x(1) = 0;
v(1) = 0;

%
% Derived conditions at t = h sec. (because the method is not self starting...)
%
t(2) = t(1) + h;
ve = vdc + vac*sin(omega*t(1));
f = (1/m) * ( ve^2*epsa/2/(g0-x(1))^2 - c*v(1) - k*x(1) );
x(2) = x(1) + h*v(1) + (h^2/2)*f;
ve = vdc + vac*sin(omega*t(2));
ff = ve^2*epsa/2/(g0-x(2))^2;
v(2) = ( x(2) - x(1) + (h^2/2/m)*(ff-k*x(2)) ) / (h+h^2*c/2/m);

% Source AC voltage
vs(1) = vac*sin(omega*t(1));
vs(2) = vac*sin(omega*t(2));

%
% Start the iteration with i = 3
%
for i = 3:N,
 t(i) = t(i-1) + h;
 ve = vdc + vac*sin(omega*t(i-1));
 vs(i) = ve - vdc;
 f = (1/m) * ( ve^2*epsa/2/(g0-x(i-1))^2 - c*v(i-1) - k*x(i-1) );
 x(i) = 2*x(i-1) - x(i-2) + h^2*f;
 ff = ve^2*epsa/2/(g0-x(i))^2;
 v(i) = ( x(i) - x(i-1) + (h^2/2/m)*(ff-k*x(i)) ) / (h+h^2*c/2/m);
end

for i = 1:N,
Force(i) = (vdc+vac*sin(omega*t(i)) )^2 *epsa / (g0^2*2) / (1-x(i)/g0)^2;
 xf = x(i)/g0;
af(i) = (vdc+vac*sin(omega*t(i)) )^2 *epsa / (g0^2*2) * (1+2*xf+3*xf^2);
end

echo on

%
% Plot results
%
plot(t,x,'-');
title('Time history');
xlabel('Time [sec]');
ylabel('x [m]');

pause
clg

plot(x,v);
title('Phase Plane Portrait');
xlabel('x (position)');
ylabel('v (Velocity)');
grid

pause
clg
subplot(2,1,1);
plot(t,x,'-');
title('Time history');
xlabel('Time [sec]');
ylabel('x [m]');
subplot(2,1,2);
plot(x,v);
title('Phase Plane Portrait');
xlabel('x (position)');
ylabel('v (Velocity)');