clear all
% define parameters
g= 9.8;
m = 750;
mu = 4*pi*1e-7;
N = 324;
A = 0.021;

h0=0.008;
i0 = h0*sqrt(m*g/(mu*N^2*A*0.25));

%define state space model
A = [0  1; 2*g/h0 0];
B = [0   0; -2*g/i0   1/m];
C = eye(2);
D = zeros(2);
maglev_sys = ss(A,B,C,D);

% sepcify inputs and intial conditions
Ts = 1e-3;
time = [0:Ts:10]';
U  = zeros(length(time),2);
X0 = [0.001*1e-3; 0];

%% simulate 
Y = lsim(maglev_sys,U,time,X0);


% plot
figure;
subplot(211);
plot(time,Y(:,1));
set(gca,'fontsize',18);
xlabel('time (sec)');
ylabel(' position (rad)');
subplot(212);
plot(time,Y(:,2));
set(gca,'fontsize',18);
xlabel('time (sec)');
ylabel(' velocity (rad/s)');

% the sysyem is unstable, so the states grow to infinity

% you'll notice it looks like the states are close to until the very end
% and then they shoot up to very large values. That is a misleading
% impression created by the linear y-axis scale

% plot the states in log scale to see the continous growth over time
figure;
subplot(211);
semilogy(time,Y(:,1));
set(gca,'fontsize',18);
xlabel('time (sec)');
ylabel(' position (rad)');
subplot(212);
semilogy(time,Y(:,2));
set(gca,'fontsize',18);
xlabel('time (sec)');
ylabel(' velocity (rad/s)');