%
% 一様電界中に置かれた誘電体球
% 
% MATLAB R2015a
% 2016/7/24 Copyright (c) Takuichi Hirano
%
function di_sphere_in_efield
    global eps1 a e0 e1 p

    %eps0=8.854e-12;     % 真空の誘電率
    eps0=1.;            % 真空の誘電率

    %---- 入力パラメータ
    eps1=1.0*eps0;      % 背景領域の誘電率
    eps2=4.0*eps0;      % 球の誘電率

    e0=1.0;
    a=1.0;
    %----

    e1=((3.*eps1)./(2.*eps1+eps2)).*e0;
    p=4.*pi.*a.^3.*((eps1.*(eps2-eps1))./(2.*eps1+eps2)).*e0;

    % グラフィックを描画
    xmin=-4*a;
    xmax= 4*a;
    ymin=-4*a;
    ymax= 4*a;
    x = linspace(xmin,xmax,200);
    y = linspace(ymin,ymax,200);
    [X,Y] = meshgrid(x,y);

    % R
    R=@(x,y) sqrt(x.^2+y.^2);

    % 電位分布の計算
    %V=1.*heaviside(R(X,Y));
    V=voltage(X,Y);
    [Ex,Ey] = gradient(-V);             % 数値的に勾配を計算

    figure;
    contour(X,Y,V,[-4*a:a/5:4*a]);             % グラフを描画
%    surf(X,Y,V);              % グラフを描画
    colormap('jet');
%    zmin=-0.5; zmax=0.5;
    axis tight;
    zlim([-4*a,4*a]);
    caxis([-4*a,4*a]);

    % 等高線+ベクトル分布
    hold on
    streamslice(X,Y,Ex,Ey);
    hold off

    % 等高線+力線
%    figure;
%    contour(X,Y,V,linspace(zmin,zmax,11))
%    hold on
%    streamslice(X,Y,Ex,Ey)
%    colormap('jet');
%    hold off
end

function v = voltage(x,y)
    global eps1 a e0 e1 p

    % R
    R=@(x,y) sqrt(x.^2+y.^2);

    v=(-e0.*x+(p.*(x./R(x,y)))./(4.*pi.*eps1.*(R(x,y).^2))).*heaviside(R(x,y)-a) ... % 球の外側
     +(-e0.*x+e1.*x).*heaviside(-(R(x,y)-a));
end

