% Analytical solutions for the incompressible Stokes flow 
% Function return eta,rho,vx,vy,pr
function[eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,x,y,z,gx,gy,gz,gridtest)

% Check coordinates
if(x<0)
    x=0;
end
if(x>xsize)
    x=xsize;
end
if(y<0)
    y=0;
end
if(y>ysize)
    y=ysize;
end
if(z<0)
    z=0;
end
if(z>zsize)
    z=zsize;
end

% gridtest=1,2,3 model parameters
nx    = 1;
nz    = 1;
km    = 1.6*pi;
kn    = 3  *pi;
sigma = 1;
B     = 6.9;
mc    = 1e3;
mb    = 1e6;

% gridtest=4,5 model parameters
eta1=1; % Viscosity at x=0,y=0,z=0 Pa s
eta2=5; % Viscosity at x=0,y=ysize,z=0 Pa s
eta3=5; % Viscosity at x=xsize,y=0,z=0 Pa s
eta4=5; % Viscosity at x=0,y=0,z=zsize Pa s
% Density rho=beta1*eta+beta2;
beta1=1;
beta2=3000;
% Integration constants 2D
cc1=100; %Velocity
cc2=1; %Vx
cc3=2; %Vy
cc4=0; % Pressure constant
% Integration constants 3D
% c1
c1x=100; %Vx
c1y=100; %Vy
c1z=100; %Vz
% c2
c2x=1; %Vx
c2y=1; %Vy
c2z=1; %Vy
cP=0; % Pressure constant
% Gravity
Gx=gx;
Gy=gy;
Gz=gz;

% Exponential viscosity 3D
if(gridtest==7)
    % Koefficients for the exponential viscosity equation eta=c*exp(a*x+b*y)
    C=eta1;
    a=(log(eta3)-log(eta1))/xsize;
    b=(log(eta2)-log(eta1))/ysize;
    c=(log(eta4)-log(eta1))/zsize;
    % Limitation a*c1x+b*c1y+c*c1z=0
    if(a~=0)
        c1x=-(b*c1y+c*c1z)/a;
    elseif(b~=0)
        c1y=-(a*c1x+c*c1z)/b;
    elseif(c~=0)
        c1z=-(a*c1x+b*c1y)/c;
    end
    % Compute results
    eta=C*exp(a*x+b*y+c*z);
    rho=beta1*eta+beta2;
    a1=beta1*(Gy*a*b+Gz*a*c-Gx*(b*b+c*c))/(a*a+b*b+c*c)^2;
    a2=beta2*(Gy*a*b+Gz*a*c-Gx*(b*b+c*c))/ (a*a+b*b+c*c)^2;
    b1=beta1*(Gx*a*b+Gz*b*c-Gy*(a*a+c*c)) / (a*a+b*b+c*c)^2;
    b2=beta2*(Gx*a*b+Gz*b*c-Gy*(a*a+c*c)) / (a*a+b*b+c*c)^2;
    d1=beta1*(Gx*a*c+Gy*b*c-Gz*(a*a+b*b)) / (a*a+b*b+c*c)^2;
    d2=beta2*(Gx*a*c+Gy*b*c-Gz*(a*a+b*b)) / (a*a+b*b+c*c)^2;
    p1=beta1* (Gx*a+Gy*b+Gz*c) / (a*a+b*b+c*c);
    p2=beta2* (Gx*a+Gy*b+Gz*c) / (a*a+b*b+c*c);
    vx=a1*log(eta)+(a1-a2-c1x)/eta-a2*log(eta)/eta+c2x;
    vy=b1*log(eta)+(b1-b2-c1y)/eta-b2*log(eta)/eta+c2y;
    vz=d1*log(eta)+(d1-d2-c1z)/eta-d2*log(eta)/eta+c2z;
    pr=p1*eta+p2*log(eta)+cP;

