forked from andrewgiuliani/matplasmaopt
-
Notifications
You must be signed in to change notification settings - Fork 0
/
biotsavartAndGrad.m
121 lines (95 loc) · 5.72 KB
/
biotsavartAndGrad.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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
%% Biot-Savart integration on a generic curve
function [B1, B2, B3,...
dB1_dx, dB2_dx, dB3_dx,...
dB1_dy, dB2_dy, dB3_dy,...
dB1_dz, dB2_dz, dB3_dz]=biotsavartAndGrad(x,y,z, coilData)
coil_points = coilData.coil_points;
I = coilData.I;
coil_tangents = coilData.coil_tangents;
nCoils = coilData.C;
nfp = coilData.nfp;
ss = coilData.ss;
%% Declaration of variables
gamma = 4 * pi * 10^(-7) ;
dB = zeros(size(coil_points,1),3);
dB_dx = zeros(size(coil_points,1),3);
dB_dy = zeros(size(coil_points,1),3);
dB_dz = zeros(size(coil_points,1),3);
point = [x y z];
for i = 1: nCoils
%% Numerical integration of Biot-Savart law
diff = point-coilData.coil_field{i};
dist = sqrt(diff(:,1).^2 + diff(:,2).^2 + diff(:,3).^2);
dist_x = diff(:,1)./dist;
dist_y = diff(:,2)./dist;
dist_z = diff(:,3)./dist;
d_point_dx = repmat([1,0,0], size(coilData.coil_tangents{i},1), 1);
d_point_dy = repmat([0,1,0], size(coilData.coil_tangents{i},1), 1);
d_point_dz = repmat([0,0,1], size(coilData.coil_tangents{i},1), 1);
dB = dB + (gamma * I(i) ) / (4. * pi) * cross(coil_tangents{i}, diff)./dist.^3;
dB_dx = dB_dx + (gamma * I(i) ) / (4. * pi) * ( cross(coil_tangents{i}, d_point_dx).*dist.^3 ...
-3*dist.^2 .* dist_x .* cross(coil_tangents{i}, diff) ) ./ dist.^6;
dB_dy = dB_dy + (gamma * I(i) ) / (4. * pi) * ( cross(coil_tangents{i}, d_point_dy).*dist.^3 ...
-3*dist.^2 .* dist_y .* cross(coil_tangents{i}, diff) ) ./ dist.^6;
dB_dz = dB_dz + (gamma * I(i) ) / (4. * pi) * ( cross(coil_tangents{i}, d_point_dz).*dist.^3 ...
-3*dist.^2 .* dist_z .* cross(coil_tangents{i}, diff) ) ./ dist.^6;
if ss == 1
diff = point-coilData.coil_field{coilData.nfp*coilData.C + i};
dist = sqrt(diff(:,1).^2 + diff(:,2).^2 + diff(:,3).^2);
dist_x = diff(:,1)./dist;
dist_y = diff(:,2)./dist;
dist_z = diff(:,3)./dist;
dB = dB + (gamma * -I(i) ) / (4. * pi) * cross(coil_tangents{nfp*nCoils + i}, diff)./dist.^3;
dB_dx = dB_dx + (gamma * -I(i) ) / (4. * pi) * ( cross(coil_tangents{nfp*nCoils + i}, d_point_dx).*dist.^3 ...
-3*dist.^2 .* dist_x .* cross(coil_tangents{nfp*nCoils + i}, diff) ) ./ dist.^6;
dB_dy = dB_dy + (gamma * -I(i) ) / (4. * pi) * ( cross(coil_tangents{nfp*nCoils + i}, d_point_dy).*dist.^3 ...
-3*dist.^2 .* dist_y .* cross(coil_tangents{nfp*nCoils + i}, diff) ) ./ dist.^6;
dB_dz = dB_dz + (gamma * -I(i) ) / (4. * pi) * ( cross(coil_tangents{nfp*nCoils + i}, d_point_dz).*dist.^3 ...
-3*dist.^2 .* dist_z .* cross(coil_tangents{nfp*nCoils + i}, diff) ) ./ dist.^6;
end
for t = 1:nfp-1
diff = point-coilData.coil_field{t*coilData.C + i};
dist = sqrt(diff(:,1).^2 + diff(:,2).^2 + diff(:,3).^2);
dist_x = diff(:,1)./dist;
dist_y = diff(:,2)./dist;
dist_z = diff(:,3)./dist;
dB = dB + (gamma * I(i)) / (4. * pi) * cross(coil_tangents{t*nCoils + i}, diff)./dist.^3;
dB_dx = dB_dx + (gamma * I(i) ) / (4. * pi) * ( cross(coil_tangents{t*nCoils + i}, d_point_dx).*dist.^3 ...
-3*dist.^2 .* dist_x .* cross(coil_tangents{t*nCoils + i}, diff) ) ./ dist.^6;
dB_dy = dB_dy + (gamma * I(i) ) / (4. * pi) * ( cross(coil_tangents{t*nCoils + i}, d_point_dy).*dist.^3 ...
-3*dist.^2 .* dist_y .* cross(coil_tangents{t*nCoils + i}, diff) ) ./ dist.^6;
dB_dz = dB_dz + (gamma * I(i) ) / (4. * pi) * ( cross(coil_tangents{t*nCoils + i}, d_point_dz).*dist.^3 ...
-3*dist.^2 .* dist_z .* cross(coil_tangents{t*nCoils + i}, diff) ) ./ dist.^6;
if ss == 1
diff = point-coilData.coil_field{coilData.nfp*coilData.C +t*coilData.C + i};
dist = sqrt(diff(:,1).^2 + diff(:,2).^2 + diff(:,3).^2);
dist_x = diff(:,1)./dist;
dist_y = diff(:,2)./dist;
dist_z = diff(:,3)./dist;
dB = dB + (gamma * -I(i)) / (4. * pi) * cross(coil_tangents{nfp*nCoils +t*nCoils + i}, diff)./dist.^3;
dB_dx = dB_dx + (gamma * -I(i) ) / (4. * pi) * ( cross(coil_tangents{nfp*nCoils +t*nCoils + i}, d_point_dx).*dist.^3 ...
-3*dist.^2 .* dist_x .* cross(coil_tangents{nfp*nCoils +t*nCoils + i}, diff) ) ./ dist.^6;
dB_dy = dB_dy + (gamma * -I(i) ) / (4. * pi) * ( cross(coil_tangents{nfp*nCoils +t*nCoils + i}, d_point_dy).*dist.^3 ...
-3*dist.^2 .* dist_y .* cross(coil_tangents{nfp*nCoils +t*nCoils + i}, diff) ) ./ dist.^6;
dB_dz = dB_dz + (gamma * -I(i) ) / (4. * pi) * ( cross(coil_tangents{nfp*nCoils +t*nCoils + i}, d_point_dz).*dist.^3 ...
-3*dist.^2 .* dist_z .* cross(coil_tangents{nfp*nCoils +t*nCoils + i}, diff) ) ./ dist.^6;
end
end
end
Bfield = sum(dB,1) * 2 * pi / size(dB,1);
B1 = Bfield(1);
B2 = Bfield(2);
B3 = Bfield(3);
B_dx = sum(dB_dx,1) * 2 * pi / size(dB_dx,1);
B_dy = sum(dB_dy,1) * 2 * pi / size(dB_dy,1);
B_dz = sum(dB_dz,1) * 2 * pi / size(dB_dz,1);
dB1_dx = B_dx(1);
dB1_dy = B_dy(1);
dB1_dz = B_dz(1);
dB2_dx = B_dx(2);
dB2_dy = B_dy(2);
dB2_dz = B_dz(2);
dB3_dx = B_dx(3);
dB3_dy = B_dy(3);
dB3_dz = B_dz(3);
end