[Quasistatic Cavity Resonance for Ubiquitous Wireless Power Transfer
(PDFファイル)](https://s3-us-west-1.amazonaws.com/disneyresearch/wp-
content/uploads/20170215220933/Quasistatic-Cavity-Resonance-for-Ubiquitous-
Wireless-Power-Transfer-Paper.pdf)

こちらの論文の効率を導出するプログラムをC言語で再現したいのですが出力結果が望んだものから大きく外れた値をとってしまいます．  
個人的には，alpha,betaあたりのベクトルの入った積分が怪しいと考えています．

```

    #define _CRT_SECURE_NO_WARNINGS
    #define _USE_MATH_DEFINES
    
    #include <stdio.h>
    #include <string.h>
    #include <stdlib.h>
    #include <math.h>
    
    //関数
    double hfield_x(double I1, double w1 ,double x1,double y1 ) {   //論文(6)式(キャビティ内のあるx,y座標での磁界の大きさを求める式のｘ方向のベクトルaxの係数）
        double h_x=0;
        int m=0, l=0;
    
        for (m = -1;m <= 1;m++) {
            for (l = -1;l <= 1;l++) {
                h_x+= -I1 * pow(-1, m + l)*(y1 - m * w1) / (pow(x1 - l * w1, 2) + pow(y1 - m * w1, 2));
            }
        }
        h_x=h_x/ (2 * M_PI);
    
        return h_x;
    }
    
    double hfield_y(double I1, double w1, double x1, double y1) {   //論文(6)式(キャビティ内のあるx,y座標での磁界の大きさを求める式のｙ方向のベクトルayの係数）
        double h_y=0;
        int m=0, l=0;
    
        for (m = -1;m <= 1;m++) {
            for (l = -1;l <= 1;l++) {
                h_y += I1 * pow(-1, m + l)*(x1 - l * w1) / (pow(x1 - l * w1, 2) + pow(y1 - m * w1, 2));
            }
        }
        h_y = h_y / (2 * M_PI);
    
        return h_y;
    }
    
    double alpha(double hfx, double hfy,double myu1) {  //論文(9)式(キャビティ内に蓄積された総磁気エネルギーを求める式)　
        double ans1 = (myu1 / 2)*(hfx * hfx + hfy * hfy);
    
        return ans1;
    }
    
    double beta(double hfx, double hfy, double myu1,double nx1,double ny1) {    //論文(10)式(受電コイルの捕捉する磁束)
        double ans2 = myu1 * (hfx*nx1 + hfy * ny1);
    
        return ans2;
    }
    
    
    int main(void) {
        int x = 0, z = 0;   /*直交座標*/
        double y = 0;
        int m = -1, l = -1; /*Σ計算の刻み値*/
        double I = 0.0, sum = 0.0;  //電流読み取り用,積分用
        double Q1 = 2130.0, Q2 = 360.0; /*送受信のQ値*/
        double v = 54.0, A = 2.720, w = 4.90, h = 2.30, Aw = 0.1650;    /*基本寸法*/
        double Hx = 0.0, Hy = 0.0, hx = 0.0, hy = 0.0;  /*磁界　各ベクトル,繰り返し用*/
        double a = 0.0, b = 0.0;    /*α,β*/
        double K = 0.0, X = 0.0;    /*κ, χ*/
        double Gmax[10] = { 0.0 };  /*効率*/
        double f = 1.32*pow(10, 6); //共振周波数
        double myu = 4 * M_PI*pow(10, -7);  //透磁率
        int n = 100;    //積分の細かさ
        double i = 0, j = 0, k = 0, dx = 0, dy = 0, dz = 0, d1 = 0, d2 = 0; //積分の刻み値
        double nx = 0.0, ny = 0.0;  //単位法線ベクトル
        double L2 = 0.0019910;  //受信コイルのインダクタンス(C2=7.3pF)
        double sy1 = 0, sy2 = 0, s1 = 0, s2 = 0, sx = 0;//台形公式　保存用
        double xlen[8] = { 1.5240000e-01 ,4.5720000e-01 ,7.6200000e-01 ,1.0668000 ,1.3716000,1.6764000,1.9812000,2.2860000 };   //効率を計算するx座標
    
        L2 = 1.31e-6;
    
        dx = w / n;
        dy = w / n;
        dz = h / n;
        d1 = Aw / n;
        d2 = Aw / n;
    
    
        printf("I0を入力\n");
        scanf("%lf", &I);
    
        for (i = -w / 2;i < w / 2;i += dx) {    //alphaの計算
            sy1 = 0;
            sy2 = 0;
    
            for (j = -w / 2;j < w / 2;j += dy) {
                Hx = 0, Hy = 0;
    
                if (i == 0 && j == 0) {
                    break;
                }
                else {
    
                    s1 = (alpha(hfield_x(I, w, i, j), hfield_y(I, w, i, j), myu) + alpha(hfield_x(I, w, i, j + dy), hfield_y(I, w, i, j + dy), myu)) / 2 * dy;
                    s1 = (alpha(hfield_x(I, w, i + dx, j), hfield_y(I, w, i + dx, j), myu) + alpha(hfield_x(I, w, i + dx, j + dy), hfield_y(I, w, i + dx, j + dy), myu)) / 2 * dy;
    
                    sy1 += s1;
                    sy2 += s2;
                }
    
            }
    
            sx = (sy1 + sy2) / 2 * dx;
            a += sx;
    
        }
        a = a * h;
    
    
        for (x = 0;x <= 7;x++) {
    
            Hx = 0, Hy = 0, hx = 0, hy = 0;
            Hx = hfield_x(I, w, xlen[x], 0);
            Hy = hfield_y(I, w, xlen[x], 0);
    
            nx = 0;
            ny = 1;
    
            b = 0;
    
            for (i = xlen[x] - Aw / 2;i < xlen[x] + Aw / 2;i += d1) {   //betaの計算 
    
                Hx = 0, Hy = 0;
                sy1 = 0;
                sy2 = 0;
    
                s1 = (beta(hfield_x(I, w, i, 0), hfield_y(I, w, i, 0), myu, nx, ny) + beta(hfield_x(I, w, i+d1, 0), hfield_y(I, w, i+d1, 0), myu, nx, ny)) / 2 * d1;
    
                b += s1;
            }
            b = b * Aw;
    
            K = sqrt(2) * 2 * M_PI*f*b / (4 * sqrt(L2*a));  //κの計算(8)式
    
            X = 4 * Q1*Q2*pow(fabs(K), 2) / pow(2 * M_PI*f, 2); //χの計算(11)式
    
            Gmax[x] = X / pow(1 + sqrt(1 + X), 2);  //Gmaxの計算　目的の効率を導出する式
    
            printf("%.50f \n", Gmax[x]);
    
        }
    
        return 0;
    }
    
```
CC BY-SA 4.0
2018-12-21T05:43:06.547
51409
2019-03-10T13:50:29.463
2018-12-22T12:54:29.333
76
31542
question
c
目標の出力が出てこない
学術面にこたえることはできないので、ヒントになればと思い  
回答します。

