% This code simulates the dynamic response of parallel plate capacitor on
% the application of voltage across its plates, where one plate is attached
% to a spring and other is fixed.

clear all;
clc;

m = 1.1615e-10;     % mass of the moving plate or approximation thereof
k = 164;            % spring constant or approximation thereof
c = 0.0000;         % damping coefficient
g0 = 1.6e-6;        % initial capacitor gap
eps0 = 8.854e-12;   % permittivity of the capacitor medium
area = 25e-9;       % area of the capacitor plate
vDC = 27;
vAC = 0.0;          % amplitude of AC volatage
freq_fac = 0.1;     % freq_fac*nat_freq = applied AC voltage frequency 
omega= sqrt(k/m)*freq_fac; % frequency of applied AC signal rad/sec

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

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

t=0:10^aa:10^(aa+2.2);
y0=[0;0];           % Initial conditions [position; velocity]

[T, Y] = ode45(@(t,y)myode(t,y,m,k,c,g0,vDC,vAC,eps0,area,omega),t,y0); % Solve ODE

tmin= min(T); 
tmax= max(T);
ymin= min(Y(:,1));
ymax= max(Y(:,1));

yall= roots([k, -2*k*g0, k*g0^2, -vDC^2*eps0*area/2]);
ymean=yall(3);

Fnet = vDC^2*eps0*area./(2*(g0-Y(:,1)).^2) - c*Y(:,2) - k*Y(:,1);
KE = 0.5*m*Y(:,2).^2;

scrsz = get(0,'ScreenSize');
figure('Position',scrsz);
ti=size(Y,1);
x=[0, 2*g0]; 
h(1)=subplot(1,2,1);
h(2)=subplot(1,2,2);

for i=1:ti
    subplot(h(1));
    cla (h(1));
    hold on;
    plot(x,[-g0,-g0],'-k','linewidth',5);
    plot(x,[-Y(i,1), -Y(i,1)],'-r','linewidth',5);
    plot(x,[0,0],'--b', x,[-ymean,-ymean],'--k');
    axis([-0.2*g0, 2.2*g0, -1.1*g0, 1.2*g0]);
    
    text(0,0.6*g0,['Net force = ' num2str(Fnet(i)) ' N']);
    text(0,0.5*g0,['KE =        ' num2str(KE(i)) ' J']);
    hold off
    
    subplot(h(2));
    cla (h(2));
    hold on;
    plot(T(1:i),Y(1:i,1),'-',T(i),Y(i,1),'rd');
    plot([tmin, tmax],[ymean,ymean],'--k');
    title('Time history');
    xlabel('Time [sec]');
    ylabel('y [m]');
    axis([tmin tmax ymin ymax]);
    pause(0.05)
    hold off
end


y1=(-g0/10):g0/1000:0.7*g0;
PE0= 0.5*k*y1.^2;
PE= 0.5*k*y1.^2 - vDC^2*eps0*area./(2*(g0-y1));
PE1= 0.5*k*Y(:,1).^2 - vDC^2*eps0*area./(2*(g0-Y(:,1)));

figure(2)

hold on
plot(y1,PE0,'-k');
plot(0,0,'o','markersize',10,'markerface','b')
xlabel('position Y (m)');
ylabel('Potential Energy (J)');
pause

plot(y1,PE,'-r','linewidth',2);
plot(0,-vDC^2*eps0*area/2/g0,'o','markersize',10,'markerface','b');
xlabel('position Y (m)');
ylabel('Potential Energy (J)');
axis([min(y1), max(y1), min(PE), 0.5*max(PE0)])
pause

for i=1:ti
    clf
    hold on
    plot(y1,PE,'-r','linewidth',2);
    plot(Y(i,1),PE1(i),'o','markersize',10,'markerface','b')
    xlabel('position Y (m)');
    ylabel('Potential Energy (J)');
    pause(0.05)
end
hold off