% PARAMETERS
% fi - porosity
% alfa - coefficient of linear solid thermal expansion
% E - Young modulus
% nu - Poisson ratio
% gamma - specific heat
% kapa - thermal conductivity
% k - mobility
% Kf - fluid bulk modulus
% b - Biot coefficient
% alfaf - coefficient of linear fluid thermal expansion
% ro - mass density
% P0 - applied pressure
% T0 - initial temperature
% R0 - sphere radius
% Pf - initial fluid pressure
% 
fi = 0.2;
alfa = 0.000012;
E = 10000;
nu = 0.25;
K0 = E/(3*(1-2*nu));
beta = 3*alfa*K0;
lambda = E*nu/((1+nu)*(1-2*nu));
mu = E/(2*(1+nu));
gamma = 1000;
kapa = 0.0000025;
k = 0.0001;
Kf = 1/0.00046;
M = Kf/fi;
b = 1;
alfaf = 0.00005;
alfam = alfa*(b-fi)+fi*alfaf;
ro = 0.0026;
P0 = 1000;
T0 = 20;
R0 = 4;
Pf = 0;
%
%COEFFICIENTS
% Ct - thermal diffusion coefficient
% Cf - fluid diffusion coefficient
% J - matrix J
% S - J^(-1)
% AA - matrix A, BB - matrix B
% AAp - matrix Ap, BB - matrix Bp  for pressure responce
% DD - matrix D, FF - F, GG - G, HH - H, KK - F'
% delta - instantaneous responce for the temperature
% deltap - instantaneous responce for the pressure
%
Ct = (kapa/(ro*gamma))*(lambda+2*mu)/((lambda+2*mu)...
   +T0*beta*beta/(ro*gamma));
Cf = k*M*(lambda+2*mu)/(lambda+2*mu+b*b*M);
d = b*beta/(lambda+2*mu)-3*alfam;
A0 = 1/(Ct*Cf)-d*d*T0/(k*kapa);
a11 = 1/(Cf*A0);
a12 = -d*T0/(kapa*A0);
a21 = -d/(k*A0);
a22 = 1/(Ct*A0);
J = [a11, a12; a21, a22];
S = [1/Ct, d*T0/kapa; d/k, 1/Cf];
A1 = T0*beta/kapa-3*alfam*T0*(lambda+2*mu)/(kapa*b);
A2 = ro*gamma/kapa+3*alfam*T0*beta/(kapa*b);
A3 = 3*alfam*T0*(lambda+2*mu)/(kapa*b);
B1 = b/k+(lambda+2*mu)/(b*M*k);
B2 = 3*alfam/k+beta/(b*M*k);
B3 = (lambda+2*mu)/(b*M*k);
A1p = T0*beta/kapa+ro*gamma*(lambda+2*mu)/(kapa*beta);
A2p = -ro*gamma*b/(kapa*beta)-3*alfam*T0/kapa;
A3p = -ro*gamma*(lambda+2*mu)/(kapa*beta);
B1p = b/k-3*alfam*(lambda+2*mu)/(beta*k);
B2p = -3*alfam*b/(beta*k)-1/(M*k);
B3p = -3*alfam*(lambda+2*mu)/(beta*k);
AA = [A1, A2; B1, -B2];
BB = [A3, -A2; -B3, B2];
AAp = [A1p, A2p; B1p, -B2p];
BBp = [A3p, -A2p; -B3p, B2p];
I = [1, 0; 0, 1];
DD = [4*mu/(R0*R0*R0*(lambda+2*mu)), 0; 0, 0];
FF = [-P0/(lambda+2*mu); T0];
GG = [T0*beta/kapa, 0; b/k, 0];
HH = GG/(I+DD/AA*BB*R0*R0*R0/3)*DD/AA;
KK = GG/(I+DD/AA*BB*R0*R0*R0/3)*FF;
HHp = GG/(I+DD/AAp*BBp*R0*R0*R0/3)*DD/AAp;
KKp = GG/(I+DD/AAp*BBp*R0*R0*R0/3)*FF;
delta = -3*P0*T0*(b*alfam*M+alfa*K0)/(ro*gamma*(K0+b*b*M)+9*K0*T0*...
        (alfa*alfa*K0+alfam*M*(2*b*alfa-alfam)));
