%
%   Book: MATLAB-Based Electromagnetics (Pearson Prentice Hall)
%   Author: Branislav M. Notaros
%   Instructor Resources
%   (c) 2011
%
%   This MATLAB code or any part of it may be used only for educational purposes 
%   associated with the book
%
%
%


%   Charged hemisphere -- numerical integration


clear all;
close all;
ros = input('Enter a uniform surface charge density in Coulombs per square meter: ');
a = input('Enter the radius of the hemisphere : ');
z = input('Enter the point position at z axis outside the hemisphere: ');
EPS0 = 8.854e-12; %Permittivity of Air

%Due to symmetry there is only z component of electric field
E = 0; 
d_theta = pi/1000;

for theta = 0:d_theta:(pi/2)
    ds = 2*pi*a^2*sin(theta)*d_theta;
    dE = (ros/(4*pi*EPS0)*ds*(z - a*cos(theta)))/...
        (sqrt(z^2 + a^2 - 2*a*z*cos(theta)))^3;
    E = E + dE;
end
fprintf('Numerical result for the electric field at the point at z axis is %f V/m.\n', E);
% analytical formula
Ea = (ros*a^2)/(2*EPS0*z^2)*(a/(sqrt(z^2 + a^2)) + (z - a)/abs(z - a));
fprintf('Analytical result for the electric field at the same point is %f V/m.\n', Ea);
error = abs(Ea - E)*100/Ea;
fprintf(' Relative error in percent is %f\n', error);