%% Main file to simualte the model as descibed by van den Herik et al. 2020
%  potato_main includes options to change loss_mode and timestep, as well
%  as the set of implemented differential equations
%Other scripts needed:
%  potato_init, defines model architecture and initial conditions
%  potato_parameters, defines the used parameters
%  potato_transport_phloem_xylem, calculates transport in the phloem and xylem
%  potato_loss_retieval, calculates loss, retrieval and removal from phloem and apoplast

%% 
close all; clear all;

%initialize variables
potato_init

%Initialize time step and endtime in seconds
tstep=1; tend=3600*24;
tspan= 0 : tstep : tend; 

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

%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);
toc 

%% function called by odesolver
function dy=odesys(t,y,loss_mode)
%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 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; 
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)
dy=dy';
end