C ++ для цикла останавливается до достижения конца и перезапускается

Question

Я совал свой нос в работу с массивами в c ++, и я писал код решателя 1D Euler, который я написал в matlab, и преобразовывал его в c ++ в качестве практического упражнения.

Эта проблемазаключается в том, что цикл for должен работать до тех пор, пока счетчик не достигнет значения N_cells-1, но независимо от того, насколько высоко я установил число, он всегда будет равен 57, затем перезапускается с 2 и продолжает делать это, пока я не нажму на экран вывода.Я также запустил код с номером N_cells меньше 57, и я получил код ошибки, который я включил ниже.

Я довольно плохо знаком с массивами и заголовочными файлами в c ++, так что я уверен, что это что-топросто, но я просто не могу его найти.Я знаю, что это связано с массивом fqL, но я не знаю, что.

Ошибка при использовании числа <57: </strong>

#include "stdafx.h"
#include "Flux.h"
#include <iostream>
#include <chrono>
using namespace std;

void Flux(double * q, double y, double R, int N_cells,double * Flux)
{
    double qL[3];
    double qR[3];

    for (int i = 0; i < N_cells - 1; i++) {
        //Initialize left and right sides
        //-------------------

        qL[0] = q[0, i];
        qL[1] = q[1, i];
        qL[2] = q[2, i];

        qR[0] = q[0, i + 1];
        qR[1] = q[1, i + 1];
        qR[2] = q[2, i + 1];
        //-------------------

        //Calculate Left side Parameters
        //-------------------
        double PL;
        //double fqL[3];
        double cL2;
        double HL;
        double uL;
        PL = (y - 1)*(qL[2] - 0.5 / qL[0] * (qL[1] * qL[1]));
        double fqL[3] = { qL[1],
            (3 - y)*0.5*(qL[1] * qL[1]) / qL[0] + (y - 1)*qL[2],
            y*qL[1] * qL[2] / qL[0] - (y - 1)*0.5*(qL[1] * qL[1] * qL[1]) / (qL[0] * qL[0]) };
        cL2 = y * (y - 1)*(qL[2] / qL[0] - 0.5*(qL[1] / qL[0])*(qL[1] / qL[0]));
        HL = 0.5*(qL[1] / qL[0])*(qL[1] / qL[0]) + cL2 / (y - 1);
        uL = qL[1] / qL[0];
        //Calculate Right side Parameters
        //-------------------
        double PR;
        //double fqR[3];
        double cR2;
        double HR;
        double uR;
        PR = (y - 1)*(qR[2] - 0.5 / qR[0] * (qR[1] * qR[1]));
        double fqR[3] = { qR[1],
            (3 - y)*0.5*(qR[1] * qR[1]) / qR[0] + (y - 1)*qR[2],
            y*qR[1] * qR[2] / qR[0] - (y - 1)*0.5*(qR[1] * qR[1] * qR[1]) / (qR[0] * qR[0]) };
        cR2 = y * (y - 1)*(qR[2] / qR[0] - 0.5*(qR[1] / qR[0])*(qR[1] / qR[0]));
        HR = 0.5*(qR[1] / qR[0])*(qR[1] / qR[0]) + cR2 / (y - 1);
        uR = qR[1] / qR[0];
        //-------------------

        //Calculate Roe's Variables
        //-------------------------------- -
        double u;
        double H;
        double c;
        double rho;
        u = (sqrt(qL[1])*qL[2] / qL[1] + sqrt(qR[1])*qR[2] / qR[1]) / (sqrt(qL[1]) + sqrt(qR[1]));
        H = (sqrt(qL[1])*HL + sqrt(qR[1])*HR) / (sqrt(qL[1]) + sqrt(qR[1]));
        c = sqrt((y - 1)*(H - 0.5*u *u));
        rho = sqrt(qL[1] * qR[1]);
        //-------------------------------- -


        //-------------------------------- -
        double g[3] = { u - c, u, u + c };

        double v[3][3] = { {1, u - c, H - u * c},
        {1, u, 0.5*u*u},
        {1, u + c, H + u * c } };

        double a[3] = { 0.5 / (c*c)*((PR - PL) - c * rho*(uR - uL)),
            (qR[0] - qL[0]) - 1 * (PR - PL) / (c*c),
            0.5 / (c*c)*((PR - PL) + c * rho*(uR - uL)) };

        double SUM[3];

            SUM[0] = g[0] * a[0] * v[0][0] + g[1] * a[1] * v[1][0] + g[2] * a[2] * v[2][0];
            SUM[1] = g[0] * a[0] * v[0][1] + g[1] * a[1] * v[1][1] + g[2] * a[2] * v[2][1];
            SUM[2] = g[0] * a[0] * v[0][2] + g[1] * a[1] * v[1][2] + g[2] * a[2] * v[2][2];

        double Flux[3];
            Flux[0,i] = 0.5*(fqL[0] + fqR[0]) - 0.5*SUM[0];
            Flux[1,i] = 0.5*(fqL[1] + fqR[1]) - 0.5*SUM[1];
            Flux[2,i] = 0.5*(fqL[2] + fqR[2]) - 0.5*SUM[2];

            std::cout << i << endl;


    }



}