Ocatve script to Plot the electrostatic field lines of a system of charges

field_lines.m

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");

Introduction

This GNU octave¹ 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:

% --------- 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];

% .....

The comments are pretty verbose. But you need to know some basics of
octave. At least enough to define variables and vectors.

Footnotes

1. This script is not compatible with Matlab as it is. ↩