%% Using ND_PS_PencilE_Substep_v3
%
%%
%% PDE for 0<\tau_\mathrm{M}<\tau_1
% For the first time period 0<\tau_\mathrm{M}<\tau_1 we have no substrate in the bulk i.e \mathrm{[S]_{bulk}}. We are unable at this point to solve the coupled PDEs such that initial condition for the substrate can be set equals to 0. Which is why we remove the substrate bit from the coupled PDEs. This ends up looking as (this has been double checked in the word file called "Simulation Plan v2") under and it is this single PDE that we solve with no flux BC at the RHS and Nersntian concentration at the LHS.
% \frac{\partial m_{ox}(\chi, \tau_{m})}{\partial \tau_m}=\frac{\partial^2 m_{ox}(\chi, \tau_m)}{\partial \chi^2}
%
%% System of PDEs for \tau_1<\tau_\mathrm{M}<\tau_2
% \frac{\partial s(\chi, \tau_{m})}{\partial \tau_m}=\frac{1}{\alpha}\frac{\partial^2 s(\chi, \tau_m)}{\partial \chi^2}-\frac{\gamma\kappa^2 m_{ox}(\chi,\tau_m) s(\chi, \tau_m)}{\eta\alpha\gamma m_{ox}(\chi,\tau_m)(1+\mu s(\chi,\tau_m))+\alpha s(\chi,\tau_m)}
% \frac{\partial m_{ox}(\chi, \tau_{m})}{\partial \tau_m}=\frac{\partial^2 m_{ox}(\chi, \tau_m)}{\partial \chi^2}-\frac{\kappa^2 m_{ox}(\chi,\tau_m) s(\chi, \tau_m)}{\gamma m_{ox}(\chi,\tau_m)(1+\mu s(\chi,\tau_m))+s(\chi,\tau_m)}
%% For Initial Condition 0<\tau_\mathrm{M}<\tau_1
% s(\chi,0)=0                    The initial condition is there but the equation it supports does not do anything. Ideally it should be set to 0 to better represent the physical system but this is not incorrect either as it feeds the ones to the next pdefun call.
% m_{ox}(\chi,0)=0
%% For Initial Condition \tau_1<\tau_\mathrm{M}<\tau_2
% s(\chi,0)=s(\chi,\tau_{1})
% m_{ox}(\chi,0)=m_{ox}(\chi,\tau_{1})
%% And Boundary Condition 0<\tau_\mathrm{M}<\tau_1
% Substrate L.H.S            \frac{\partial s(0,0<\tau_\mathrm{M}<\tau_{1})}{\partial \chi}=0
% Substrate R.H.S           s(1,0<\tau_\mathrm{M}<\tau_{1})=1
% Mediator R.H.S             \frac{\partial m_{ox}(1,0<\tau_\mathrm{M}<\tau_{1})}{\partial \chi}=0
%% Mediator L.H.S 0<\tau_\mathrm{M}<\tau_1
% Mediator L.H.S              m_{ox}(0,0<\tau_\mathrm{M}<\tau_{1}))=m_{\epsilon} where, Value on the R.H.S is a constant value.
%% Mediator L.H.S \tau_1<\tau_\mathrm{M}<\tau_2 for new \mu and \gamma
% Mediator L.H.S              m_{ox}(0,(0,\tau_{1}<\tau_\mathrm{M}<\tau_{2}))=m_{\epsilon}
%

function [I_num, tspan, c_s, c_mox, c_red, chi1] = ND_PS_EK_PencilE_SubStep_v3(k, N1, N2, m_layer, Area, tspan1, tspan2, l, D_m, I_Exp1, I_Exp2, callVideo)


% Legend
% k(1)=kappa
% k(2)=eta
% k(3)=gamma1
% k(4)=mu1
% k(5)=alpha

