所望の力から電流を求めるコード書けた

Author: MASAHIRO KUBOTA

Ampere2.m

@@ -1,4 +1,4 @@
-function [F, T] = Ampere2(I2_x, I2_y, I2_z, I1_x, I1_y, I1_z, a, N, x, y, z, p1, q1, l1, p2, q2, l2, splitB, splitA)
+function [F, T] = Ampere2(I2_x, I2_y, I2_z, I1_x, I1_y, I1_z, a, N, x, y, z, quat1, quat2, splitB, splitA)
 %原点にあるN巻き，半径a，電流Iのコイルによって，(x, y, z)にあるN巻き，半径a，電流I_x, I_y, I_zのコイルに発生するアンペール力（力，トルク）
 %   Detailed explanation goes here
 
@@ -8,23 +8,22 @@
 phi = 0;
 F = [0, 0, 0];
 T = [0, 0, 0];
-quat = quaternion([p2,q2,l2],'euler','XYZ','point');
 
 
 %z coil
 while i < splitA
 
     i = i + 1;
     phi = phi + d_phi;
-    rot_point_z = rotatepoint(quat, [a * cos(phi), a * sin(phi), 0]);
-    B = magnetic_flux_three_coil(rot_point_z(1) + x, rot_point_z(2) + y, rot_point_z(3) + z, I1_x, I1_y, I1_z, a, N, p1, q1, l1, splitB);
+    rot_point_z = rotatepoint(quat2, [a * cos(phi), a * sin(phi), 0]);
+    B = magnetic_flux_three_coil(rot_point_z(1) + x, rot_point_z(2) + y, rot_point_z(3) + z, I1_x, I1_y, I1_z, a, N, quat1, splitB);
     %disp(norm(B))
-    d_I = I2_z * N * rotatepoint(quat, [-sin(phi), cos(phi), 0]) * d_phi;
+    d_I = I2_z * N * rotatepoint(quat2, a * [-sin(phi), cos(phi), 0]) * d_phi;
     %disp(I2_z)
     %disp(B)
     d_F = cross(d_I, B);
     %d_F = Ampere(I, phi, B);
-    d_T = cross(rotatepoint(quat, [a*cos(phi), a*sin(phi), 0]), d_F);
+    d_T = cross(rotatepoint(quat2, [a*cos(phi), a*sin(phi), 0]), d_F);
     %disp(d_B)
     F = F + d_F;
     T = T + d_T;
@@ -37,13 +36,13 @@
 
     i = i + 1;
     phi = phi + d_phi;
-    rot_point_y = rotatepoint(quat, [a * sin(phi), 0, a * cos(phi)]);
-    B = magnetic_flux_three_coil(rot_point_y(1) + x, rot_point_y(2) + y, rot_point_y(3) + z, I1_x, I1_y, I1_z, a, N, p1, q1, l1, splitB);
-    d_I = I2_y  * N* rotatepoint(quat, [cos(phi), 0, -sin(phi)] )* d_phi;
+    rot_point_y = rotatepoint(quat2, [a * sin(phi), 0, a * cos(phi)]);
+    B = magnetic_flux_three_coil(rot_point_y(1) + x, rot_point_y(2) + y, rot_point_y(3) + z, I1_x, I1_y, I1_z, a, N, quat1, splitB);
+    d_I = I2_y  * N* rotatepoint(quat2, a * [cos(phi), 0, -sin(phi)] )* d_phi;
 
     d_F = cross(d_I, B);
     %d_F = Ampere(I, phi, B);
-    d_T = cross(rotatepoint(quat, [a*sin(phi), 0, a*cos(phi)]), d_F);
+    d_T = cross(rotatepoint(quat2, [a*sin(phi), 0, a*cos(phi)]), d_F);
     %disp(d_B)
     F = F + d_F;
     T = T + d_T;
@@ -58,14 +57,14 @@
 
     i = i + 1;
     phi = phi + d_phi;
-    rot_point_x = rotatepoint(quat, [0, a * cos(phi), a * sin(phi)]);
-    B = magnetic_flux_three_coil(rot_point_x(1) + x, rot_point_x(2) + y, rot_point_x(3) + z, I1_x, I1_y, I1_z, a, N, p1, q1, l1, splitB); 
-    d_I = I2_x * N* rotatepoint(quat, [0, -sin(phi), cos(phi)]) * d_phi;
+    rot_point_x = rotatepoint(quat2, [0, a * cos(phi), a * sin(phi)]);
+    B = magnetic_flux_three_coil(rot_point_x(1) + x, rot_point_x(2) + y, rot_point_x(3) + z, I1_x, I1_y, I1_z, a, N, quat1, splitB); 
+    d_I = I2_x * N* rotatepoint(quat2, a * [0, -sin(phi), cos(phi)]) * d_phi;
     %d_I = I2_x * N* [0, -sin(phi), cos(phi)]  * d_phi;
     d_F = cross(d_I, B);
     %d_F = Ampere(I, phi, B);
     %disp(d_F)
