function [tgrid,xgrid,RHInt,SigInt] = IsothermalDryingST(Ops,...% cdts
                                      TherPar,PowPar,... % parameters
                                      depth,period,...   % depth & period
                                      nx,nt,...          % resolution
                                      flP)               % plot flag

% Calculates the time-dependent interstitial solvent vapor mass fraction Y 
% coupled to its counterpart Ye at equilibrium with the mass fraction of 
% solvent within grains per dry solid mass in a Smart Tray (ST) with 
% uniform powder depth. Solvent vapor is diffused through the powder and 
% is absorbed at a rate proportional to (Y-Ye). Ye is set by the powder's 
% isotherm, which relates the solvent mass fraction Om within the porous 
% solid powder and ratio RH of solvent partial and saturation pressure 
% Om=Om1 x RH + Om0. Temperature and pressure are invariant.
%
% The problem is rephrased in terms of Psi = (Y-Ys)/(Yi-Ys) and 
% Sig = (Ye-Ys)/(Yi-Ys), where Ys and Yi are, respectively, the vapor 
% mass fraction imposed at the free surface (for example in a moisture 
% chamber) and a reference interstitial vapor mass fraction within 
% the entire sample. 
%
% For a uniformly-wetted sample initially at equilibrium, Yi is the uniform 
% vapor mass fraction in the sample initially, and therefore Sig=Psi=1. 
%
% Time t is dim(-) with tAp and distance x from surf with sqrt{D tAp/tor}
% where tAp is the characteristic time detected by a vapor sorption 
% analyzer, if the powder conforms to the linear isotherm above. D is
% solvent vapor diffusion in the carrier gas at p and T; tor is tortuosity
% of the powder. 
%
% Michel Louge      MYL3@cornell.edu        Nov 30, 2023
%
% References:
%
% Louge, M. Y., Valance, A., Xu, J., Ould el-Moctar, A., & Chasle, P.(2022)
% Water vapor transport across an arid sand surface - Non-linear thermal 
% coupling, wind-driven pore advection, subsurface waves, and exchange with
% the atmospheric boundary layer. Journal of Geophysical Research: 
% Earth Surface, 127, e2021JF006490. https://doi.org/10.1029/2021JF006490
%
% Louge, M. Y., Mandur, J, Grigorov, P., Blincoe, W., Lamberto, D., Bower,
% C., & Meyer, R. F. (2024) Non-invasive, continuous, quantitative 
% detection of solvent content in vacuum drying, AAPS J., under review
%
%	input	-	Ops     =	row vector of operational parameters
%                           [RHi RHs T p] where
%                               RHi = uniform initial equilib powder RH (-)
%                               RHs = RH imposed at the free surface (-)
%                               T   = uniform temperature (K)
%                               p   = uniform pressure (Pa)
%           -   TherPar =   row vector of thermodynamic parameters
%                           [pA0 TAa TA0 Rha Dst nd Mwa Mca] where
%                               pA0 = Antoine pressure (Pa)
%                               TAa = Antoine activation temperature (K)
%                               TA0 = Antoine offset temperature (K)
%                               Rha = fundamental gas constant (J/mole.K)
%                               Dst = vapor diffusion at stand cdts (m^2/s)
%                               nd  = coefficient of diffusion temp dep (-)
%                               Mwa = molecular weight of solvent (kg/mole)
%                               Mca = mol weight of carrier gas (kg/mole)
%           -   PowPar  =   row vector of powder properties
%                           [rhp nu tor Om1 tAp] where
%                               rhp = particle material density (kg/m^3)
%                               nu  = powder solid volume fraction (-)
%                               tor = powder tortuosity (-)
%                               Om1 = slope isotherm Om vs RH at equil (-)
%                               tAp = apparent 1st-order kinetic time (sec)
%           -   depth   =   uniform powder thickness (m)
%           -   period  =   calculation period (sec)
%           -   nx      =   # evenly-spaced mesh pts along dist from surf
%           -   nt      =   # of equal time intervals in overall period
%           -   flP     =   if == 1, plot results; otherwise not
%
%	output 	-	tgrid   = grid of identical cols along which t (sec) varies
%           -   xgrid   = grid of ident rows with dist from free surf (m)
%           -   RHInt   = gridded RH linear interpolant, see example below
%           -   SigInt  = gridded Sig linear interpolant
%
%   global variable
%   
%           -   eta     =   Om1 rhp nu Rha T/[Mwa psat(1-nu)] where
%                           psat=pA0*exp[-TAa/(T-TA0)] = satur pres (Pa)
%           
% Example: define parameters, run prgm, make a surface plot of RH
%
% Ops =[0.5 0.65 273.15+50 1.013e5];
% TherPar = [1.24e10 3841.2 45.2 8.314 2.576e-5 0.861 0.018 0.029];
% tAp=424;PowPar = [1550 0.208 pi/2 0.11055 tAp];
% [tgrid,xgrid,RHInt] = IsothermalDryingST(Ops,TherPar,PowPar,...
%                                          0.008,30000,20,100,1);
%
% figure;surf(tgrid/3600,xgrid*1000,RHInt(tgrid,xgrid))
% ylabel('distance from free surface (mm)')
% xlabel('time (hr)')
%