% Constants Used in the Nernst Equation   
n           = 2;                                    % No. of Electrons Involved
F           = 96485.3329;                           % C/mol
R           = 8.314;                                % J/(mol.K) or CV/(mol.K)
Temp        = 273.15+21;                            % K
m           = 0;                                    % Cartesian Co-ordinates
N           = N1;
color       = 'blue';

% FIRST Glucose STEP    
chi1        = linspace(0, 1, N1); 
tau_M1      = (D_m*tspan1(end))/l^2; 
tau1        = linspace(0, tau_M1, N1);                              % Dimensionless Time 
E0          = 0.180;                                                % V (Vs. SCE Reference)
epsilon1    = ((n*F)/(R*Temp))*(0.35-E0);                           % Dimensionless Potential 
c_M_ox      = 1/(1+exp(-epsilon1));                                 % Mox BC 
ic_arg_1    = {@(x)zeros(size(N1)) ; @(x)zeros(size(N1))};
IC          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC          =  @(xl, yl, xr, yr, t)PDE_ND_PS_EK_BC1(xl, yl, xr, yr, c_M_ox); 

sol1        = pdepe(m, @(x, t, c, DcDx)PDE_PS_EK_StubStep1(x, t, c, DcDx), ... 
                    IC, BC, chi1, tau1);
               
% Concentration Profiles c(i, j, k)(Solutions) for FIRST Glucose STEP 
c1            = sol1(:, :, 1);                                           % Substrate Conc.  
c2            = sol1(:, :, 2);                                           % Mox Conc.
c3            = 1 - c2;                                                  % Mred Conc. 

% Calculation of Dimensionless Mox Flux at x=0
for counter   = 2:N1
    [~, J_num_M(counter)]  = pdeval(m, chi1, c2(counter,:),  0);          % ND M Flux at chi=0
    [~, J_num_S(counter)]  = pdeval(m, chi1, c1(counter,:),  1);          % ND S Flux at chi=1
end

j_num_1       = -(J_num_M.*m_layer*D_m*1)/l;
I_num1        = n * F * Area .* j_num_1; 


% SECOND Glucose STEP 
chi2        = linspace(0, 1, N2); 
tau_M2      = (D_m*tspan2(end))/l^2; 
tau2        = linspace(tau_M1, tau_M2, N2);      
ic_arg_2    = {@(x)zeros(size(N1)) ; @(x)interp1(chi2, sol1(N1, :, 2), x)};
IC2         =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_2); 
BC          =  @(xl, yl, xr, yr, t)PDE_ND_PS_EK_BC2(xl, yl, xr, yr, c_M_ox); 

sol2        = pdepe(m, @(x, t, c, DcDx)PDE_PS_EK_StubStep2(x, t, c, DcDx, k), ... 
                    IC2, BC, chi2, tau2);

% Concentration Profiles Step-2 c(i, j, k)(Solutions) for SECOND Glucose STEP 
c14         = sol2(:, :, 1);                                             % Substrate Conc.  
c25         = sol2(:, :, 2);                                             % Mox Conc.
c36         = 1 - c25;                                                    

% Calculation of Dimensionless Mox Flux at x=0
for counter   = 2:N2
    [~, J_num_M_02(counter)]  = pdeval(m, chi2, c25(counter,:),  0);      % ND M Flux at chi=0
    [~, J_num_S_l2(counter)]  = pdeval(m, chi2, c14(counter,:),  1);      % ND S Flux at chi=1
end

j_num_2       = -(J_num_M_02.*m_layer*D_m*1)/l;
I_num2        = n * F * Area .* j_num_2; 
I_num         = [I_num1, I_num2]; 
I_num         = I_num'; 
tspan1 = linspace(0, round(((tau_M1*(l^2))/D_m)), N1); 
tspan2 = linspace(round(((tau_M1*(l^2))/D_m)), round(((tau_M2*l^2)/D_m),2), N2); 
tspan  = [tspan1,tspan2];
c_mox  = [c2; c25];
c_red  = [c3; c36];
c_s    = [c1; c14];

