变量意外更改值

我的代码有问题。调试时，我发现变量"de"在第一次循环 m 后更改了其值，并导致我的错误结果。"de"应该是一个恒定的直通处理。

感谢您帮助我更正代码！这是我的代码

#include <iostream>
#include <stdio.h>
#include <fstream>
#include <math.h>
using namespace std;
//Physical constants
#define ER 5.7599 //Rydberg energy in meV unit
#define h  658.2119 //Planck constant in meV.fs unit
#define pi 3.1415926 //pi number = 3.141592653589793238463
//Input constant
double Delta=30.0, T2=100.0, Chi0=0.4, emax=300.0, del_t=33.3333333, tmax=500.0; //time in fs unit, energy in meV unit
const int M = 300, N = 100;
//Energy step, time step
double de = emax/N;
double dt = (tmax+3*del_t)/M;
//
double y[N][3]; //y[][0] is f, y[][1] is Rep, y[][2] is Imp
double k[N][3], F[N][3], Y[N][3];
int n, i;
//Define used funcs
//init func for arrays [][3]
void initArray(double y[][3]) {
    for (i = 0; i < N; i++){
        for (int j = 0; j <= 2; j++) {y[i][j] = 0.0;}}
}
//func g
double g (int a, int b){
    double g_va;
    if (a == b) g_va = 0.0;
    else
        g_va = sqrt(ER)/pi/sqrt(double(a*de))*log(fabs((sqrt(double(a))+sqrt(double(b)))/(sqrt(double(a))-sqrt(double(b))))); //"log" in C++ is natural logarithm
    return g_va;}
//RHS(t, f, Rep, Imp)
void RHS(double t, double y[][3], double F[][3]){
    double Es, ReIp, ImIp;
    for (n=1;n<=N;n++){
        Es = ReIp = ImIp = 0.0;
        //Calculate some energy integrals
        for (int s=1;s<=N;s++){
            //for electron density
            Es += 2.0*de*g(n, s)*y[s-1][0];
            //polarization density ReIp and ImIp
            ReIp += de*g(n, s)*y[s-1][1];
            ImIp += de*g(n, s)*y[s-1][2];
        }
        //dipole energy
        double Epo = 0.5*h*sqrt(pi)*Chi0/del_t*exp(-t*t/del_t/del_t);
        //Rabi frequency
        double ReRabi = (Epo + ReIp)/h;
        double ImRabi = ImIp/h;
        //Re-norm energy for electron-hole
        double En = de*n - Es;
        F[n-1][0] = -2.0*(y[n-1][1]*ImRabi - y[n-1][2]*ReRabi);
        F[n-1][1] = (En-Delta)*y[n-1][2]/h - y[n-1][1]/T2 + (2.0*y[n-1][0] - 1)*ImRabi;
        F[n-1][2] = -(En-Delta)*y[n-1][1]/h - y[n-1][2]/T2 - (2.0*y[n-1][0] - 1)*ReRabi;
    }
}
//RK4 method
void RK4 (double t, double y[][3], void (*RHS)(double t, double y[][3], double F[][3])){
    RHS(t, y, F); //calculate k1
    for (n=1;n<=N;n++){
        for (i=0;i<=2;i++){
            k[n-1][i] = F[n-1][i];
            Y[n-1][i] = y[n-1][i] + 0.5*F[n-1][i];} //calculate y for k2
    }
    RHS(t + 0.5*dt, Y, F); //calculate k2
    for (n=1;n<=N;n++){
        for (i=0;i<=2;i++){
            k[n-1][i] += 2.0*F[n-1][i];} //append 2*k2 to k1
            Y[n-1][i] = y[n-1][i] + 0.5*F[n-1][i]; //calculate y for k3
    }
    RHS(t + 0.5*dt, Y, F); //calculate k3
    for (n=1;n<=N;n++){
        for (i=0;i<=2;i++){
            k[n-1][i] += 2.0*F[n-1][i];} //append 2*k3 to k1
            Y[n-1][i] = y[n-1][i] + F[n-1][i]; //calculate y for k4
    }
    RHS(t + dt, Y, F); //calculate k4
    for (n=1;n<=N;n++){
        for (i=0;i<=2;i++){
            k[n-1][i] += F[n-1][i]; //append k4 to k1
            y[n-1][i] += (dt/6)*k[n-1][i];} //calculate y(t+dt)
    }
}
int main (){
    double t = -3.0*del_t;
    double Ne, p, ReP, ImP, P;
    initArray (y); initArray(F); initArray(k); initArray(Y);
    FILE *output1; FILE *output2;
    output1=fopen("output1.dat", "w"); output2=fopen("output2.dat", "w");
    //prepare to print data
    fprintf(output1,"%10.8stt%10.8stt%10.8stt%10.8stt%10.8stt%10.8sn",
                    "time(fs)","erg(meV)","f","Rep","Imp","p");
    fprintf(output2,"%10.8stt%10.8stt%10.8sn","time(fs)","Edensity","Pdensity");
    for (int m=0;m<=M;m++){
    Ne = ReP = ImP = 0.0;
    for (n=1;n<=N;n++){
        //Calculate density
        p = sqrt(y[n-1][1]*y[n-1][1] + y[n-1][2]*y[n-1][2]);
        Ne += de*sqrt(n*de)*y[n-1][0]; //total e or h density
        ReP += de*sqrt(n*de)*y[n-1][1];
        ImP += de*sqrt(n*de)*y[n-1][2];
        //print y[], p
        fprintf(output1,"%10.5ftt%10.5ftt%10.5ftt%10.5ftt%10.5ftt%10.5fn", //maximum 10 digits before dot and 5 digits after dot, right-ident
                        t,de*n,y[n-1][0],y[n-1][1],y[n-1][2],p);
        //print to screen
        cout << t << "tt" << de*n << "tt" << y[n-1][0] << "tt" << p << "n";
        }
    //print out Ne and P
    P = sqrt(ReP*ReP + ImP*ImP); //total polarisation density
    fprintf(output2,"%10.5ftt%10.5ftt%10.5fn",t,Ne,P);
    //Apply RK4 method
    if (m == M) goto stop;
    else
    RK4 (t, y, &RHS); //after RK4 process, array y is updated and will become new initial value for next loop of t
    t = t + dt;
    }
    stop:
    //close output files
    fclose(output1); fclose(output2);
    cout << "completion!" << "n";
    return 0;
}