%% Chapter 10. Physics
%% Motion in Central Force Field

syms x y z real
syms f(r) phi(r)
rvec = [x, y, z]
r = sqrt(dot(rvec, rvec))
curl(f(r)*rvec)

gradient(phi(r))

simplify(divergence((1/r^3)*rvec))

%
%% A Newtonian Vector Field

figure
a = 1; b=1;
[X, Y, Z] = meshgrid(-1.25:0.5:1.25, -1.25:0.5:1.25, -1.25:0.5:1.25);
quiver3(X, Y, Z, -a.*(X-1)./((X-1).^2+Y.^2+Z.^2).^(1./2)...
    -b.*(X+1)./((X+1).^2+Y.^2+Z.^2).^(1./2),...
    -a.*Y./((X-1).^2+Y.^2+Z.^2).^(1./2)-b.*Y./((X+1).^2+Y.^2+Z.^2).^(1./2),...
    -a.*Z./((X-1).^2+Y.^2+Z.^2).^(1./2)-b.*Z./((X+1).^2+Y.^2+Z.^2).^(1./2))
axis equal

%
%% An Elecromagnetic Field

b = -1
figure;
[X, Y, Z] = meshgrid(-1.25:0.5:1.25, -1.25:0.5:1.25, -1.25:0.5:1.25);
quiver3(X, Y, Z, -a.*(X-1)./((X-1).^2+Y.^2+Z.^2).^(1./2)...
    -b.*(X+1)./((X+1).^2+Y.^2+Z.^2).^(1./2),...
    -a.*Y./((X-1).^2+Y.^2+Z.^2).^(1./2)-b.*Y./((X+1).^2+Y.^2+Z.^2).^(1./2),...
    -a.*Z./((X-1).^2+Y.^2+Z.^2).^(1./2)-b.*Z./((X+1).^2+Y.^2+Z.^2).^(1./2))
axis equal