deltap = -P0*M*(9*alfa*alfam*K0*T0+b*ro*gamma)/(ro*gamma*(K0+b*b*M)+...
         9*K0*T0*(alfa*alfa*K0+alfam*M*(2*b*alfa-alfam)));
%    
% ROOTS FINDING
% X3 - roots for the temperature responce
% X3p - roots for the presure responce
x11 = pi/R0;
x12 = 0;
x21 = 0;
x22 = pi/R0;
X0 = [x11, x12; x21, x22];
X0p = [x11, x12; x21, x22];
n0 = 100;
X1 = zeros(2,2);
X3 = zeros(2,2);
X1p = zeros(2,2);
X3p = zeros(2,2);
for i = 1:n0
    X1 = fsolve(@(X)(funm(R0*X,@sin)/funm(R0*X,@cos)-(HH-X*X)\(HH*R0*X)),X0);
    X1p = fsolve(@(X)(funm(R0*X,@sin)/funm(R0*X,@cos)-(HHp-X*X)\(HHp*R0*X)),X0);
    X3 = cat (3, X3, X1);
    X3p = cat (3, X3p, X1p);
    X0 = X0+[pi/R0,0;0,pi/R0];
    X0p = X0p+[pi/R0,0;0,pi/R0];
end
%      
% SOLUTION
% vv - temperature distribution
% pp - pressure distribution
% tfn - non-dimensional time for temperature
% tf0n - non-dimensional time for pressure
t = 1;
nn = 50;
d_t = t/nn;
tr1=zeros(2,2);
tr1p=zeros(2,2);
tr=zeros(2,2);
trp=zeros(2,2);
tf1=zeros(1,nn);
tf2=zeros(1,nn);
tf02=zeros(1,nn);
tfn=zeros(1,nn);
tf0n=zeros(1,nn);
tf3=zeros(1,nn);
tf03=zeros(1,nn);
vv=zeros(1,nn);
pp=zeros(1,nn);
for j=1:nn
    tf1(j) = d_t*(j);
    tf2(j) = R0*R0*tf1(j)/(Ct*3);
    tf02(j) = R0*R0*tf1(j)/(Cf*3);
    tr=zeros(2,2);
    trp=zeros(2,2);
    % Summation of the equation (46) (tr - for temperature, trp - for pressure)
    for i = 1:n0
        sn = X3(:,:,i+1);
        snp = X3p(:,:,i+1);
        tr = tr+expm(-J*(sn*sn)*tf2(j))*((HH-sn*sn)/...
            funm(sn*R0,@cos)-HH)/(sn*sn-3*HH+R0*R0*HH*HH); 
        trp = trp+expm(-J*(snp*snp)*tf02(j))*((HHp-snp*snp)/...
            funm(snp*R0,@cos)-HHp)/(snp*snp-3*HHp+R0*R0*HHp*HHp);
    end
    tr1 = cat (3, tr1, tr);
    tr1p = cat (3, tr1p, trp);
    tr2 = AA\tr1(:,:,j)*KK*2;
    tr2p = AAp\tr1p(:,:,j)*KKp*2;
    AK = AA\KK;
    AKp = AAp\KKp;
    vv(j)=-tr2(2)+T0;
    pp(j)=-tr2p(2)/P0;
    if j==1
        tf3(j)=0;
        tf03(j)=0;
    else
        tf3(j)=tf2(j-1);
        tf03(j)=tf02(j-1);
    end
    tf0n = 3*Cf*tf03/(R0*R0);
    tfn = 3*Ct*tf3/(R0*R0);
end
%plot(tfn,vv);
plot(tf0n,pp);
        
    