%poles and corresponding equipotential, force curves
clf;close all;
[a,b,c]=sphere(10);
s = 'mis-aligned magnets:equipotential surfaces and B-field';
m1 = 5; m2 = 1; mu1 = 5; mu2 = 5; gee = 32.;
xn(1) = 0; xs(1)= 0; xn(2)= 0+sin(pi/60); xs(2)= 0;
yn(1) = -2; ys(1)= -1; yn(2)= 1+cos(pi/60); ys(2)= 1;
axis([-3 3 -3 3]); hold on
xgrid = -3:.02:3; ygrid = -3:.02:3;
for i=1:2
contour(.1*a+xn(i),.1*b+yn(i),.1*c)
contour(.1*a+xs(i),.1*b+ys(i),.1*c)
end
% forces on each magnet
% first compute pairwise distances
rnn=(((xn(1)-xn(2))^2+(yn(1)-yn(2))^2))^(-3/2);
rns=(((xn(1)-xs(2))^2+(yn(1)-ys(2))^2))^(-3/2);
rsn=(((xs(1)-xn(2))^2+(ys(1)-yn(2))^2))^(-3/2);
rss=(((xs(1)-xs(2))^2+(ys(1)-ys(2))^2))^(-3/2);
% Forces on magnet 2 (movable) due to magnet 1 (fixed)
Fxn = (xn(2)-xn(1))*rnn-(xn(2)-xs(1))*rsn;
Fyn = (yn(2)-yn(1))*rnn-(yn(2)-ys(1))*rsn;
Fxs =-(xs(2)-xn(1))*rns+(xs(2)-xs(1))*rss;
Fys =-(ys(2)-yn(1))*rns+(ys(2)-ys(1))*rss;
Ftoty = mu1*mu2*(Fyn + Fys) - m2*gee
%find mu1 so magnetic repulsion balances gravity
if Ftoty ~= 0
if (mu2*(Fyn + Fys)) ~=0
mu1 = (m2*gee)/(mu2*(Fyn + Fys))
Ftoty = 0;
else
fprintf('impossible to balance')
break
end
end
Ftotx = mu1*mu2*(Fxn + Fxs);
%plot force through midpoint and compute torque
Tx = cross([Fxn-Fxs,Fyn-Fys,0],.5*[xn(2)-xs(2),yn(2)-ys(2),0]);
fprintf('net force and torque on magnet 2 \n')
fprintf(' force %12.5f %12.5f torque %12.5f \n',Ftotx,Ftoty,Tx(3))
quiver((xn(2)+xs(2))/2,(yn(2)+ys(2))/2,3*Ftotx,3*Ftoty,'r')
% Plot equipotential curves:
[xx,yy] = meshgrid(xgrid,ygrid);
p1 = ((xx-xn(1)).^2+(yy-yn(1)).^2).^(-1/2);
p2 = ((xx-xn(2)).^2+(yy-yn(2)).^2).^(-1/2);
p3 = ((xx-xs(1)).^2+(yy-ys(1)).^2).^(-1/2);
p4 = ((xx-xs(2)).^2+(yy-ys(2)).^2).^(-1/2);
pp = mu1*(p1 - p3) + mu2*(p2 - p4);
levels = 10.^(-4:4);
levels = [-levels,levels];
contour(xx,yy,pp,levels), title(s), colormap(1e-6*[1 1 1]);
x1=xx(1:10:301,1:10:301);y1=yy(1:10:301,1:10:301);
p13=p1(1:10:301,1:10:301).^3; p23=p2(1:10:301,1:10:301).^3;
p33=p3(1:10:301,1:10:301).^3; p43=p4(1:10:301,1:10:301).^3;
Bx = mu1*(-(x1-xn(1)).*p13 + (x1-xs(1)).*p33)...
+ mu2*(-(x1-xn(2)).*p23 + (x1-xs(2)).*p43);
By = mu1*(-(y1-yn(1)).*p13 +(y1-ys(1)).*p33)...
+ mu2*(-(y1-yn(2)).*p23 +(y1-ys(2)).*p43);
quiver(x1,y1,Bx,By)
hold off