%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Diffusion of dyes through porous medium by molecular diffusion and 
% scarce-range phenomenon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Last at 03:30 am on Friday, 19 of January 2024
%
% 1 - Problem Analysis
%----------------------------------------------------------------------
% Assumptions:
%       --Let's consider specifically the penetration or diffusion of an 
% artificial or natural-based dye, such as erythrosine and red gardenia in 
% a porous media, namely the cherry's flesh.
%       --It was deemed appropriate to assume that the accumulation of 
% solutes in the skin, represented by (r_S-r_O), and accordingly the 
% resistance to species transfer, was insignificant in comparison to the 
% flux through the porous surface and the flux itself, given by the cherry 
% skin thickness of 0.0004 m. The cherries’ flesh was considered 
% homogeneous and isotropic, and the diffusion process exhibited radial 
% symmetry with unidirectional diffusion. Due to the great agitation of the 
% sweetener solution surrounding the internal and external surfaces, the 
% concentration in these areas immediately acquires the concentration of 
% the liquid. Furthermore, the thickness of the epidermis was found to be 
% much smaller than the radii of the inner and outer faces, thus the 
% difference in surface areas could be disregarded and the thickness 
% assumed as (r_S-r_O).
%       --By fitting a trend line to the redness parameter, a*, in 
% comparison to the present concentration of erythrosine or red gardenia 
% dye, an alternative version of equation (2) can be obtained, equation 
% (3), thereby establishing a correlation that enabled the experimental 
% data to be fitted to the general equation governing species transport.
% [da/dtheta=((2/R)*(da/dR))+((d^2a)/(dR^2))]                                                                                                                                  (3)
% Subjected to the following initial and boundary conditions:
% I.C.:  At theta=0; a==1;                                      at A<R<1
%                                                                       (4)
% B.C.1: At theta>0; a==0;                                      at   R==A
%                                                                       (5) 
% B.C.2: At theta>0; (da/dR+((D_S/D_F)*((r_0)/(r_s-r_0))*a)=0;  at   R = 1  
%                                                                       (6)
% Where 
% theta==((D_ij*t)/(r_o)^2);  R==(r/r_0);  A==(r_A/r_0);
% a==((a_(t))-(a_s))/((a_i)-(a_s));  
% Note that A, is a non-dimensional radius, given by the relation between 
% r_A and r_0, (a_t) is the redness within the cherry’s flesh as 
% time progressed, (a_s) the redness in the cherry's skin in contact 
% with the sweetening solution, and (a_i) the redness within the 
% cherry’s flesh, at the beginning of the experiment.
%
% Then, an infinite series can be expressed, representing the 
% evolution of the dimensionless parameter (a_V)_(theta) in 
% relation to dimensionless time theta, modified to 
% that presented by Maldonado and González Pacheco (2022).
% To determine effective diffusion coefficients and eigenvalues lamba_m, 
% an iterative nonlinear regression method was employed using the least 
% squares method in conjunction with Microsoft Excel Solver® software.
%
% 2 - Problem definition
%-------------------------------------------------------------------------
%Data (for a particular case study):
%       -- Cherry's external radius, r_s, in m : 0.0104755416667
%       -- thick epidermis, (r_s-r_0), in m: 0.0004
%       -- Diffusion coeficient in cherry's flesh in m^2/s: 10^-8-
%       10^-12
%       -- Diffusion coeficient in cherry's skin in m^2/s: 10^-9-
%       10^-14
%       --Constante c=1
%
%Outputs:
%       --Concentration profile of the chemical species, erythrosine or 
%         red gardenia dye, as appropriate, as it advances in the radial 
%         direction within the cherry and on the skin, as well as as a 
%         function of time.
%
% 3 - Process outline (algorithm and program)
%-------------------------------------------------------------------------
%
% The MATLAB pdepe PDE solver solves systems of several parabolic and 
% elliptic PDEs in 1D. The coefficient f(x,t,u,dudx) is a flow term and 
% s(x,t,u,dudx) is a source term.
%
function diffusion1dsphericalcoordinates
function [c,f,s]=pdefun(x,t,u,dudx)
    %Input variables:
    c=1;
    f=dudx;
    s=0;
end

%Condition at the beginning of the domain of t:
    function u0=icfun(x)
        u0=1;
    end
%Boundary conditions in the domain of r:
    function [pL,qL,pR,qR]=bcfun(xL,uL,xR,uR,t)
        pL=uL;% Here it is assumed that the dimensionless concentration is 
        % equal to zero, given that we are at the border of the fruit where
        % the concentration of dye is the same as in the epidermis of the 
        % cherry and equal to that present in the liquid.
        qL=0;% Here "q" on the left is zero due to similarity with the 
        % general equation for the calculation of boundary conditions.
        pR=-3.837654;% This expression is presented since it arises from 
        % proposing a flow balance through the cherry epidermis and 
        % emphasizing the negligible resistance that it offers in the 
        % transport of species.
        qR=1;% Here "q" on the right is iqual 1 due to similarity with the 
        % general equation for the calculation of boundary conditions.

    end
% Space domain 
% The diametral distance of the cherry is approximately 27 mm, and in r=r_a 
% is where the values of dye concentration, therefore, the parameter a* 
% are the lowest at time = 0. As the experiment progresses, more intense 
% coloration values should be observed in r=r_a and even more so in 
% the peripheral area, that is, in the area of the skin of the cherry, 
% since an influx of dye should occur through the cherry flesh and another 
% through its skin, evidencing a double diffusion boundary at r=r_s. This 
% same assumption should be verified by implementing a third Robin-type 
% boundary condition.
x=linspace(0,12,15); % x is equal to r, then, the radial distance es 12 mm.
%Time domain
% Here it is stated that the time of the experiment will be from an initial 
% time, at t=0 to approximately 1 hour, that is, t=60, and about 30 units 
% will be given in the time interval t.
t=linspace(0,60,15);
%Due to it is spherical geometry:
m=2;
%The solution will then be:
sol=pdepe(m,@pdefun,@icfun,@bcfun,x,t);%use of the "sol" command, where
%pdepe returns the solution in a 3-D array sol.
u=sol(:,:,1);%For this problem, u has only one component, so
              % sol is an array, but in general you can extract the kth 
              % component of the solution with the command u = sol(:,:,k). 
              % Then the first component is extracted from the solution 
              % "sol".
%
% 4 - Visualization of the results (application and evaluation)
% ------------------------------------------------------------------------
%
figure(1) 
surf(x,t,u)
axis([0 12 0 60 0 80]);
title ('Erythrosine 357ppm 40°C');
xlabel('Radial Distance (mm)');
ylabel ('Time (min)');
zlabel('Redness, a*');
%***********************************************************************
figure(2)
waterfall(x,t,u)
axis([0 10 0 100 0 65]);
title ('');
xlabel('');
ylabel (')');
zlabel('');
%
figure(3)
mesh(x,t,u)
axis([0 10 0 100 0 65]);
title ('');
xlabel('');
ylabel ('');
zlabel('');
%Contour lines (Azimut=0 Y Elev=90)
figure(4)
contour(x,t,u,50)
axis([0 10 0 100 0 65]);
title ('');
xlabel('');
ylabel ('');
zlabel('');
%
figure(5)
contourf(x,t,u,50)
axis([0 10 0 100 0 65]);
title ('');
xlabel('');
ylabel ('');
zlabel('');
end
