LM法を用いて非線形フィティングを行いたい

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define P 501
int N = (P + 1) / 2;

// 関数の定義
double f(double q, double J, double ganma) {
    double A = 2 * J * (1 - cos(q));
    double C = ganma / 2;
    double E = pow((A + C), 2) - pow(C, 2);
    double F = sqrt(E);

    return F;
}

// ヤコビアン行列を計算する関数
void calculateJacobian(double x[], double J[], double ganma[], int n, double J_val, double ganma_val) {
    for (int i = 0; i < n; i++) {
        double q = x[i];

        double A = (1-cos(q));
        double B = 4*J_val*A*A;
        double C = ganma_val*A;
        double G = +J_val*A;

        double E =4*J_val*J_val*A*A+2*J_val*ganma_val*A;
        double F = pow(E,0.5);

        J[i] = (A+B)/F;
        ganma[i] = G/F;
    }
}

// Levenberg-Marquardt法
void levenbergMarquardt(double x[], double y[], int n, double* J, double* ganma) {
    double lambda = 0.0000001; // 初期のパラメータ更新係数
    double epsilon = 1e-6;   // 収束判定の閾値
    int maxIterations = 100000; // 最大反復回数

    double J_val = *J;
    double ganma_val = *ganma;

    double residuals[P];    // 残差ベクトル
    double Jacobian[P];     // ヤコビアン行列
    double ganma_Jacobian[P];

    int iterations = 0;
    double prevCost = 0.0;

    while (iterations < maxIterations) {
        // 残差ベクトルの計算
        for (int i = 0; i < n; i++) {
            double q = x[i];
            residuals[i] = y[i] - f(q, J_val, ganma_val);
        }

        // ヤコビアン行列の計算
        calculateJacobian(x, Jacobian, ganma_Jacobian, n, J_val, ganma_val);

        // システム行列の計算
        double H[2][2] = {{0}};
        double b[2] = {0};

        for (int i = 0; i < n; i++) {
            H[0][0] += Jacobian[i] * Jacobian[i];
            H[0][1] += Jacobian[i] * ganma_Jacobian[i];
            H[1][0] += ganma_Jacobian[i] * Jacobian[i];
            H[1][1] += ganma_Jacobian[i] * ganma_Jacobian[i];

            b[0] += Jacobian[i] * residuals[i];
            b[1] += ganma_Jacobian[i] * residuals[i];
        }

        // パラメータの更新
        double H_sum = H[0][0] + H[1][1];
        double det = H[0][0] * H[1][1] - H[0][1] * H[1][0];

        double delta[2];
        delta[0] = (H_sum * b[0] - (H[0][1] * b[1])) / det;
        delta[1] = ((-H[1][0] * b[0]) + (H[0][0] * b[1])) / det;

        J_val += delta[0];
        ganma_val += delta[1];

        // 収束判定
        double currentCost = 0.0;
        for (int i = 0; i < n; i++) {
            currentCost += pow(residuals[i], 2);
        }

        if (fabs(currentCost - prevCost) < epsilon) {
            break;
        }

        if (currentCost < prevCost) {
            lambda /= 10.0;
        } else {
            lambda *= 10.0;
            J_val -= delta[0];
            ganma_val -= delta[1];
        }

        prevCost = currentCost;
        iterations++;
    }

    *J = J_val;
    *ganma = ganma_val;
}

int main() {
    FILE* fp;
    FILE* fp_1;
    FILE* fp_2;

    double x[(P + 1) / 2], y[(P + 1) / 2];

    // データポイント
    if ((fp = fopen("J.txt", "w")) == NULL) {
        printf("Cannot open the file\n");
        exit(1);
    }
    fp_1 = fopen("q.txt", "r");
    if (fp_1 == NULL) {
        printf("ファイルを開くことができませんでした．\n");
    }
    for (int l = 0; l < N; l++) {
        fscanf(fp_1, "%lf", &(x[l]));
    }
    fp_2 = fopen("omega.txt", "r");
    if (fp_2 == NULL) {
        printf("ファイルを開くことができませんでした．\n");
    }
    for (int i = 0; i < N; i++) {
        fscanf(fp_2, "%lf", &(y[i]));
    }

    // パラメータの初期値
    double J = 0.30;
    double ganma = 0.090;

    // Levenberg-Marquardt法でパラメータを最適化
    levenbergMarquardt(x, y, N, &J, &ganma);

    // 結果表示
    printf("J: %.6f\n", J);
    printf("ganma: %.6f\n", ganma);

    for (int i = 0; i < (P + 1); i++) {
        double p = M_PI / (P + 1) * i;
        double fc = f(p, J, ganma);
        double fc_1 = f(p, J, 0);
        fprintf(fp, "%lf %lf %lf\n", p, fc, fc_1);
    }

    fclose(fp);
    fclose(fp_1);
    fclose(fp_2);

    return 0;
}