function pde = Kelloggdata
%%  KELLOGGDATA data of Kellogg Problem
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',0,'exactu',@exactu,'g_D',@exactu,'Du',@Du,'d',@d);

    function K =  d(p)  % Diffusion constant
    idx = (p(:,1).*p(:,2) >0);
    K  = ones(size(p,1),1);
    K(idx) = 161.4476387975881;
    end

    function u =  exactu(p) % exact solution
    gamma = 0.1;
    sigma = -14.9225565104455152;
    rho = pi/4;
    r = sqrt(sum(p.^2,2));
    theta = atan2(p(:,2),p(:,1));
    theta = (theta>= 0).*theta + (theta<0).*(theta+2*pi);
    mu = (theta>= 0 & theta<pi/2).*cos((pi/2-sigma)*gamma).*cos((theta-pi/2+rho)*gamma)...
        +(theta>= pi/2 & theta<pi).*cos(rho*gamma).*cos((theta-pi+sigma)*gamma)...
        +(theta>= pi & theta<1.5*pi).*cos(sigma*gamma).*cos((theta-pi-rho)*gamma)...
        +(theta>= 1.5*pi & theta<2*pi).*cos((pi/2-rho)*gamma).*cos((theta-1.5*pi-sigma)*gamma);
    u =  r.^gamma.*mu;
    end

    function uprime =  Du(p)
    % the gradient of exact solution
    gamma = 0.1;
    sigma = -14.92256510455152;
    rho = pi/4;
    theta = atan2(p(:,2),p(:,1));  % jiaodu
    theta = (theta>= 0).*theta +(theta<0).*(theta+2*pi);
    t = 1+(p(:,2).^2)./(p(:,1).^2);
    r = sqrt(sum(p.^2,2));
    rg = r.^gamma;

    ux1 = (p(:,1)>= 0.0 & p(:,2)>= 0.0).*(rg.*gamma./r.*cos((pi/2-sigma)*gamma)./r.*p(:,1)...
        .*cos((theta-pi/2+rho)*gamma)...
        +rg.*cos((pi/2-sigma)*gamma).*sin((theta-pi/2+rho)*gamma)...
        *gamma.*p(:,2)./(p(:,1).^2)./t);

    uy1 = (p(:,1)>= 0.0 & p(:,2)>= 0.0).*(rg*gamma./r.*cos((pi/2-sigma).*gamma)...
        .*cos((theta-pi/2+rho).*gamma)./r.*p(:,2)...
        -rg.*cos((pi/2-sigma).*gamma).*sin((theta-pi/2+rho)*gamma)...
        *gamma./p(:,1)./t);

    ux2 = (p(:,1)<= 0.0 & p(:,2)>= 0.0).*(r.^(-1.9).*p(:,1)*gamma.*cos(rho.*gamma)...
        .*cos((theta-pi+sigma).*gamma)...
        +rg.*cos(rho*gamma).*sin((theta-pi+sigma).*gamma)*gamma...
        .*p(:,2)./(p(:,1).^2)./t);
    uy2 = (p(:,1)<= 0.0 & p(:,2)>= 0.0).*( r.^(-1.9).*p(:,2)*gamma.*cos(rho.*gamma)...
        .*cos((theta-pi+sigma)*gamma)...
        -rg.*cos(rho*gamma).*sin((theta-pi+sigma).*gamma)*gamma./p(:,1)./t);

    ux3 = (p(:,1)<= 0.0 & p(:,2)<= 0.0).*(r.^(-1.9).*p(:,1).*gamma.*cos(sigma.*gamma)...
        .*cos((theta-pi-rho).*gamma) ...
        +rg.*cos(sigma.*gamma).*sin((theta-pi-rho).*gamma).*gamma...
        .*p(:,2)./(p(:,1).^2)./t);
    uy3 = (p(:,1)<= 0.0 & p(:,2)<= 0.0).*(r.^(-1.9).*p(:,2)*gamma.*cos(sigma*gamma)...
        .*cos((theta-pi-rho).*gamma) ...
        -rg.*cos(sigma*gamma).*sin((theta-pi-rho).*gamma)*gamma./p(:,1)./t);

    ux4 = (p(:,1)>= 0.0& p(:,2)<= 0.0).*(r.^(-1.9).*p(:,1).*gamma.*cos((pi/2-rho)*gamma)...
        .*cos((theta-3*pi/2-sigma)*gamma) ...
         +rg.*cos((pi/2-rho).*gamma).*sin((theta-3*pi/2-sigma)*gamma)...
         *gamma.*p(:,2)./(p(:,1).^2)./t);

    uy4 = (p(:,1)>= 0.0 & p(:,2)<= 0.0).*(r.^(-1.9).*p(:,2)*gamma.*cos((pi/2-rho)*gamma)...
        .*cos((theta-3*pi/2-sigma)*gamma)...
        -rg.*cos((pi/2-rho)*gamma).*sin((theta-3*pi/2-sigma)*gamma)...
        *gamma./p(:,1)./t);     

    uprime(:,1) =  ux1+ux2+ux3+ux4;
    uprime(:,2) =  uy1+uy2+uy3+uy4;
    end
end