Last active
May 20, 2023 21:01
-
-
Save SidhBhat/d5d8df26228bf3e8659dfbcd786d82af to your computer and use it in GitHub Desktop.
Ocatve script to Plot the electrostatic field lines of a system of charges
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
clear all; | |
clf; | |
global Q; | |
global K; | |
global plot_len; | |
global line_density; | |
% --------- User Configuration --------- % | |
%% Charges are defined as structures 'Q' containing two fields 'q' and 'x' | |
%% specifying the charge and position position is given as a two dimensional | |
%% row vector. Multiple charges can be defined by creation an array of these | |
%% structures. | |
%% Define Charges | |
%% The default setup is an electric dipole. | |
%% The separation between charges | |
sep = 0.1; | |
Q(1).q = 1; | |
Q(1).x = [sep,0]; | |
Q(2).q = -1; | |
Q(2).x = [-sep,0]; | |
% % Uncomment the below lines for additional test charges | |
% Q(3).q = 3; | |
% Q(3).x = [-sep,-sep]; | |
% | |
% Q(4).q = -1; | |
% Q(4).x = [-sep,sep]; | |
% | |
% Q(5).q = 5; | |
% Q(5).x = [0.5,-0.6]; | |
% Coulombs constant, It really does not matter what you set here, as long | |
% as it is non zero | |
K = 1; | |
% This defines the plot limits. Specifically the box if a square box with | |
% x,y limits between -plot_len to plot_len. | |
plot_len = 1; | |
% The number of lines originating from each charge | |
line_density = 20; | |
% The size of the smallest charge | |
marker_size = 20; | |
% --------- End Configuration --------- % | |
%% Everything below this line is not meant to be normally edited by | |
%% the user unless you know what you are doing. | |
global number = uint8(length(Q)); % Number of charges | |
global theta = 2*pi / line_density; % The angle between lines originating from a charge | |
% A critical variable, influencing the resolution of the field lines | |
global gaps = 0.01; | |
% normalise a vector. | |
function v = normalise(vec) | |
v = vec/norm(vec); | |
endfunction | |
% Initialise the reference vector for each Q | |
% and the exclude line list. | |
for i = 1:number; | |
Q(i)._exclude_ = 0; | |
if (norm(Q(i).x)) | |
Q(i)._r_ = normalise(Q(i).x); | |
else | |
Q(i)._r_ = [1,0]; | |
endif | |
endfor | |
% calculate field at point r | |
function E = e_field(r) | |
global Q; | |
global K; | |
global number; | |
E = [0,0]; | |
for i = 1:number; | |
E += K*Q(i).q*(r - Q(i).x)/norm(r - Q(i).x)^3; | |
endfor | |
endfunction | |
% angle between two vectors | |
function ang = vector_angle(r1,r2) | |
k = cross([r1,0],[r2,0])(3); | |
ang = acos(dot(r1,r2)/norm(r1)/norm(r2)); | |
if(k < 0) | |
ang = 2*pi - ang; | |
endif | |
endfunction | |
% This function returns true if the point R is close to a charge | |
% this is defined as 50% of the 'gaps' variable | |
function val = charge_proxi(R) | |
global Q; | |
global gaps; | |
global number; | |
global theta; | |
global r; | |
val = false; | |
for i = 1:number; | |
if(norm(R - Q(i).x) <= 0.5*gaps) | |
val = true; | |
ang = round(vector_angle(Q(i)._r_, R - Q(i).x)/theta); | |
Q(i)._exclude_(end + 1) = uint8(ang); | |
return; | |
endif | |
endfor | |
endfunction | |
% this function returns true if a line indexed by n should be plotted | |
% according to Q(q)'s exclude list. | |
function val = plot_line(q,n) | |
global Q; | |
val = !any(Q(q)._exclude_ == n); | |
endfunction | |
% returns true if r is within the boundaries of the plot window | |
function val = within_bounds(r) | |
global plot_len; | |
val = true; | |
if(abs(r(1)) > plot_len || abs(r(2)) > plot_len) | |
val = false; | |
endif | |
endfunction | |
% returns the rotation matrix for a counterclockwise rotation by n*theta. | |
function rot = R(n) | |
global line_density; | |
global theta; | |
rot = [ cos(n*theta), -sin(n*theta) ; sin(n*theta), cos(n*theta) ]; | |
endfunction | |
% The integration distance for calculation a field line. | |
dels = gaps/10; | |
% Start plot | |
hold("on"); | |
for j = 1:number; | |
for i = 1:line_density; | |
if(plot_line(j,i)) | |
X.x = 0; X.y = 0; | |
X.x(end) = (x = Q(j).x + 0.501*gaps*Q(j)._r_*(R(i)'))(1); | |
X.y(end) = x(2); | |
E = normalise(e_field(x)); | |
if(Q(j).q < 0) | |
E = -E; | |
endif | |
while(!charge_proxi(x) && within_bounds(x)) | |
X.x(end+1) = (x += E*dels)(1); | |
X.y(end+1) = x(2); | |
E = normalise(e_field(x)); | |
if(Q(j).q < 0) | |
E = -E; | |
endif | |
endwhile | |
plot(X.x,X.y,"-k","markersize", 5); | |
clear X E x; | |
endif | |
endfor | |
endfor | |
% smallest charge. | |
smallest_Q = abs(Q(1).q); | |
for i = 1:number; | |
if(smallest_Q > abs(Q(i).q)) | |
smallest_Q = abs(Q(i).q); | |
endif | |
endfor | |
% Plot markers for the charges | |
for i = 1:number; | |
if(Q(i).q > 0) | |
plot(Q(i).x(1),Q(i).x(2),".r","markersize",abs(Q(i).q)/smallest_Q*marker_size); | |
else | |
plot(Q(i).x(1),Q(i).x(2),".b","markersize", abs(Q(i).q)/smallest_Q*marker_size); | |
endif | |
endfor | |
% define limits and axis. | |
plot_len += 0.1; | |
axis("off","equal",[-plot_len,plot_len,-plot_len,plot_len]); | |
grid("off"); | |
xlabel("x"); | |
ylabel("y"); |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Introduction
This GNU octave1 script plots the electrostatic field lines
of a system of charges that you define.
Usage
Using this script is simple. Just edit the section titled
User Configuration
:The comments are pretty verbose. But you need to know some basics of
octave. At least enough to define variables and vectors.
Footnotes
This script is not compatible with Matlab as it is. ↩