surface charge columb law

How Matlab used the columb law to calculate the electric field for points charges , surfaces

does the surface is divided into points charges? What is the approach for this fuction ?

function SurfaceChargeField

%   Program to find V and E of surface (rectangular and circular) 
%   Charge Distribution
%=============================================
% Adjust These Parameters
% Evaluation of electric field at P(x,y,z)
PLoc=[0 0 1]   ;    
  
% The discretization element side d
ElementSize=0.001; 
Eps0=10e-9/(36*pi);
Eps=Eps0;    % Medium is free  space 

%=====================Surface Parameters
% Charge distribution is assumed to be on a rectangle, Disk or Partial Disk
SurfParam.SurfaceShape= 2;  % 1 for rectangle  Patch
                            % 2 for disk, annular disk or partial disk
                           
SurfParam.Rho0=5*10e-9 ;  % Unifrom Rho    
SurfParam.UseRandChargeDist = 0;   % 0 for uninform charge distribution and 1 for random distribution
if SurfParam.SurfaceShape ==1    % Rectangular Patch Shape
   SurfParam.X=[0 1];    % X Min and Max
   SurfParam.Y=[-1 1];    % Y Min and Max
elseif SurfParam.SurfaceShape ==2 
    SurfParam.DiskRad =2 ;      % Outer Radius of disk
    SurfParam.InnerDiskRad=1;   % Inner disk rad (Use 0 for full disk)
   
    SurfParam.Phi=[0 2*pi];     % Start and End Angles For Partial Disk
end
%============================================================
%==========================================================
% Use function rho_s_fun below to adjust the charge density
dx=ElementSize;
dy=ElementSize;   
ds=dx*dy;         % Square Element Areal
format short   % prints only 6 decimal places
%====================Check Parameters======================
if SurfParam.SurfaceShape== 1
    if (SurfParam.X(2) <=SurfParam.X(1) || SurfParam.Y(2) <=SurfParam.Y(1))
        display('Error Check Dimension Values');
        return;
    end
elseif SurfParam.SurfaceShape== 2
    if (SurfParam.DiskRad <=SurfParam.InnerDiskRad)
        display('Error DiskRad is Tool Small');
        return;
    end
    %======= Area Around the Disk
    SurfParam.X=SurfParam.DiskRad*[-1 1];
    SurfParam.Y=SurfParam.DiskRad*[-1 1];
      
