File: ef230-2023-08/www/examples/coulomb.m Download
%! Electric forces in a system of particles - Coulomb's Law
% University of Tennessee : EF 230 Fall, 2009 : Will Schleter

function coulomb
clear all; close all; clc;

% Set values
% q(i) - charge of particle i
% qs - charge scale (10^qs) C
% p(i) - position [x y z] of particle i
% ps - position scale (10^ps) m
[q qs p ps] = getchargedata();

% init
k=9; ks=9; % constant and scale for Coulomb's Law
fs = (ks+qs+qs)-(ps+ps); % scaling for the result
pd = [];
figure(1);
hold on;
axis equal
title('Coulomb''s Law Example');
vs=.5;

% Loop through all particles
for a=1:length(q);
    f=[0 0 0]; % init force on particle a
    pa=p(a,:); % position of a as [x y z]
    for b=1:length(q);
        if a==b, continue, end; % skip if same particle
        pb=p(b,:); % position of b as [x y z]
        r=norm(pb-pa); % distance between a and b
        u=(pb-pa)./r; % unit vector from a to b
        if sign(q(a)) == sign(q(b))
            u=-u; % flip direction if repelling
        end
        fab=k.*abs(q(a).*q(b))./(r.^2); % magnitude of the force between a and b
        f=f+fab.*u; % add b's force on a to total force on a
    end
    F(a,:)=f; % save results in F
    fprintf('Force on %u is (%g %g %g)x10^%d N (%.3g @ %.3g deg)\n',a,f,fs,norm(f),atan2(f(2),f(1))*180/pi);
    pd(end+1,:)=[pa f]; % save plot data
    text(pa(1)+.02,pa(2),pa(3),sprintf('%u(%d)',a,q(a))); % show charge
    plot3(pa(1),pa(2),pa(3),'*');
end
quiver3(pd(:,1),pd(:,2),pd(:,3),pd(:,4),pd(:,5),pd(:,6),vs);

function [q qs p ps] = getchargedata()
% user interface for entering data or selecting a data set
choices = {'User Entered','21-2','21-3','21-4'};
c=menu('Examples',choices);
switch c
    case 1
        prompt={'Number of particles','Charge magnitude (C)','Position magnitude (m)'};
        ui=inputdlg(prompt);
        n=str2num(ui{1});
        qs=str2num(ui{2});
        ps=str2num(ui{3});
        for i=1:n, prompt{i}=num2str(i); end;
        title='q,x,y,z for each particle';
        ui=inputdlg(prompt,title);
        for i=1:n
            tmp = sscanf(ui{i},'%f,%f,%f,%f');
            q(i)=tmp(1);
            p(i,:)=tmp(2:4);
        end
    case 2
        q(1)=25;
        q(2)=-75;
        qs=-9;
        p(1,:)=[0 0 0];
        p(2,:)=[3 0 0];
        ps=-2;
    case 3
        q(1)=1;
        q(2)=-3;
        q(3)=5;
        qs=-9;
        p(1,:)=[2 0 0];
        p(2,:)=[4 0 0];
        p(3,:)=[0 0 0];
        ps=-2;
    case 4
        q(1)=2;
        q(2)=2;
        q(3)=4;

        p(1,:)=[0 .3 0];
        p(2,:)=[0 -.3 0];
        p(3,:)=[.4 0 0];
        qs=-6;
        ps=0;
    otherwise
        error('No choice made');
end