function pep3

%% Parameters : 

% Plotting options : 
linecon = 2.5;   % linewidth for control diagram
line    = 2;     % linewidth for state-space 3D and 2D diagram

marker_size_arrow = 9;  % marker size (up arrow) in control diagram 
marker_size_conce = 20; % marker size for ini/fin concentrations

% Control :
alpha0  = 0;
alpha1  = 1;
alphain = 100;

% Density parameters : 
roTFA = 1489;                   % [g/L]
roACN = 786;                    % [g/L]

% Cube parameters :
maxPEP = 15;
maxTFA = 0.3;
minTFA = -0.005;
maxACN = 30;

% Initial conditions [Model] :
c10 = 1;                        % [g/L] PEP1
c20 = 0.1*roTFA/100;            % [g/L] TFA  -> B * roTFA / 100
c30 = 30*roACN/100;             % [g/L] ACN  -> E * roACn / 100   
V0  = 5;                        % [L]
A   = 0.0522;                   % [m2]

% Final conditions :
c1f = 5;                        % [g/L] PEP1
c2f = 0.02*roTFA/100;           % [g/L] TFA  -> B * roTFA / 100;
c3f = 3*roACN/100;              % [g/L] ACN  -> E * roACn / 100;

x0 = [c10 100*c20/roTFA 100*c30/roACN];
xf = [c1f 100*c2f/roTFA 100*c3f/roACN];

% Gather data:

% Time-Optimal :
[statetf timetf time1 time2 time3] = pep3_opt;
PEPopt = statetf(:,1);
TFAopt = 100*statetf(:,2)/roTFA;
ACNopt = 100*statetf(:,3)/roACN;
time1opt = time1; time2opt = time2; time3opt = time3;
a0opt = alpha0*ones(length(time1opt),1);
a1opt = alpha1*ones(length(time2opt),1);
aiopt = alphain*ones(length(time3opt),1);
timeopt = [time1opt;time2opt;time3opt];
aopt = [a0opt;a1opt;aiopt];

% C-CVD
[statetf time1 time2] = pep3_ccvd;
PEPccvd = statetf(:,1);
TFAccvd = 100*statetf(:,2)/roTFA;
ACNccvd = 100*statetf(:,3)/roACN;
time1ccvd = time1;
time2ccvd = time2;
a0ccvd = alpha0*ones(length(time1ccvd),1);
a1ccvd = alpha1*ones(length(time2ccvd),1);
timeccvd = [time1ccvd;time2ccvd];
accvd = [a0ccvd;a1ccvd];


% VVD
[timetf statetf alfa] = pep3_vvd;

PEPvvd = statetf(:,1);
TFAvvd = 100*statetf(:,2)/roTFA;
ACNvvd = 100*statetf(:,3)/roACN;
timetfvvd = timetf;
aVVD = alfa*ones(length(timetfvvd),1);

% Plot:

figure
plot3(PEPopt, TFAopt, ACNopt,'b',PEPccvd, TFAccvd, ACNccvd,'g:',PEPvvd, ...
      TFAvvd, ACNvvd,'r-.','LineWidth',line)  
hold on
plot3(x0(1),x0(2),x0(3),'go','Markersize',marker_size_conce,...
                             'MarkerEdgeColor','g',...
                             'LineWidth',line)
                         
plot3(xf(1),xf(2),xf(3),'rx','Markersize',marker_size_conce,...
                             'MarkerEdgeColor','r',...
                             'MarkerFaceColor','r',...
                             'LineWidth',line)
set(gca, 'box','on')
axis([0 maxPEP -0.005 maxTFA 0 maxACN])
xlabel('c1'); ylabel('c2'); zlabel('c3')
legend('Time-Optimal','C-CVD','VVD','Location','NorthEast');

set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 11 8.25])
print -depsc 'sol3_R1R2_states.eps'


figure
plot([time1opt;time2opt;time3opt],[a0opt;a1opt;aiopt],'b',[time1ccvd;time2ccvd],[a0ccvd;a1ccvd],'g:',timetfvvd, aVVD, 'r-.','LineWidth',linecon)

hold on
plot(time3opt(end), 1.15,'^','MarkerEdgeColor','b',...
                          'MarkerFaceColor','b',...
                          'markersize',marker_size_arrow)

% Options 
legend('Time-Optimal','C-CVD','VVD','Location','SouthEast')
xlabel('time'); ylabel('alpha')
axis([0 12 -0.05 1.18])
     
set(gcf, 'PaperUnits', 'centimeters') 
set(gcf, 'PaperPosition', [3 10 11 8.25])
print -depsc 'sol3_R1R2_control.eps'
     