-    d_T = cross(rotatepoint(quat, [0, a*cos(phi), a*sin(phi)]), d_F);
+    d_T = cross(rotatepoint(quat2, [0, a*cos(phi), a*sin(phi)]), d_F);
     %disp(d_B)
     F = F + d_F;
     %disp("d_F is")

README.md

@@ -3,18 +3,28 @@ READMEに，使えるファイルを書かないとなー
 
 READMEにどのファイルどうつかったらどういうことができるかまとめないとなー
 
+Ampere2.mは原点にあるN巻き，半径a，電流Iのコイルによって，(x, y, z)にあるN巻き，半径a，電流I_x, I_y, I_zのコイルに発生するアンペール力（力，トルク）を計算．
+coil_FT.mでいろんな相対姿勢でビオサバールとアンペールで電磁力とトルクを計算する．
+
+dynamics2record.mで任意の場所から10^-5Nの上限でレコード盤軌道にもって行く制御
+
+dynamic2record_movie.mで1つの衛星がレコード盤軌道に投入されるシミュレーション動画．dynamics2record.mを改造．フィードバックゲインに注意しろ！
+
+force2current.mで所望の力を電流に変換する．
 
-plot_magnetic_field_FT2で特定の位置姿勢，電流，分割電流数での相手のコイルにかかる力とトルクを計算
-coil_FT.mで，色んな姿勢パターンでコイルに加わる力とトルクを算出．
+plot_magnetic_field_FT2.mで特定の位置姿勢，電流，分割電流数での相手のコイルにかかる力とトルクを計算
+plot_magnetic_field_FT2_dipole.mで特定の位置姿勢，電流，分割電流数での相手のコイルにかかる力とトルクをダイポール近似で計算
 
 test3.mでトルクがほぼ0になる値をプロット
 test8.mで地球の磁場を計算
 test9.mはchatGPTで生成したコード．matlabで動画を生成し，保存するためのコード．
+test10.mでコイルjの電流要素から見た，コイルkの電流要素の位置ベクトルを表示（I_ij確かめるために使った）.plot_magnetic_field_FT2.mと同じ．計算の仕方が違う．逆問題を解くための式のバグ取りのため．
+test12.mはtest10.mを改造して，特定の力を出すために必要な電流を出力するためのコード
+test11.mは特定の力を出すために必要な電流を出力して，その値をplot_magnetic_field_FT2.mで検証している．最大電流も定義してるから，実際に出せる力がこのコードでわかる．
+
+
+
 
-coil_FT.mでいろんな相対姿勢でビオサバールとアンペールで電磁力とトルクを計算する．
-plot_magnetic_field_FT2.mで決められた相対姿勢，電流で発生る電磁力とトルクを計算する．
 
-dynamics2record.mで任意の場所から10^-5Nの上限でレコード盤軌道にもって行く制御
 
-dynamic2record_movie.mで1つの衛星がレコード盤軌道に投入されるシミュレーション動画．dynamics2record.mを改造．
 

calculateI.asv

