function MotorSim(ControlChange)

% Two-phase step motor simulation.
% This file can be used to investigate the effects of linearization.
% If ControlChange = 0 then the nonlinear and linearized simulations should
% match exactly (in steady state). As ControlChange increases, the linearized simulation
% diverges more and more from the nonlinear simulation.

if ~exist('ControlChange', 'var')
    ControlChange = 0;
end

Ra = 1.9; % Winding resistance
L = 0.003; % Winding inductance
lambda = 0.1; % Motor constant
J = 0.00018; % Moment of inertia
B = 0.001; % Coefficient of viscous friction

dt = 0.0005; % Integration step size
tf = 1; % Simulation length 

x = [0; 0; 0; 0]; % Initial state
xlin = x; % Linearized approximation of state
w = 2 * pi; % Control input frequency

dtPlot = 0.005; % How often to plot results
tPlot = -inf;

% Initialize arrays for plotting at the end of the program
xArray = [];
xlinArray = [];
tArray = [];

dx = x - xlin; % Difference between true state and linearized state

tStart = cputime;

% Begin simulation loop
for t = 0 : dt : tf
    if t >= tPlot + dtPlot
        % Save data for plotting
        tPlot = t + dtPlot - eps;
        xArray = [xArray x];
        xlinArray = [xlinArray xlin];
        tArray = [tArray t];
    end
    % Nonlinear simulation
    ua0 = sin(w*t); % nominal winding A control input
    ub0 = cos(w*t); % nominal winding B control input
    ua = (1 + ControlChange) * sin(w*t); % true winding A control input
    ub = (1 + ControlChange) * cos(w*t); % true winding B control input
    xdot = [-Ra/L*x(1) + x(3)*lambda/L*sin(x(4)) + ua/L;
        -Ra/L*x(2) - x(3)*lambda/L*cos(x(4)) + ub/L;
        -3/2*lambda/J*x(1)*sin(x(4)) + 3/2*lambda/J*x(2)*cos(x(4)) - B/J*x(3);
        x(3)];
    x = x + xdot * dt;
    x(4) = mod(x(4), 2*pi);
    % Linear simulation
    w0 = -6.2832; % nominal rotor speed
    theta0 = -6.2835 * t + 2.3679; % nominal rotor position
    ia0 = 0.3708 * cos(2*pi*(t-1.36)); % nominal winding a current
    ib0 = -0.3708 * sin(2*pi*(t-1.36)); % nominal winding b current
    du = [ua - ua0; ub - ub0];
    F = [-Ra/L 0 lambda/L*sin(theta0) w0*lambda/L*cos(theta0);
        0 -Ra/L -lambda/L*cos(theta0) w0*lambda/L*sin(theta0);
        -3/2*lambda/J*sin(theta0) 3/2*lambda/J*cos(theta0) -B/J -3/2*lambda/J*(ia0*cos(theta0)+ib0*sin(theta0));
        0 0 1 0];
    G = [1/L 0; 0 1/L; 0 0; 0 0];
    dxdot = F * dx + G * du;
    dx = dx + dxdot * dt;
    xlin = [ia0; ib0; w0; theta0] + dx;
    xlin(4) = mod(xlin(4), 2*pi);
end

TCPU = cputime - tStart;
disp(['CPU time = ', num2str(TCPU)]);

% Compare linear and nonlinear simulation.
close all;
figure; hold on;
plot(tArray,xArray(1,:),'b-','LineWidth',1.0);
plot(tArray,xlinArray(1,:),'r--','LineWidth',1.0); 
set(gca,'FontSize',12); set(gcf,'Color','White'); set(gca,'Box','on');
xlabel('Seconds'); ylabel('Winding A Current (Amps)');
legend('Nonlinear', 'Linearized');

figure; hold on;
plot(tArray,xArray(2,:),'b-','LineWidth',1.0);
plot(tArray,xlinArray(2,:),'r--','LineWidth',1.0); 
set(gca,'FontSize',12); set(gcf,'Color','White'); set(gca,'Box','on');
xlabel('Seconds'); ylabel('Winding B Current (Amps)');
legend('Nonlinear', 'Linearized');

figure; hold on;
plot(tArray,xArray(3,:),'b-','LineWidth',1.0);
plot(tArray,xlinArray(3,:),'r--','LineWidth',1.0); 
set(gca,'FontSize',12); set(gcf,'Color','White'); set(gca,'Box','on');
xlabel('Seconds'); ylabel('Rotor Speed (Rad / Sec)');
legend('Nonlinear', 'Linearized');

figure; hold on;
plot(tArray,xArray(4,:),'b-','LineWidth',1.0);
plot(tArray,xlinArray(4,:),'r--','LineWidth',1.0); 
set(gca,'FontSize',12); set(gcf,'Color','White'); set(gca,'Box','on');
xlabel('Seconds'); ylabel('Rotor Position (Radians)');
legend('Nonlinear', 'Linearized');

% Put all four plots in a single figure
figure;
set(gcf,'Color','White'); 

subplot(2,2,1); hold on;
plot(tArray,xArray(3,:),'b-','LineWidth',1.5);
plot(tArray,xlinArray(3,:),'r--','LineWidth',1.5); 
set(gca,'FontSize',12); set(gca,'Box','on');
ylabel('Speed (Rad / Sec)');
legend('Nonlinear', 'Linearized');

subplot(2,2,2); hold on;
plot(tArray,xArray(4,:),'b-','LineWidth',1.5);
plot(tArray,xlinArray(4,:),'r--','LineWidth',1.5); 
set(gca,'FontSize',12); set(gca,'Box','on');
ylabel('Position (Radians)');

subplot(2,2,3); hold on;
plot(tArray,xArray(1,:),'b-','LineWidth',1.5);
plot(tArray,xlinArray(1,:),'r--','LineWidth',1.5); 
set(gca,'FontSize',12); set(gca,'Box','on');
xlabel('Seconds'); ylabel('Current A (Amps)');

subplot(2,2,4); hold on;
plot(tArray,xArray(2,:),'b-','LineWidth',1.5);
plot(tArray,xlinArray(2,:),'r--','LineWidth',1.5); 
set(gca,'FontSize',12); set(gca,'Box','on');
xlabel('Seconds'); ylabel('Current B (Amps)');