if callVideo
    [myVideo] = Video(k, N, chi1, c_s, c_mox, c_red, tspan, I_num);
end

% This is the case finder function for Step-1
[J_ana_ND, I_ana, Case_state] = case_finder2(k, n, F, c_M_ox, m_layer, l, D_m, Area); 
% Outputs 
fprintf('Case_State after Glucose Addition: %d \n',Case_state) 
fprintf('log(kappa) is %f \n',log10(k(1))) 
fprintf('log(eta) is %f \n',log10(k(2)))
fprintf('log(gamma_1) is %f \n',log10(k(3)))
fprintf('log(mu_1) is %f \n',log10(k(4)))

z = 1;

end

%

function [cc, ff, ss] = PDE_PS_EK_StubStep1(x, t, c, DcDx)
%                      S;             Mox;   
cc              = [    1;              1]; 

ff              = [    1;              1].*DcDx;

ss              = [    0;              0]; 

end

%% pdepe Function after glucose additon 

function [cc, ff, ss] = PDE_PS_EK_StubStep2(x, t, c, DcDx, k)

%                      S;          Mox;    
cc              = [ k(5);           1]; 

ff              = [    1;           1].*DcDx;
S_kin           = ((k(3)/k(2))*(k(1)^2)*c(1)*c(2))./...
                  (k(3).*c(2).*(1+(k(4).*c(1)))+c(1));
M_kin           = ((k(1)^2).*c(1).*c(2))./...
                  (k(3).*c(2).*(1+(k(4).*c(1)))+c(1));
ss              = [-S_kin;  -M_kin]; 

end


%% Initial Condition --> Initial Concentrations of Species 

function c0     = PDE_ND_PS_EK_IC(x, ic_arg_1)
   
%                         Substrate;               Mred;  
c0              = [ic_arg_1{1}(x).';  ic_arg_1{2}(x).'];

end

%% Boundary Condition-1

function [pl, ql, pr, qr] = PDE_ND_PS_EK_BC1(xl, yl, xr, yr, c_M_ox)

%                Substrate;         Mediator;
pl     = [0               ;    yl(2)-c_M_ox];     

ql     = [1               ;               0];

pr     = [yr(1)-0         ;               0];

qr     = [0               ;               1]; 


end

%% Boundary Condition-2

function [pl, ql, pr, qr] = PDE_ND_PS_EK_BC2(xl, yl, xr, yr, c_M_ox)

pl     = [0               ;    yl(2)-c_M_ox];     

ql     = [1               ;               0];

pr     = [yr(1)-1         ;               0];

qr     = [0               ;               1]; 


end

%% Case Finder Function for Substrate Step 1