% Linear viscosity 3D
elseif(gridtest==6)
    % Koefficients for the linear viscosity equation eta=a*x+b*y+c
    e=eta1;
    a=(eta3-eta1)/xsize;
    b=(eta2-eta1)/ysize;
    c=(eta4-eta1)/zsize;
    eta=a*x+b*y+c*z+e;
    rho=beta1*eta+beta2;
    % Limitation a*c1x+b*c1y+c*c1z=0
    if(a~=0)
        c1x=-(b*c1y+c*c1z)/a;
    elseif(b~=0)
        c1y=-(a*c1x+c*c1z)/b;
    elseif(c~=0)
        c1z=-(a*c1x+b*c1y)/c;
    end
    a1=beta1*(Gy*a*b+Gz*a*c-Gx*(b*b+c*c))/(a*a+b*b+c*c)^2;
    a2=beta2*(Gy*a*b+Gz*a*c-Gx*(b*b+c*c))/ (a*a+b*b+c*c)^2;
    b1=beta1*(Gx*a*b+Gz*b*c-Gy*(a*a+c*c)) / (a*a+b*b+c*c)^2;
    b2=beta2*(Gx*a*b+Gz*b*c-Gy*(a*a+c*c)) / (a*a+b*b+c*c)^2;
    d1=beta1*(Gx*a*c+Gy*b*c-Gz*(a*a+b*b)) / (a*a+b*b+c*c)^2;
    d2=beta2*(Gx*a*c+Gy*b*c-Gz*(a*a+b*b)) / (a*a+b*b+c*c)^2;
    p1=beta1* (Gx*a+Gy*b+Gz*c) / (a*a+b*b+c*c);
    p2=beta2* (Gx*a+Gy*b+Gz*c) / (a*a+b*b+c*c);
    vx=-(a1/2+a2-c1x)*log(eta)+1/4*a1*eta*eta+a2*eta-1/4*a1-a2+c2x;
    vy=-(b1/2+b2-c1y)*log(eta)+1/4*b1*eta*eta+b2*eta-1/4*b1-b2+c2y;
    vz=-(d1/2+d2-c1z)*log(eta)+1/4*d1*eta*eta+d2*eta-1/4*d1-d2+c2z;
    pr=1/2*p1*eta*eta+p2*eta-1/2*p1-p2+cP;

% Exponential viscosity 2D
elseif(gridtest==5)
    % Koefficients for the exponential viscosity equation eta=c*exp(a*x+b*y)
    c=eta1;
    a=(log(eta3)-log(eta1))/xsize;
    b=(log(eta2)-log(eta1))/ysize;
    % Gravity related constants
    rho=beta1*c*exp(a*x+b*y)+beta2;
    A=rho*(a*gy-b*gx)/(a^2+b^2)^2;
    a1=beta1*(a*gy-b*gx)/(a^2+b^2)^2;
    a2=beta2*(a*gy-b*gx)/(a^2+b^2)^2;
    A1=rho*(b*gy+a*gx)/(a^2+b^2);
    b1=beta1*(b*gy+a*gx)/(a^2+b^2);
    b2=beta2*(b*gy+a*gx)/(a^2+b^2);
    % Compute results
    eta=c*exp(a*x+b*y);
    vx=b*a1*log(eta)+b*(a1-a2-cc1)/eta-b*a2*log(eta)/eta +cc2;
    vy=-a*a1*log(eta)-a*(a1-a2-cc1)/eta+a*a2*log(eta)/eta +cc3;
    vz=0;
    pr=b1*eta+b2*log(eta)+cc4;
    
% Linear viscosity 2D
elseif(gridtest==4)
    % Koefficients for the linear viscosity equation eta=a*x+b*y+c
    c=eta1;
    a=(eta3-eta1)/xsize;
    b=(eta2-eta1)/ysize;
    % Gravity related constants
    rho=beta1*(a*x+b*y+c)+beta2;
    A=rho*(a*gy-b*gx)/(a^2+b^2)^2;
    a1=beta1*(a*gy-b*gx)/(a^2+b^2)^2;
    a2=beta2*(a*gy-b*gx)/(a^2+b^2)^2;
    A1=rho*(b*gy+a*gx)/(a^2+b^2);
    b1=beta1*(b*gy+a*gx)/(a^2+b^2);
    b2=beta2*(b*gy+a*gx)/(a^2+b^2);
    % Compute results
    eta=a*x+b*y+c;
    vx=b*a1*eta*eta/4+b*a2*eta+b*(cc1-a2-a1/2)*log(eta)+cc2;
    vy=-a*a1*eta*eta/4-a*a2*eta-a*(cc1-a2-a1/2)*log(eta)+cc3;
    vz=0;
    pr=b1*eta*eta/2+b2*eta+cc4;
    
% Inclusion
elseif (gridtest==1)
    x=x-xsize/2;
    y=y-ysize/2;
    % Matrix
    eta=1;
    rho=1;
    % Inclusion
    cdist=(x^2+y^2)^0.5;
    if(cdist<0.5*xsize/2)
        eta=mc;
    end
    g.xmin=-xsize/2;
    g.xmax= xsize/2;
    g.ymin=-ysize/2;
    g.ymax= ysize/2;
    sol      = eval_anal_Dani( x, y, g );
    vx=sol.vx;
    vy=sol.vz;
    vz=0;
    pr=sol.P;

% solcx
elseif(gridtest==2)
    rho = cos( nx*pi*x )*sin( nz*pi*y );
    eta = 1;
    if (x >= 0.5*xsize)
        eta = mb;
    end
    sol      = eval_sol_cx( x, y );
    vx=sol.vx;
    vy=sol.vz;
    vz=0;
    pr=sol.P;

% solkz
elseif (gridtest==3)
    rho = -sigma*sin(km*y)*cos(kn*x);
    eta = exp(2*B*y);
    sol      = eval_sol_kz( x, y );
    vx=sol.vx;
    vy=sol.vz;
    vz=0;
    pr=sol.P;

end
