```function res = sensn(F,ploton,dispon)
% SENSN simulates transfer through n layers (generalization of sens to n layers)
%
% PHYSICAL FORMULATION
% pi = ki.Ci (Ci in mass/volume)
% p1=p2= => K12 = C1/C2 = k2/k1 at equilibrium
%
% continuity at interface: D1.r1.dC1/dx = D2.r2.dC2/dx
% => D1*.dp1/dx = D2*.dp2/dx with Di* = (ri/ki).Di
%
% transport eq: dCi/dt = Di.d2Ci/dx2
% with: Di* = Di/ki and pi = ki.Ci:
% ri/ki.dpi/dt = Di*.d2pi/dx2
%
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% The current formulation (june 2006) assumes:
%   that all conc are expressed in mass/volume (density already corrected)
%     e.g. p = keq.Ceq with keq = k/r and Ceq = C*r/C0max
%   the reference layer is chosen as the layer with the lowest Di*/li value (mass transport resistance)
%   the equivalent dilution factor respectively to the reference layer must be found with the following approach
%       crit   = inline('sum(r.*(l0+cumsum(l)).^(g+1)-(l0+[0 cumsum(l(1:end-1))]).^(g+1))/l0-L','l0','l','r','g','L');
%       l0     = lref*fzero(crit,100,[],l/lref,r,g,L);
%           where  L is the mass dilution factor between P and F
%                  l0 is the equivalent thickness for the food
%                  g is the geometry number (0=cartesian, 1=cylindrical, 2=spherical)
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% SAFE FOOD PACKAGING PORTAL 1.0 - 17/01/06 - INRA\Olivier Vitrac - rev. 28/06/06
%
%
% Follow this link to see the scaled dimensionless formulation (compare with sens.m to compare with the monolayer code).
% We do not provide the code of pdesolve1D.m It is a general PDE 1D solver (derived from the one used in the PDE toolbox, Matlab v. 5.3)
% Contact us, if you want a copy of pdesolve1D.m (you can use equivalently parabolic.m in the current PDE toolbox, Note however that the macroscopic mass balance is not performed as required in the BC)
% definitions
global timeout
timeout		= 300; % s
options		= odeset('RelTol',1e-4,'AbsTol',1e-4,'Stats','yes','Initialstep',1e-5,'Maxstep',.05,'Maxorder',5);
Fdefault	= 	struct(...
'Bi'		, 	1e3,...	Biot [hm.L1/D]
'k'			,	[1 1 1 1],...[0.5 3 2],...	ki, i=1 (layer in contact with the liquid)
'D'         ,   [1e-16 1e-14 1e-14 1e-14],... diffusion coefficient
'k0'        ,   1,... 0 = liquid
'l'         ,   [50 20 10 120]*1e-6,...[50 20 10 120]*1e-6,... m
'L'			,	200/1800,...	dilution factor (respectively to iref)
'C0'		,	[0 500 500 500],...	initial concentration in each layer
'options'	,	options...
); % if iref is mission, it is indentified
Fdefault.t	= 0:.00001:.005; %[0:.00001:.005 .01:.01:.1 .11:.1:5]';
method		= 'cubic';
BiKminode	= 0.1;
MaxStepmax	= 10;
nstepchoice = 200;
n_default	= 1e4;
if Fdefault.Bi<10, Fdefault.t = Fdefault.t/(Fdefault.Bi/10); end
ploton_default = false;
dispon_default = false;
nmesh_default  = 50; % number of nodes for a layer of normalized thickness 1
geometry_default = 0; % slab

% arg check
initon = false;
if ~nargin, initon = true; end
if nargin<1, F = []; end
if nargin<2, ploton = []; end
if nargin<3, dispon = []; end
if isempty(F), F = Fdefault; end
if ~isfield(F,'autotime'), F.autotime = 'on'; end
if ~isfield(F,'n'), F.n = n_default; end
if ~isfield(F,'geometry'), F.geometry = geometry_default; end
if ~isfield(F,'nmesh'), F.nmesh = nmesh_default; end
if ~nargout, ploton=true; dispon=true; end
if isempty(ploton), ploton_default; end
if isempty(dispon), dispon_default; end

% physical check
m = Inf;
for prop = {'D' 'k' 'C0'};
if strcmp(prop{1},'C0')
F.(prop{1}) = F.(prop{1})(F.(prop{1})>=0);
else
F.(prop{1}) = F.(prop{1})(F.(prop{1})>0);
end
m = min(m,length(F.(prop{1})));
end
F.m = m;
% init
if initon
res = Fdefault;
if dispon, disp('... init mode'), end
return
end
if strcmp(lower(F.autotime),'on')
ti		= linspace(min(F.t),max(F.t),F.n)';
else
ti		= F.t;
end
MaxStepmax = max(MaxStepmax,ceil(F.t(end)/nstepchoice));

% renormalization (fields X1->Xn)
F.k = F.k/F.k0; F.k0 = 1;   % by convention
a = F.D(1:m)./F.k(1:m);
if isfield(F,'iref')
iref = F.iref;
else
[crit,iref] = min( a./F.l(1:m) );
F.iref = iref;
end
F.a = a./a(iref);
F.lref  = F.l(1:m)./F.l(iref);
F.lrefc = cumsum(F.lref);
F.C0    = F.C0(1:m);
F.l     = F.l(1:m);
F.D     = F.D(1:m);

% mesh generation
xmesh = unique([linspace(0,sum(F.lref),F.lrefc(end)*F.nmesh*F.m) F.lrefc])'; %mesh1D(300,[0 1],[.4 .2]); %

% update mesh to geometry
if F.geometry>0 && xmesh(1)==0
pos0 = F.L;
F.lrefc = F.lrefc + pos0;
xmesh = xmesh + pos0;
else
pos0 = 0;
end
F.pos0 = pos0;
F.pos = [pos0-eps F.lrefc]; % only used for extrapolation

% equilibrium concentration
C0eq = trapz(xmesh,F.C0(ilayer(F,xmesh))'.*xmesh.^F.geometry)/( (1/F.L).^(F.geometry+1) + trapz(xmesh,1./F.k(ilayer(F,xmesh))'.*xmesh.^F.geometry) );
F.peq = F.k0 * C0eq;

% solving
if dispon, disp('SENSN: solving...'), end
errpde = false;
tic, [sol,Fca,Fcb] = pdesolve1D(F.geometry,@pdediff_PDE,@pdediff_IC,@pdediff_BC,xmesh,F.t,F.options,F);
F.x = xmesh;

% extraction of the solution
if ~errpde
if dispon, disp(['SENSN: end in ' num2str(toc) ' s']), end
nt      = length(F.t);
p		= sol(:,:,1)*F.peq;
if F.m>1
C   = p./repmat(F.k(ilayer(F,xmesh)),nt,1);
else
C   = p./F.k;
end
V       = trapz(xmesh,xmesh.^F.geometry);
V0      = (1/F.L).^(F.geometry+1);
Cm		= trapz(F.x,C.*repmat(xmesh',size(C,1),1).^F.geometry,2)/V;
C0       = trapz(xmesh,F.C0(ilayer(F,xmesh))'.*xmesh.^F.geometry)/V;
C1		= C(:,1);
fc		= pdeloss(F.x,C,0,F.geometry);
Cmi		= interp1(F.t,Cm,ti,method);
fci		= interp1(F.t,fc,ti,method);
C1i		= interp1(F.t,C1,ti,method);
p1i		= interp1(F.t,p(:,1),ti,method);
fi		= F.Bi*p1i-F.Bi*(F.L/F.lrefc(end))*fci; % S.Bi*((Ua-0) - S.L * -Fc(1) )
if fi(end)<0, fi = fi-min(fi); end
res.p = p;
res.peq = F.peq;
res.V  = V;
res.C0 = C0;
res.F  = F;
res.x  = F.x-pos0;
res.Cx = C;
res.tC = F.t;
res.t  = ti;
res.C  = Cmi;
res.CF = (res.C0-res.C)*V/V0;
res.fc = fci;
res.f  = fi;
res.fD = fi/F.Bi;

else
if dispon, disp(['SENSN is TIMEOUT => check your parameters']), end
res = [];
end

% plots
if ploton
figure(1)
subplot(231), hold on, plot(F.x,C), xlabel('x'), ylabel('C')
subplot(232), hold on, plot(ti,Cmi,'b-'), xlabel('Fo')
subplot(233), hold on, plot(sqrt(ti),Cmi,'b-'), xlabel('Fo^{1/2}')
subplot(234), hold on, plot(ti,fi,'b-'), ylabel('j')
subplot(236), hold on, plot(Cmi,fi,'b-'), xlabel('C'), ylabel('j')
figure(2)
subplot(121), hold on, plot(ti/F.Bi,Cmi,'k:'), xlabel('Fo/Bi'), ylabel('C')
subplot(122), hold on, plot(Cmi,1/F.Bi*fi,'k:'), xlabel('C'), ylabel('j/Bi')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% P H Y S I C A L   F O R M U L A T I O N
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PDE problem :  [C,F,S] = PDEFUN(X,T,U,DUDX,P1,P2...)
function [c,f,s] = pdediff_PDE(x,t,U,dUdx,S)
global timeout
if toc>timeout, eval('erron'), end
i = find(floor(S.pos-x)<0,1,'last'); %i = find(floor([-eps S.lrefc]-x)<0,1,'last');
c = 1/S.k(i);
f = S.a(i)*dUdx;
s = 0;

% IC : U = ICFUN(X,P1,P2...)
function U0 = pdediff_IC(x,S)
%i = min(floor(interp1q(cumsum([0 S.lref])',(0:S.m)',x-S.lrefc(1)))+1,S.m); % OLD LINE
% i = min(floor(interp1q(cumsum([0 S.lref])',(0:S.m)',x-S.pos0))+1,S.m); % NEW LINE (24/02/06 OV)
i = ilayer(S,x);
U0 = S.C0(i).*S.k(i)/S.peq;
if size(U0,1)==size(x,2), U0 = U0'; end

% BC : [PL,QL,PR,QR] = BCFUN(XL,UL,XR,UR,T,P1,P2...)
function [pa,qa,pb,qb,purge] = pdediff_BC(xa,Ua,xb,Ub,t,Fc,S)
global Hsmooth
% pa = 0;
% qa = 1;
% pb = S.Bi*((Ub-0) - S.L * - Fc(2) );
% qb = 1;
pa = S.Bi*((Ua-0) - S.L * -Fc(1) );   % K=S.k(1)/S.k0 with S.k0=1
qa = -1;
pb = 0;
qb = 1;

% General 1D PDE formulation (without advection)%
%  c(x,t,u,Du/Dx) * Du/Dt = x^(-m) * D(x^m * f(x,t,u,Du/Dx))/Dx + s(x,t,u,Du/Dx)
%  --------------                            --------------       --------------
%   capacitance                                   flux                source
% (diagonal positive nxn matrix)                nx1 vector           nx1 vector
%
%  u(t,x) can be defined in R^n on a<=x<=b during t0 <= t <= tf
%  m = 0, 1, or 2 => {slab, cylindrical, or spherical} symmetry
%  subjected to BC:  p(x,t,u) + q(x,t) * f(x,t,u,Du/Dx) = 0
%                               -------
%                           diagonal nxn matrix
%
```
`CONTACT US FOR ANY REQUESTION, BUG REPORT`

SAFE FOOD PACKAGING PORTAL
Date: 26-June-2006