-
Notifications
You must be signed in to change notification settings - Fork 1
/
demo_dipole.m
76 lines (55 loc) · 1.22 KB
/
demo_dipole.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
%% Demo of dipole equations
%% Verify Yung 1998 example (Table 1)
ma = [100;0;0];
mb = [100;0;0];
rab = [1;0;0];
f1 = dipole_forcetorque(ma,mb,rab)
ma = [0;-100;0];
mb = [100;0;0];
rab = [1;0;0];
dipole_forcetorque(ma,mb,rab)
f2 = dipole_forcetorque(ma,mb,rab)
ma = [100/sqrt(2);0;100/sqrt(2)];
mb = [0;0;100];
rab = [1;0;0];
dipole_forcetorque(ma,mb,rab)
f3 = dipole_forcetorque(ma,mb,rab)
f3*sqrt(2)
ma = [0;0;100];
mb = [0;0;100];
rab = [1;0;0];
dipole_forcetorque(ma,mb,rab)
f4 = dipole_forcetorque(ma,mb,rab)
%% Offset dipoles in repulsion
ma = [0;0;+1];
mb = [0;0;+1];
z = linspace(0,0.5);
rab = [0;0;1]*z + [0.1;0;0];
[f5,t5] = dipole_forcetorque(ma,mb,rab);
figure(1); hold on; clf
yyaxis left
plot(z,f5(3,:))
xlabel('Displacement')
ylabel('Force')
yyaxis right
plot(z,t5(2,:))
ylabel('Torque')
%% Rotating dipoles
ma = [0;0;1];
R = 0.5;
t = linspace(0,pi);
mb = [-cos(t);zeros(size(t));sin(t)];
rab = [R*cos(t);zeros(size(t));R*sin(t)];
[f6,t6] = dipole_forcetorque(ma,mb,rab);
figure(1); hold on; clf
yyaxis left
plot(t,f6(3,:))
ylabel('Force')
yyaxis right
plot(t,t6(2,:))
ylabel('Torque')
grid on
xlabel('Angle, deg.')
xticks([0 pi/4 pi/2 3*pi/4 pi])
xticklabels({'0','\pi/4','\pi/2','3\pi/4','\pi'})
xlim(t([1 end]))