@@ -0,0 +1,64 @@
+function I_jk = calculateI(x_j, x_k, q_j, q_k, split, a, j_coil, k_coil)
+%２つのコイルによって発生する電磁力の周回積分項を計算
+%   Detailed explanation goes here
+
+%{
+if j_coil == 1
+    E = [0,pi/2,0];
+    q_x = quaternion(E,'euler','XYZ','point');
+    q_j = q_x*q_j;
+end
+
+if j_coil == 2
+    E = [-pi/2,0,0];
+    q_x = quaternion(E,'euler','XYZ','point');
+    q_j = q_x*q_j;
+end
+
+if k_coil == 1
+    E = [0,pi/2,0];
+    q_x = quaternion(E,'euler','XYZ','point');
+    q_k = q_x*q_k;
+end
+
+
+if k_coil == 2
+    E = [-pi/2,0,0];
+    q_x = quaternion(E,'euler','XYZ','point');
+    q_k = q_x*q_k;
+end
+
+
+
+
+d_phi = 2*pi/split;
+d_theta = 2*pi/split;
+
+I_jk = zeros(1,3);
+phi = 0;
+
+for l = 1:split
+    phi = phi + d_phi;
+    theta = 0;
+    II_ij = zeros(1,3);
+    for m = 1:split
+        theta = theta + d_theta;
+        r_jk = calculater_jk(x_j, x_k, q_j, q_k, theta, phi, a);
+
+        %磁場の向きdispすると向きが合ってるかのデバッグができる．
+        II_ij = II_ij + 1/norm(r_jk)^3*cross([-a*sin(theta), a*cos(theta), 0], r_jk)*d_theta;
+        
+        %disp(1/norm(r_jk)^3*cross(r_jk, [-a*sin(phi), a*cos(phi), 0]))
+    end
+
+    q_jk = quatmultiply(q_k, quatconj(q_j));
+    %disp(rad2deg(phi))
+    %disp(r_jk)
+    %disp(II_ij)
+    %disp(cross(II_ij, rotatepoint(q_jk, [-a*sin(theta), a*cos(theta), 0])*d_theta))
+    rt = rotatepoint(q_jk, [-a*sin(phi), a*cos(phi), 0]);
+    I_jk = I_jk + cross(rt, II_ij)*d_phi;
+    %disp(x_k - x_j + rotatepoint(q_jk, [-a*sin(theta), a*cos(theta), 0]))
+    %disp(I_jk)
+    
+end

calculateI.m

@@ -0,0 +1,67 @@
+function I_jk = calculateI(x_j, x_k, q_j, q_k, split, a, j_coil, k_coil)
+%２つのコイルによって発生する電磁力の周回積分項を計算
+%   Detailed explanation goes here
+
+
+if j_coil == 1
+    E = [0,pi/2,0];
+    q_x = quaternion(E,'euler','XYZ','point');
+    q_j = q_j*q_x;
+end
+
+if j_coil == 2
+    E = [-pi/2,0,0];
+    q_x = quaternion(E,'euler','XYZ','point');
+    q_j = q_j*q_x;
+end
+
+if k_coil == 1
+    E = [0,pi/2,0];
+    q_x = quaternion(E,'euler','XYZ','point');
+    q_k = q_k*q_x;
+end
+
+
+if k_coil == 2
+    E = [-pi/2,0,0];
+    q_x = quaternion(E,'euler','XYZ','point');
+    q_k = q_k*q_x;
+end
+
+
+
+
+d_phi = 2*pi/split;
+d_theta = 2*pi/split;
+
+I_jk = zeros(1,3);
+phi = 0;
+
+for l = 1:split
+    phi = phi + d_phi;
+    theta = 0;
+    II_ij = zeros(1,3);
+    for m = 1:split
+        theta = theta + d_theta;
+        r_jk = calculater_jk(x_j, x_k, q_j, q_k, theta, phi, a);
+
+        %磁場II_ijの向きdispすると向きが合ってるかのデバッグができる．
+        II_ij = II_ij + 1/norm(r_jk)^3*cross(rotatepoint(q_j, [-a*sin(theta), a*cos(theta), 0]), r_jk)*d_theta;
+        %disp(r_jk)
+        %disp(1/norm(r_jk)^3*cross(r_jk, [-a*sin(phi), a*cos(phi), 0]))
+    end
+
+    q_jk = quatmultiply(q_k, quatconj(q_j));
+    
+    %disp(rad2deg(phi))
+    %disp(r_jk)
+    %disp(II_ij)
+    %disp(cross(II_ij, rotatepoint(q_jk, [-a*sin(theta), a*cos(theta), 0])*d_theta))
+    %disp(q_k)
+    rt = rotatepoint(q_k, [-a*sin(phi), a*cos(phi), 0]);
+    I_jk = I_jk + cross(rt, II_ij)*d_phi;
+    
+    %disp(x_k - x_j + rotatepoint(q_jk, [-a*sin(theta), a*cos(theta), 0]))
+    %disp(I_jk)
+    
+end

calculater_jk.m

@@ -0,0 +1,22 @@
+function r_jk = calculater_jk(x_j, x_k, q_j, q_k, theta, phi, a)
+%コイルjの電流要素から見た，コイルkの電流要素の位置ベクトル
+%   Detailed explanation goes here
+%a = 0.015;
+%x_j = [0.1, 0.1, 0.1];
+%E_j = [pi/3,0,pi/4];
+%q_j = quaternion(E_j,'euler','XYZ','point');
+%theta = deg2rad(120);
+%dl_j = rotatepoint(q_j, [-a*sin(theta), a*cos(theta), 0]);
+x_jR = x_j + rotatepoint(q_j, [a*cos(theta), a*sin(theta), 0]);
+%x_jR0 = x_j + rotatepoint(q_j, [a, 0, 0]);
+
+%x_k = [-0.1, -0.1, -0.1];
+%E_k = [-pi/3,0,-pi/4];
+%q_k = quaternion(E_k,'euler','XYZ','point');
+%phi = deg2rad(-120);
+%dl_k = rotatepoint(q_k, [-a*sin(phi), a*cos(phi), 0]);
+x_kR = x_k + rotatepoint(q_k, [a*cos(phi), a*sin(phi), 0]);
+%x_kR0 = x_k + rotatepoint(q_k, [a, 0, 0]);
+
+r_jk = x_kR - x_jR;
+end

draw_coil.m

@@ -13,4 +13,5 @@ function draw_coil(x, q, a)
 hold on
 plot3(coil2_y(1,:) + x0, coil2_y(2,:) + y0, coil2_y(3,:) + z0)
 plot3(coil2_z(1,:) + x0, coil2_z(2,:) + y0, coil2_z(3,:) + z0)
+hold on
 end

dynamic2record_movie.m

@@ -174,6 +174,7 @@
 
 %各時刻での電磁力の大きさをグラフ化
 figure
+str='';
 plot((1:N)/60, u_data)
 xlabel('時間（min）');
 ylabel('力（N）');