%***********************************************************************
%  This code is written by Chiu-Yen Kao
%  Copyright is reserved  !!



function [u,v] = pp_to_uv(pp)
global xmin xmax ymin ymax Mx My pmeshx pmeshy
global gx gy nx ny nxs nys nxe nye dx dy x y trap_mtr meshx meshy
global s1 s2 nn sigma1 sigma2 r1 r2 alpha1 alpha2 k1 k2 mu rho1 rho2 beta1 beta2 g
global c1 c2

global runTime

u = zeros(nx,ny);
v = zeros(nx,ny);

int_pp = sum(sum(pp.*trap_mtr));
s1 =(c1*exp(-g*runTime)+ sum(sum(pmeshx.*pp.*trap_mtr)))/int_pp;
s2 =(c2*exp(-g*runTime)+ sum(sum(pmeshy.*pp.*trap_mtr)))/int_pp;

%s1 = 0;
%s2 = 0;


for i = nxs : nxe
    for j = nys : nye
        x1 = x(i);
        x2 = y(j);
        u(i,j)= - mu*x1+(alpha1*x1.^nn./(k1^nn+x1.^nn)+sigma1*s1/(rho1+s1))./(1+x2/r2)+beta1;
        v(i,j)= - mu*x2+(alpha2*x2.^nn./(k2^nn+x2.^nn)+sigma2*s2/(rho2+s2))./(1+x1/r1)+beta2;            
    end   
end

% extrapolate u and v
for i = nxs-1 : -1 : 1
    for j = nys : nye
        u(i,j) = 2.*u(i+1,j) - 1.*u(i+2,j);
        v(i,j) = 2.*v(i+1,j) - 1.*v(i+2,j);
    end
end
for i = nxe+1 : nx
    for j = nys : nye
        u(i,j) = 2.*u(i-1,j) - 1.*u(i-2,j) ;
        v(i,j) = 2.*v(i-1,j) - 1.*v(i-2,j);
    end
end
for j = nys-1 : -1 : 1
    for i = 1 : nx
        u(i,j) = 2.*u(i,j+1) - 1.*u(i,j+2) ;
        v(i,j) = 2.*v(i,j+1) - 1.*v(i,j+2) ;
    end
end
for j = nye+1 : ny
    for i = 1 : nx
        u(i,j) = 2.*u(i,j-1) - 1.*u(i,j-2) ;
        v(i,j) = 2.*v(i,j-1) - 1.*v(i,j-2) ;
    end
end