-
Notifications
You must be signed in to change notification settings - Fork 1
/
cuboid_force_x_y.m
91 lines (71 loc) · 2.57 KB
/
cuboid_force_x_y.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
%% cuboid_force_x_y
%
% Calculate the forces between two parallel cuboid magnets, magnetised in the
% x- and y-directions respectively.
% \START
% \begin{mfunction}{cuboid_force_x_y}
function forces_xyz = cuboid_force_x_y(size1,size2,offset,J1,J2)
J1 = J1(1);
J2 = J2(2);
if ( abs(J1)<eps || abs(J2)<eps )
forces_xyz = [0; 0; 0];
return;
end
component_x = 0;
component_y = 0;
component_z = 0;
for ii = [1 -1]
for jj = [1 -1]
for kk = [1 -1]
for ll = [1 -1]
for pp = [1 -1]
for qq = [1 -1]
ind_sum = ii*jj*kk*ll*pp*qq;
u = -offset(3,:) + size2(3)*jj - size1(3)*ii;
v = offset(2,:) + size2(2)*ll - size1(2)*kk;
w = offset(1,:) + size2(1)*qq - size1(1)*pp;
r = sqrt(u.^2+v.^2+w.^2);
tmp = zeros(size(u));
atan_term_u = tmp; atan_term_v = tmp; atan_term_w = tmp;
log_ru = tmp; log_rw = tmp; log_rv = tmp;
ind = abs(u) > eps;
if any(ind), atan_term_u(ind) = atan(v(ind).*w(ind)./(r(ind).*u(ind))); end
ind = abs(v) > eps;
if any(ind), atan_term_v(ind) = atan(u(ind).*w(ind)./(r(ind).*v(ind))); end
ind = abs(w) > eps;
if any(ind), atan_term_w(ind) = atan(u(ind).*v(ind)./(r(ind).*w(ind))); end
ind = abs(r-u) > eps;
if any(ind), log_ru(ind) = log(r(ind)-u(ind)); end
ind = abs(r+w) > eps;
if any(ind), log_rw(ind) = log(r(ind)+w(ind)); end
ind = abs(r+v) > eps;
if any(ind), log_rv(ind) = log(r(ind)+v(ind)); end
cx = ...
+ v.*w.*log_ru ...
- v.*u.*log_rw ...
- u.*w.*log_rv ...
+ 0.5*u.^2.*atan_term_u ...
+ 0.5*v.^2.*atan_term_v ...
+ 0.5*w.^2.*atan_term_w;
cy = ...
- 0.5*(u.^2-v.^2).*log_rw ...
+ u.*w.*log_ru ...
+ u.*v.*atan_term_v ...
+ 0.5*w.*r;
cz = ...
- 0.5*(u.^2-w.^2).*log_rv ...
+ u.*v.*log_ru ...
+ u.*w.*atan_term_w ...
+ 0.5*v.* r;
component_x = component_x + ind_sum.*cx;
component_y = component_y + ind_sum.*cy;
component_z = component_z + ind_sum.*cz;
end
end
end
end
end
end
forces_xyz = J1*J2/(4*pi*(4*pi*1e-7))*[ component_z; component_y; -component_x ];
end
% \end{mfunction}