close all; clear all;

%Time step and endtime in seconds
tstep=10; tend=3600*24;
tspan= 0 : tstep : tend;

%Set leakage-retrieval-mode
%1 = No-Leakage, 2 = loss-only, 3 = loss-Retrieval
loss_mode = 1;

%Xylem inclusion?
%0 = No Xylem, 1 = Xylem
Xyl = 0;

%Loading rate range (nmol/s)
v_load_ =  linspace(0.1, 25,25)*1e-9;

%Characteristics
sc.Vmax1 = 22.5e-9; 
sc.Vmax2 = 2.5e-9;
sc.Km1   = 75;
sc.Km2   = 75;
sc.l1    = 0.25;
sc.l2    = 0.25;

%Storage unloading rates
unload_sink_1 =[]; unload_sink_2 =[];

%initialize architecture
potato_init

for i = 1:length(v_load_)


%Solve ODE system defined in function odesys
opt=odeset('RelTol', 1e-13, 'AbsTol', 1e-14);
tic
[t,y]=ode15s(@odesys, tspan, C_ini,opt,loss_mode, Xyl,v_load_(i),sc);
toc

%calculate unloading rates both sinks
wp=y(end,phloem_c_water)'; sp=y(end,phloem_c_sucrose)';
pp=y(end,phloem_c_pressure)';vp=wp/rho; cp=sp./vp;

unload_1  = -sc.Vmax1 / N_1_un .* (cp(Loc_1_un,:)./( cp(Loc_1_un,:) + sc.Km1));
unload_2  = -sc.Vmax2 / N_2_un .* (cp(Loc_2_un,:)./( cp(Loc_2_un,:) + sc.Km2));

unload_sink_1 = [unload_sink_1 -sum(unload_1(:,[ end]))'];
unload_sink_2 = [unload_sink_2 -sum(unload_2(:,[ end]))'];

end

% create partitioning coefficient plot
PC1 = unload_sink_1 ./ (unload_sink_1+unload_sink_2) * 100;

plot(v_load_*1e9, PC1, 'Linewidth',1.5);
xlabel('Loading rate (nmol/s)'); ylabel('PC1 (%)');
ylim([50 100]);

%% function called by odesolver
function dy=odesys(t,y,loss_mode,Xyl, vload, sc)
%Update variables every timestep
    potato_init;

    %phloem
    wp=y(phloem_c_water); sp=y(phloem_c_sucrose);
    pp=y(phloem_c_pressure);vp=wp/rho; cp=sp./vp;
    %xylem
    wx=y(xylem_c_water);px=y(xylem_c_pressure);
    %apoplast
    wpa=y(phloem_a_water);spa=y(phloem_a_sucrose);
    vpa=wpa/rho; cpa=spa./vpa;

%Calcualte basic phloem and xylem transport rates
    potato_transport_phloem_xylem

%Calculate unloading at sinks
    potato_unloading

%Calculate leakage-retrieval and apoplast removal rates along the stem
    potato_loss_retrieval

%% Differential equations
dy(phloem_c_water)     = (JWp_in - JWp_out + JWlat_p)*rho;
dy(phloem_c_sucrose)   = JSp + JS_lat  + JS_load;
dy(phloem_c_pressure)  = e_phloem * 1./(wp) .* dy(phloem_c_water)';
dy(phloem_a_sucrose)   = -JS_lat - Resp_apo;

if Xyl == 0;
dy(xylem_c_water)      = 0*(JWx_in - JWx_out + JWlat_x)*rho;  %water xylem (g/s)
dy(xylem_c_sucrose)    = 0;
dy(xylem_c_pressure)   = 0*e_xylem * 1./(wx) .* dy(xylem_c_water)';%pressure xylem (MPa/s)
elseif Xyl == 1;
dy(xylem_c_water)      = (JWx_in - JWx_out + JWlat_x)*rho;  %water xylem (g/s)
dy(xylem_c_sucrose)    = 0;
dy(xylem_c_pressure)   = e_xylem * 1./(wx) .* dy(xylem_c_water)';%pressure xylem (MPa/s)
end

dy=dy';
end
