C言語初心者。静磁場の計算。台形公式を用いた重積分がうまくいかない

永久磁石がある点に作る磁束密度をC言語にて計算したいと考えています。

計算がうまくいかず、どこに間違いがあるのか見つけられていないので、
電磁気ないしC言語に見識のある方のお力添えを頂きたく存じます。

表面磁場 M = 0.24 T
磁石長さ L = 1 mm (y軸方向)
磁石幅 W = 1 mm (z軸方向)
の磁石がある点に作る磁場を求めます。

磁石表面を長さ方向(y軸方向)にn分割してdy、幅方向(z軸方向)にm分割してdzとし、
それぞれの領域が作る磁場をyとzに関する重積分で重畳して求めようとしています。

この時、磁石の表面磁束密度を M で一様だと仮定すると、分割された各要素の
持つ磁極の強さは dy*dz*M [Wb] となります。

したがって、磁場のクーロンの法則より、各要素から r [m]離れた点に磁石の要素が作る磁場dHは、

dH = (1/4πµo) * {(dy*dz*M) / r^2}

磁束密度dB [T]に直して、

dB = (1/4π) * {(dy*dz*M) / r^2)}

ここで、
要素の座標が ( xs0 , ys0 , zs0 ) (添え字はソース点の意味でsです。読みにくくてすみません)
磁束密度を求める点が ( x , y , z )
だとすると、2点間の距離 r = √{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2} なので

dB = (1/4π) * [(dy*dz*M) /{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}]

さらに、x軸方向の磁束密度を求めたいため

dB = (1/4π) * [(dy*dz*M) /{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}]
* (x-xs0) / √{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}

としてそれぞれ求め、これを
y軸方向に 0 → 10^(-3) [m]
z軸方向に 0 → 10^(-3) [m]
の範囲で積分して計算しています。

ここまでで間違いはございますでしょうか。
簡単な静磁場の計算ですが、得手でなくお恥ずかしい限りです。

積分には台形公式を用いています。以下、実際のソースコードになります。

このプログラムで計算すると、解が異常に小さく出たり、分割数を変えると解の大きさが
1桁単位で変化したりと、明らかに間違っており、おそらくどこか基礎的なことで躓いている
ことは予想できるのですが、具体的にどこが間違っているのかが分かりません。

見識のある方のお力添えを頂きたく存じます。

#define _USE_MATH_DEFINES
#include<stdio.h>
#include<math.h>
#define f(ys0,zs0) ((x-xs0)*M)/(4*M_PI)/((x-xs0)*(x-xs0)+(y-ys0)*(y-ys0)+(z-zs0)*(z-zs0))/sqrt((x-xs0)*(x-xs0)+(y-ys0)*(y-ys0)+(z-zs0)*(z-zs0))

int main(void) {
int n, m, i;
double L1, L2, W1, W2, M, xs0, xs1, x, ys0, ys1, y;
double zs0, zs1, z, dy, dz, dSy1, dSy2, Sy1, Sy2, dV, V;

//要素分割数
n = 1000;
m = 1000;
//-- パラメータの設定
L1 = 0.000;//yの下限値(磁石の長さ)
L2 = 0.001;//yの上限値
W1 = 0.000;//zの下限値(磁石の幅)
W2 = 0.001;//zの上限値
M = 0.240;//磁石の磁界強度 [T] = [Wb/m^2]
xs0 = 0.000;//ソース点
x = 0.001;//目標点
y = 0.000;//目標点
z = 0.000;//目標点
//--
//微分量の定義
dy = (L2 - L1) / n;
dz = (W2 - W1) / m;
//初期化
V = 0;

//-- 積分
for (i = 1; i < (n-1); i++) {
//y軸ソース位置の指定
ys0 = L1 + dy * i;
ys1 = ys0 + dy;
//ループ毎の初期化
Sy1 = 0;
Sy2 = 0;
for (i = 1; i < (m-1); i++) {
//z軸ソース位置の指定
zs0 = W1 + dz * i;
zs1 = zs0 + dz;
//y軸方向の微小面積計算
dSy1 = (f(ys0, zs0) + f(ys0, zs1)) * dz / 2;
dSy2 = (f(ys1, zs0) + f(ys1, zs1)) * dz / 2;
printf("dSy1 = %lf\ndSy2 = %lf\n", dSy1, dSy2);
//y軸方向の面積総和
Sy1 += dSy1;
Sy2 += dSy2;
}
//z軸方向の微小体積積分
dV = (Sy1 + Sy2)*dy / 2;
//z軸方向の体積総和
V += dV;
}
//-- 積分終了

printf("磁束密度 %lf [T]\n", V);
getchar();
}

No.1 ベストアンサー 回答者： Tacosan 回答日時： 2019/01/12 00:09

ネストしている for で同じ変数を使っているのでぐっちゃぐちゃになっている?

この回答へのお礼

解決致しました。ご指摘頂き誠にありがとうございます。非常に初歩的なところだったので、反省しきりです。

お礼日時：2019/01/12 17:38