global eta

%   internal variables (see above as well)
%
%           -   depdim  =   dim(-) depth = depth/sqrt{D tAp/tor}
%           -   Lre     =   ref length for distance x from free surface (m)
%           -   m       =   exponent in pdepe
%           -   Mra     =   ratio carrier mol weight/solvent mol weight
%           -   perdim  =   dim(-) period = period/tAp
%           -   psat    =   saturation pressure (Pa)
%           -   pst     =   standard pressure (Pa)
%           -   Sig     =   array (Ye-Ys)/(Yi-Ys) 
%                               (rows=fixed time;cols=fixed dist)
%           -   tdim    =   dim(-) row vec of times
%           -   Tst     =   standard temperature (K)
%           -   tstep   =   uniform time step (sec)
%           -   xdim    =   dim(-) row vec of distances from free surface
%           -   xstep   =   uniform depth resolution (m)
%           -   Ye      =   vapor mass frac at equil with solvent mass frac
%                           in the solid per dry solid mass
%

%% Unwrap input parameters
%
% Operational parameters
%
RHi = Ops(1) ; % initial equilibrium RH throughout the powder (-)
RHs = Ops(2) ; % imposed RH at the powder's free surface (-)
T   = Ops(3) ; % uniform temperature (K)
p   = Ops(4) ; % uniform pressure (Pa)
%
% thermodynamic parameters
%
pA0 = TherPar(1) ; % Antoine pressure (Pa)
TAa = TherPar(2) ; % Antoine activation temperature (K)
TA0 = TherPar(3) ; % Antoine offset temperature (K)
Rha = TherPar(4) ; % fundamental gas constant (J/mole.K)
Dst = TherPar(5) ; % vapor diffusion in carrier gas @ standard cdts (m^2/s)
nd  = TherPar(6) ; % temp depend exponent in rho D = rhost Dst (T/Tst)^nd
Mwa = TherPar(7) ; % solvent molecular weight (kg/mole), e.g. water
Mca = TherPar(8) ; % carrier gas molecular weight (kg/mole), e.g. air
%
pst = 1.013e5  	 ; % standard pressure (Pa)
Tst = 298.15     ; % standard temperature (298.15K)
%
% powder properties
%
rhp = PowPar(1) ; % material density of solid particles (kg/m^3)
nu  = PowPar(2) ; % powder solid volume fraction (-)
tor = PowPar(3) ; % powder tortuosity (-)
Om1 = PowPar(4) ; % slope of powder solvent isotherm (-)
tAp = PowPar(5) ; % apparent 1st-order kinetics time (sec)
%
%% Combine parameters
%
psat = pA0*exp(-TAa/(T-TA0))            ; % saturation pressure (Pa)
eta  = rhp*nu/(1-nu)*Om1/Mwa*Rha*T/psat ; % parameter eta (-)
D    = Dst*(pst/p)*(T/Tst)^(1+nd)       ; % diffusion coefficient (m^2/s)
Lre  = sqrt(D/tor*tAp)                  ; % reference length (m)
Mra  = Mca/Mwa                          ; % mol weight ratio carrier/solvnt
%
Yi   = RHi/(Mra*p/psat - RHi*(Mra-1))   ; % initial solvent vapor mass frac
Ys   = RHs/(Mra*p/psat - RHs*(Mra-1))   ; % surface solvent vapor mass frac
%
%% Resolution in time and depth
%
perdim = period/tAp ; % convert calculation period to dim(-) value
depdim = depth/Lre  ; % convert powder thickness to dim(-) value
%
tstep  = period/(nt-1) ; % time step (sec)
xstep  = depth/(nx-1)  ; % depth resolution (m)
%
[tgrid,xgrid] = ndgrid(0:tstep:period,0:xstep:depth) ; % interpolant grid
% tgrid contains t only; it has identical cols, along which t changes
% xgrid contains x only; it has identical rows, along which x changes
%
%% Solve PDEs
% see https://www.mathworks.com/help/matlab/ref/pdepe.html for details
%
m = 0                          ;   % coefficient for cartesian pdepe
xdim = linspace(0,depdim,nx)   ;   % define row vec mesh along depth (-)
tdim = linspace(0,perdim,nt)   ;   % define time mesh row vec (-)

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,xdim,tdim) ; % solve coupled PDEs
% sol(time,distance from surface,1=Psi or 2=Sigma)
% i.e., t varies along cols @ given x, while x varies along rows @ given t
%
%% Create interpolant for equilibrium RH vs time and dist from free surface
%
Sig = sol(:,:,2)        ;  % extract Sig(time,distance from free surface)
                            