このぐらいの込み入ったプログラムを書くのであれば、cUnitなどを利用して単体テストを  
積み上げることをお勧めします。  
プログラムをできるだけ小さい関数の集まりにして関数ごとの動作が正しいことを積み上げて行くやり方です。また、全体の動作の中で「どこで値が間違っているのか」を知るためには、gdbなどのデバッカの利用が効果的です。

学術論文を読み込み、提示されたプログラムをよみこみ、どこにバグがあるかを答えるという作業は大変に「重い」作業です。回答してくださる方は特に質問者様のようにその分野に明るいわけではないので。。。回答を得やすくするために、質問の範囲を狭めてはいかがでしょうか。。
CC BY-SA 4.0
2018-12-22T21:46:50.553
51450
2018-12-22T21:46:50.553
10174
51409
answer
えっと、回答ではないですが・・・  
このソースコンパイル通らないですよ？？

```

    s51409.c:21:22: error: use of undeclared identifier 'M_PI'
        h_x = h_x / (2 * M_PI);
    
```

[予約済み識別子](https://www.wdic.org/w/TECH/%E4%BA%88%E7%B4%84%E6%B8%88%E3%81%BF%E8%AD%98%E5%88%A5%E5%AD%90)  
問題となる名前（アンダースコア＋識別子は予約語）が見受けられます。気をつけましょう。

```

    s51409.c:1:9: warning: macro name is a reserved identifier [-Wreserved-id-macro]
    #define _CRT_SECURE_NO_WARNINGS
            ^
    s51409.c:2:9: warning: macro name is a reserved identifier [-Wreserved-id-macro]
    #define _USE_MATH_DEFINES
    
```

こちらの環境はLinux mint clang version 9.0.0 (trunk 354955)  
また、私は大きなプログラムの場合は２分法のようなデバッグをします。  
処理を半分に分け何処で挙動がおかしくなるか調べる→おかしくなった部分をまた半分にして・・・的なｄ＾＾
CC BY-SA 4.0
2019-03-10T13:40:02.120
53344
2019-03-10T13:50:29.463
2019-03-10T13:50:29.463
25438
25438
51409
answer