end
tic;
                X= SurfParam.X(1) : dx : SurfParam.X(2);
                Y= SurfParam.Y(1) : dy : SurfParam.Y(2);
                [XLoc,YLoc]=meshgrid(X,Y);
                YLoc=flipud(YLoc);   % Make First Row For Maximum Y
                SurfParam.ArraySize=size(XLoc);
                ZLoc=zeros(SurfParam.ArraySize);
                
                % RVec is (r - r')
                RVec=cat(3,(PLoc(1)-XLoc),(PLoc(2)-YLoc), (PLoc(3)-ZLoc));   % | r - rprime |
                RMag=sqrt(RVec(:,:,1).^2+RVec(:,:,2).^2+RVec(:,:,3).^2);
                rho_s=rho_s_fun(XLoc,YLoc,SurfParam);
  %=========== Calculate Potential V
                 VNumerator=ds*rho_s;
                 VDenominator=4*pi*Eps*RMag;  
                 ElemV=VNumerator./VDenominator;
                 VField=sum(sum(ElemV));
                 
  %=========== Calculate E Field
                
                ENumerator=ds*repmat(rho_s,1,1,3).*RVec;    
                EDenominator=4*pi*Eps*repmat(RMag.^3,1,1,3);  % | r - rprime | ^ 3
                ElemField=ENumerator./EDenominator;     
                EField=[sum(sum(ElemField(:,:,1))),sum(sum(ElemField(:,:,2))),sum(sum(ElemField(:,:,3)))];
                    
                
% ===============Display results
figure
imagesc(rho_s);
title('Charge Distribution');
colormap(gray)
colorbar
disp('The computed potential  in V is');
disp(VField);
disp('The computed electric field  in V/m is');
disp(EField);
toc
end
function rho_s=rho_s_fun(XLoc,YLoc,SurfParam)
   SurfaceShape=SurfParam.SurfaceShape;
   Rho0=SurfParam.Rho0;
    
   if  SurfaceShape==1
    rho_s=Rho0*ones(SurfParam.ArraySize) ;  % uniform charge density
  
   %=================================================
   elseif SurfaceShape==2 
  %For Disk set rho_s to be 0 outsie the disk area
    DiskRad=SurfParam.DiskRad;
    InnerDiskRad=SurfParam.InnerDiskRad;
    ElementRho=sqrt(XLoc.^2+YLoc.^2);
    
    PhiMin=SurfParam.Phi(1);
    PhiMax=SurfParam.Phi(2);
    ElemPhi=wrapTo2Pi(atan2(YLoc,XLoc));  % Notice atan2 function gives angle from -pi to pi
    
    
    rho_s=zeros(size(XLoc));
    
    DiskElem=(ElementRho >= InnerDiskRad & ElementRho <= DiskRad  & ...
        ElemPhi >= PhiMin & ElemPhi <= PhiMax);
      % We choose elements inside the disk
       rho_s=Rho0*DiskElem ;                              %  Uniform Charge
 %====================================================== 
        % Example of Nonuniform charge density 
       
        % rho_s=Rho0*ElementRho.^2.*DiskElem;           
%===========================================
%  Rand Charge Distribution ( TO show elements)
if SurfParam.UseRandChargeDist==1
    rho_s=rho_s.*randn(size(rho_s));
end
%============================
      
  end
end

Answer 1

[Simran]

Hi @abdullah,

As I understand, you want to know how, the Coulomb’s law, which is used to calculate electric field of point charges, is used in the provided MATLAB function, to calculate electric field for surface charged distribution.

The main approach you mentioned is that the surface is divided into many small elements, each of which can be treated as a point charge with a small amount of charge. By calculating the contribution from each small element and summing them up, we can approximate the total electric field and potential due to the entire surface.

We can see the same in these parts of the code:

1.) Dividing the Surface into Small Parts:

dx = ElementSize;

dy = ElementSize;

ds = dx * dy;

[XLoc, YLoc] = meshgrid(X, Y);

Here the surface is divided into small rectangular elements of size dx/dy, with each element considered a small "point charge".

You can refer to the below documentation for the more information :

https://www.mathworks.com/help/matlab/ref/meshgrid.html

2.) Charge on Each Element:

rho_s = rho_s_fun(XLoc, YLoc, SurfParam);

Here the function ”rho_s_fun” assigns a charge density ”rho_s” to each element based on its position and the specified distribution.

3.) Application of Coulomb’s Law to Each Element:

3.1) Potential Calculation:

VNumerator = ds * rho_s;

VDenominator = 4 * pi * Eps * RMag;

ElemV = VNumerator ./ VDenominator;

VField = sum(sum(ElemV));

Here for each element, we calculate its contribution to the potential at point “P“ using: [V = ρs × ds/4πϵ0×R ] where “R” is the distance from the element to point “P”. The sum of all these contributions gives us the total potential field “Vfield”.

3.2) Electric Field Calculation:

ENumerator = ds * repmat(rho_s, 1, 1, 3) .* RVec;

EDenominator = 4 * pi * Eps * repmat(RMag.^3, 1, 1, 3);

ElemField = ENumerator ./ EDenominator;

EField = [sum(sum(ElemField(:,:,1))), sum(sum(ElemField(:,:,2))), sum(sum(ElemField(:,:,3)))];

Here for each element, we calculate its contribution to the electric field at point “P” using: [E = ρs × ds × R/4πϵ0×R^3] where “R” is the vector from the element to point “P”. Sum of all these vector contributions get the total electric field, “Efield”.

The concepts used in the code are based on “standard numerical methods”, which allows us to approximate the effects of a surface charge by summing the contributions from each discrete element and the “electromagnetic theory”.

Below are the links attached that maybe relevant to the discussion:

For more info about numerical integration- https://www.mathworks.com/help/matlab/numerical-integration-and-differential-equations.html

For more info about meshgrid- https://www.mathworks.com/help/matlab/ref/meshgrid.html