function [J_ana_ND, I_ana, Case_state] = case_finder2(k, n, F, c_M_ox, m_layer, l, D_m, Area)

    ae1         = c_M_ox;
    
    % Case-1
    if k(1) < 1 && k(3)*(1+k(4)) < 1
        J_ana_ND  = (k(1).^2)*ae1; 
        Case_state  = 1; 
    % Case-2
    elseif k(1) >= 1 && k(1) < 1+(k(2)/k(3)) && k(1) > k(3)*(1+k(4))
        J_ana_ND  = k(1)*ae1;
        Case_state  = 2; 
    % Case-3
    elseif k(1) > 1+(k(2)/k(3)) && k(1)*sqrt(k(2)) > (1+k(3))*(sqrt(1+0.5*k(4)))
        J_ana_ND  = ae1+(k(2)/k(3));
        Case_state  = 3; 
    % Case-4
    elseif k(1)*sqrt(k(2)) < (1+k(3))*(sqrt(1+0.5*k(4))) && k(1) > ((1+k(4))*(sqrt(k(2))))/(sqrt(1+0.5*k(4)))
        J_ana_ND  = (sqrt(k(2)).*k(1))./(k(3).*sqrt(1+(0.5*k(4))));
        Case_state  = 4;    
    % Case-5
    elseif k(3)*(1+k(4)) > 1 && k(1) < ((1+k(4))*(sqrt(k(2))))/(sqrt(1+0.5*k(4))) && k(1) < k(3)*(1+k(4))
        J_ana_ND  = (k(1).^2)/(k(3)*(1+k(4))); 
        Case_state  = 5; 
    % Case-6
    elseif k(1)>=sqrt(2*k(2)/k(3))&&k(3)>=2*k(2)&&2*k(1)^2<=1+(k(3)/k(2))+(k(2)/k(3))&&...
            k(3)*(1+0.5*k(4))<=0.5&&k(3)<=(k(1)^2)/(2*k(2)*(1+k(4))^2)
        J_ana_ND  = sqrt(2*ae1*k(1)^2*k(2)/k(3)); 
        Case_state  = 6; 
    % Case-7
    elseif k(1)>=sqrt(2/(k(3)+k(3)*k(4)))&&k(3)*(1+k(4))>=2&&(k(3)+k(2))>=k(1)*sqrt((2*k(3))/(1+k(4)))&&...
        k(2)*(1+k(4))>=2*k(3)*(1+0.5*k(4))&&k(1)^2>=(2*k(3))*(1+k(4))
    J_ana_ND  = sqrt((2*ae1*k(1)^2)/(k(3)*(1+k(4)))); 
    Case_state  = 7; 
    else
        disp("Get Outta Here");
        Case_state  = 'nothing to see here move on lads'; 
    end
    
    j_ana       = (J_ana_ND*m_layer*D_m*1)/l; 
    I_ana       = 2 * 96485.3329 * Area .* j_ana;          
             
end

%% Video Output Function 

function [myVideo] = Video(k, N, chi1, c_s, c_mox, c_red, tspan, I_num)

% Initialize Video
G                 = figure(1); 
myVideo           = VideoWriter(sprintf('κ%.2f η%.2f γ%.2f µ%.2f α%.2f', k(1), k(2), k(3), k(4), k(5))); 
myVideo.FrameRate = 10;                                % can adjust this, 5 - 10 works well for me
myVideo.Quality   = 100; 
open(myVideo);

u = uicontrol(G, 'Style','slider','Position',[0 40 10 360],...
             'Min', 2,'Max', 2*N,'Value', 2);


for ii = 2:2*N
    
    % Plot Species Conc. in the Layer 
    subplot(2,1,1);
    yyaxis left
    plot(chi1,c_mox(ii,:),chi1,c_red(ii,:), 'LineWidth', 2);
    ylabel('{\it m}')                                             
    ylim([0,1]);
    yyaxis right 
    plot(chi1, c_s(ii,:), 'LineWidth', 2);
    ylabel('{\it s}')                                            
    ylim([0,1]);
    xlabel('\chi'); 
    legend('{\it m}_{ox}','{\it m}_{red}', '{\it s}'); 
    title('Non-Dimensionalized Conc.');
    
    subplot(2,1,2);
    addpoints(animatedline(tspan,I_num,'marker', '.', 'markersize', 9,...
    'linestyle', '--', 'MaximumNumPoints', 1), tspan(ii), I_num(ii)*1e06);
%     ylim([0,10]);
    title('Dimensionalized Response');
    text(8, 8,{'',''})
    drawnow 
    xlabel('{\it t} / s'); ylabel('{\it i} / \muA');
    

    u.Value = ii; 
    uText   = uicontrol('Style','Edit','Position',[0,00,40,40], ...
                        'String',num2str(tspan(ii),3));
    uText   = uicontrol('Style','Edit','Position',[400,160,110,30], ...
        'String',sprintf('Response = %.2f %s',I_num(ii-1)*1e06,'µA'));

    pause(0.001)
    M      = getframe(G);  
    writeVideo(myVideo, M);
     
end

end
