%% Calculates transport in the phloem and xylem


%axial area and radius
A_ax_p  =   A_ax_p0 * exp((pp-0)/e_phloem);    
r_p     =   (A_ax_p/pi).^0.5;
%Conductivity and viscosity
psi     = sp*V_./(sp*V_+wp/rho);                      %volume fraction sucrose in phloem 
mu_p    = mu_x*exp((4.68*0.956*psi)./(1-0.956*psi));  %dynamic visc,MPa s
k_p     = (ap_p*A_ax_p)./mu_p;                        %phloem axial hydraulic conductivity, m4 MPa-1 s-1 
k_x     = (r_x^2)/(8*mu_x)* A_ax_x;                   %xylem axial hydraulic conductivity,  m4 MPa-1 s-1 
%Osmotic pressure and water potential
m       =   cp ./ (rho *(1- cp*V_));
pip     = - rho * R * T * (0.9988*m + 0.089*m.*m); %non-linear osmotic pressure, MPa
phi_st  =   pp + pip;                              %phloem water potential, MPa
phi_xyl =   px;                                     %xylem water potential, MPa
%Boundary conditions xylem           
JW_trans = -0.02/rho;
Jw_root  =   k_x*(px(end)-0)/x_un;       %water uptake, m3/s

%% Transport equations loading zone
i=2:N_load; 
JWp_in  = [0; k_p(i-1)/ x_load .* (pp(i-1)-pp(i)) ];           %Water flow phloem, m3/s, 
JWx_in  = [JW_trans; k_x/ x_load  .* (px(i-1)-px(i)) ];        %Water flow xylem, m3/s, 
JS_load = v_load*ones(N_load,1);                   %phloem loading rate  (mol/s)
A_rad_load   =   2*pi*r_p(1:N_load)*x_load  ;                  %Radial area phloem
JWlat_p = Lr*(phi_xyl(1:N_load)-phi_st(1:N_load)).*A_rad_load; %Lateral water flow, m3/s
%% Transport equations LD
i=N_load+1:N_load+N_LD; 
JWp_in  = [JWp_in;  k_p(i-1)/ x_LD .* (pp(i-1)-pp(i)) ];      %Water flow phloem, m3/s, 
JWx_in  = [JWx_in;  k_x/ x_LD .* (px(i-1)-px(i)) ];           %Water flow xylem, m3/s, 
JS_load     =   [JS_load ; zeros(N_LD,1)];                    %phloem loading rate  (mol/s)
A_rad_LD   =   2*pi*r_p(i)*x_LD  ;                            %Radial area phloem
JWlat_p = [JWlat_p ;Lr*(phi_xyl(i)-phi_st(i)).*A_rad_LD];     %Lateral water flow, m3/s
%% Transport equations Unloading zone
%Axial water flow
i=N_load+N_LD+1:N_load+N_LD+N_un; 
JWp_in  = [JWp_in; k_p(i-1)/ x_un .* (pp(i-1)-pp(i)) ];       %Water flow phloem, m3/s, 
JWp_out = [JWp_in(2:end); 0];                                 
JWx_in  = [JWx_in;  k_x/ x_un .* (px(i-1)-px(i)) ];           %Water flow xylem, m3/s, 
JWx_out = [JWx_in(2:end); Jw_root];
JSp_in      = JWp_in.*[0; cp((2:N_load+N_LD+N_un)-1)];        %Sucrose flow phloem, mol/s
JSp_out     = JWp_out.* cp ;        
JSp         = JSp_in-JSp_out;
s_unload    =   -v_unload * (cp(i) - 50);       %phloem unloading rate (mol/m/s)
JS_load     =   [JS_load; s_unload];

A_rad_un   =   2*pi*r_p(i)*x_un  ;   
JWlat_p = [JWlat_p; Lr*(phi_xyl(i)-phi_st(i)).*A_rad_un];     %Lateral water flow, m3/s
JWlat_x = [ -JWlat_p];    