Ye  = (Yi-Ys)*Sig+Ys    ;  % extract Y at equilib with interstitial RH
RH  = Mra*p/psat*Ye./(1+Ye*(Mra-1)) ; % interstitial eq. RH detected by ST

RHInt  = griddedInterpolant(tgrid,xgrid,RH)  ; % linear interpolant of RH
SigInt = griddedInterpolant(tgrid,xgrid,Sig) ; % linear interpolant of Sig

%% Plot results
%
    if flP == 1
scrsz   = get(0,'ScreenSize')   ;   % finds screen size of current machine
figure('Position',...
 [1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; % reserves space on screen for figure

surf(tgrid/3600,xgrid*1e3,RH)
title(['RH vs t and x with eta = ',num2str(eta),...
       ', Damkohler = ',num2str(depdim),...
       ', depth = ',num2str(1e3*depth),' mm',...
       ', critical depth = ',num2str(1e3*sqrt(2*D/tor*tAp/eta)),' mm'])
ylabel('distance from free surface (mm)')
xlabel('time (hr)')
view(-120,45)
    end

%% --------------------------------------------------------------
function [c,f,s] = pdex4pde(~,~,u,DuDx)

%	input	-	u       =	col vec [Psi;Sig] (-)
%           -   DuDx    =   col vec [dPsi/dx;dSig/dx] (-)
%
%	output	-	c       =   column vector in time derivative (-)
%           -   f       =   flux column vector (-)
%           -   s       =   source term column vector (-)
%
%   internal variables
%
%           -   dif     =   difference between Psi and Sigma 
%
% https://www.mathworks.com/help/matlab/ref/pdepe.html

global eta

c = [1; 1]              ;	% see pdepe help for details
f = [1; 0] .* DuDx      ;   % diffusion flux

dif = u(1) - u(2)       ;   % difference between Psi and Sigma
s = [-eta*dif; dif]     ;   % source term

%% --------------------------------------------------------------
function u0 = pdex4ic(~)

%	output	-	u0      =   initial col vec [Psi;Sig] (-)

u0 = [1;1] ; % initial cdts

%% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(~,ul,~,~,t)

%	input	-	ul      =	col vec [Psi;Sig] at free surface (-)
%           -   t       =   time/tAp (-)
%
%	output	-	pl      =   column vector for free surface BC (-)
%           -   ql      =   col vec for free surface BC (-)
%           -   pr      =   col vec for base BC (-)
%           -   qr      =   col vec for base BC (-)
%
% Boundary conditions: index l at the free surface; index r at the base
% see https://www.mathworks.com/help/matlab/ref/pdepe.html for details
%

pl = [ul(1); ul(2)-exp(-t)] ; % sets values @ surface Psi=0 & Sig=exp(-t)
%
% note that the BC for Psi=0 at the surface is at odds with a non-zero IC;
% this OK as the surface vapor mass fraction should instantaneously adjust.
ql = [0; 0]                 ;

pr = [0; 0]                 ;
qr = [1; 1]                 ;   % set the flux to zero at the base

