Quasistatic Cavity Resonance for Ubiquitous Wireless Power Transfer (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;
}