%
% 点電荷のポテンシャル
% 
% MATLAB R2015a
% 2016/7/15 Copyright (c) Takuichi Hirano
%
function charges
    %eps0=8.854e-12;     % 真空の誘電率
    eps0=1.;     % 真空の誘電率

    % グラフィックを描画
    x = linspace(-1.,1.,50);
    y = linspace(-1.,1.,50);
    [X,Y] = meshgrid(x,y);

    n=2;
    q(1)=1.;
    posx(1)=-0.5;
    posy(1)=0.;

    q(2)=-1.;
    posx(2)=0.5;
    posy(2)=0.;

    % R
    R=@(X,Y) sqrt(X.^2+Y.^2);

    % 電位分布の計算
    V=zeros;
    for i=1:n
        V =V +q(i)./(4.*pi.*eps0.*R(X-posx(i),Y-posy(i)));
    end
    [Ex,Ey] = gradient(-V);             % 数値的に勾配を計算

    figure;
    surf(X,Y,V,...
        'EdgeColor','none');              % グラフを描画
    colormap('jet');
    zmin=-0.5; zmax=0.5;
    zlim([zmin zmax]);
    caxis([zmin zmax]);

    % 等高線+ベクトル分布
    % Clip data
    clip=@(E,V,zmin,zmax) E.*((zmin<V)&(V<zmax));
    Ex=clip(Ex,V,zmin,zmax);
    Ey=clip(Ey,V,zmin,zmax);

    figure;
    contour(X,Y,V,linspace(zmin,zmax,11))
    hold on
    quiver(X,Y,Ex,Ey)
    colormap('jet');
    hold off

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


