%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script to calculated phloem unloading parameters
% K. H. Jensen 1 Jun 2016.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Header
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function main
close all
clc
clear all
format short
warning off

% GLOBAL PARAMETERS
u=23e-6; % flow speed
r=3.6/2*1e-6; % SE radius
Lpd = 500e-9; % PD length
t=Lpd; % PD lentght
eta = 1.7e-3; % sap viscosity
rho = 1e3;
c = 500; % sap sugar concentration 
D=5e-10; % diffusivity of sucrose
e=75e-9; % funnel PD opening
L = 350e-6; % unloading zone lenth
A = 2*pi*r*L; % unloading zone wall area
h = 1e-5; % distance from PSE to PPP center
Lp=5e-14; % membrane permeability
Lpmin = 1e-14; % membrane permeability min range
Lpmax = 1e-12; % membrane permeability max range
ai = 7.5e-9; % pd inner radius
w = 2.8e-9; % simple pd slit width
ao=ai+w; % pd outer radius
N=240; % total number of PD

%% BASIC CALCULATIONS
disp('==== Basic calculations ====')
A1= pi*(ao.^2-ai.^2)
Q=pi*u*r^2
Q1= Q/N
I = Q*c
upd = Q./(N*A1)
Re = rho*upd.*w/eta

%% SCENARIOS
% TRANSPORT THROUGH SIMPLE PD
disp('==== Simple PD ====')

[dpInternalSimple dpEntranceSimple] = pressureCalculatorSimple(eta,u,r,N,ai,ao,w,Lpd)
dpInternalSimple+dpEntranceSimple


% TRANSPORT THROUGH FUNNEL PD
%disp('==== Funnel PD ====')
%[dpInternalFunnel dpEntranceFunnel] = pressureCalculatorFunnel(eta,u,r,N,ai,ao,w,Lpd,e)
                                      
% TRANSPORT THROUGH FUNNEL PD on PSE/PPP INTERFACE
disp('==== Funnel PD at PSE/PPP interface ====')
N=216;
[dpInternalFunnelPSEPPP dpEntranceFunnelPSEPPP] = pressureCalculatorFunnel(eta,u,r,N,ai,ao,w,Lpd,e)

% TRANSPORT THROUGH ANNULAR PD on PSE/PPP INTERFACE
disp('==== Funnel PD at PSE/PPP interface ====')
N=216;
[dpInternalAnnulus] = pressureCalculatorAnnulus(eta,u,r,N,ai,ao,w,Lpd,e)


% DIFFUSION
Nsimple=240;
Nfunnel=216;
disp('==== Diffusion ====')
%[deltacWall deltacPD deltacPDint deltacWallAndPDSimple deltacWallAndPDFunnel] = diffusionCalculator(u,r,D,ao,Nsimple,Nfunnel,c,A,h,t)
[deltacWall deltacPD deltacPDint deltacWallAndPDSimple deltacWallAndPDFunnel] = diffusionCalculator(u,r,D,ao,Nsimple,Nfunnel,c,A,h,t,e)

echo('test')
% MEMBRANE PRESSURE DIFFERENTIAL
disp('==== Membrane ====')
[deltaPmemb deltaPmembMin deltaPmembMax] = deltapMembrane(u,r,A,Lp,Lpmin,Lpmax)
end

function [dpInternalAnnulus] = pressureCalculatorAnnulus(eta,u,r,N,ai,ao,w,Lpd,e)
Q=pi*u*r^2;
Q1 = Q/N;
A1= pi*(ao.^2-ai.^2);
upd = Q./(N*A1);
x=linspace(0,Lpd,10000);
r2 = e-(e-ao)*x/Lpd;
r1 = ai*ones(1,10000);
dR = Q1*8*eta/pi*1./((r2.^2-r1.^2).*(r2.^2+r1.^2-(r2.^2-r1.^2)./log(r2./r1)));
dpInternalAnnulus = trapz(x,dR);
end




function [deltaPmemb deltaPmembMin deltaPmembMax] = deltapMembrane(u,r,A,Lp,Lpmin,Lpmax)
Q=pi*u*r^2;
deltaPmemb = Q/(A*Lp);
deltaPmembMin = Q/(A*Lpmin);
deltaPmembMax = Q/(A*Lpmax);
end

function [deltacWall deltacPD deltacPDint deltacWallAndPDSimple deltacWallAndPDFunnel] = diffusionCalculator(u,r,D,ao,Nsimple,Nfunnel,c,A,h,t,e)
rhoSimple = Nsimple/A;
rhoFunnel = Nfunnel/A;

Q=pi*u*r^2;
deltacWall = c*Q*h/(A*D);
deltacPD=c*Q./(2*Nsimple*D*ao);
deltacPDint = c*Q*t/(Nsimple*pi*ao^2*D);
deltacWallAndPDSimple = Q*c*h./(A*D)*(1+1./(2*rhoSimple*ao*h)*(1+t/ao*2/pi));
deltacWallAndPDFunnel = Q*c*h./(A*D)*(1+1./(2*rhoFunnel*ao*h)*(1+t/e*2/pi));
end

function [dpInternalSimple dpEntranceSimple] = pressureCalculatorSimple(eta,u,r,N,ai,ao,w,Lpd)
Q=pi*u*r^2;
A1= pi*(ao.^2-ai.^2);
upd = Q./(N*A1);
dpInternalSimple = 12*eta*upd*Lpd./(w.^2);
dpEntranceSimple = 32/pi*eta*upd./w;
end

function [dpInternalFunnel dpEntranceFunnel] = pressureCalculatorFunnel(eta,u,r,N,ai,ao,w,Lpd,e)
Q=pi*u*r^2;
Q1 = Q/N;
A1= pi*(ao.^2-ai.^2);
upd = Q./(N*A1);
%dpInternalFunnel = 8*eta*Q1*Lpd/e*1/(3*ao^3);
dpInternalFunnel = 8*eta*Q1*Lpd*(e^2+ao*(e+ao))/(3*pi*e^3*ao^3);
dpEntranceFunnel = 32/pi*eta*upd./w;
end



