Теряется точность в сравнении с exel

#include <iostream>
#include <cmath>
#include <conio.h>
using namespace std;
 
int main(){
setlocale(LC_CTYPE,"rus");
int A = 235;//Массовое число элемента.
long double sigma_a,sigma_f,sigma_s;// Микросечения захвата, деления и расеяния.
long double sigma_tr;// Транспортное микросечение.
long double sigma_a_ock =1.21 * 1e-28;//
long double v;// Среднее число нейтронов.
long double  t;
long double M_dv = 40;// Масса ДВ.
long double M_shkr;// Критическая масса шара из ДВ..
long double k_bes;//
long double p;//
long double R_kr;// Критический радиус шара из ДВ.
long double Pi = 3.14159;//Число Пи.
long double Ed;//
switch(A){
    case 233:sigma_a = 1.32e-28,sigma_f = 1.27e-28,sigma_s = 8.7e-28;
    sigma_tr = sigma_s + sigma_a;
    v = 2.51;
    M_shkr = 16.45;
    k_bes =2.760000000000000000000000000;
    p = 18700;
    R_kr = 0.059000000000000000000000000;
    Ed = 191;
    break;
 
    case 235:sigma_a = 1.32 * 1e-28,sigma_f = 1.27 * 1e-28,sigma_s = 8.7 * 1e-28;
    sigma_tr = sigma_s + sigma_a;
    v=2.44;
    M_shkr = 35.888;
    k_bes=2.64000000000000000000000000000;
    p = 18700;
    R_kr = 0.0770000000000000000000000000;
    Ed = 214.;
    break;
    
    case 239:sigma_a = 2.02e-28,sigma_f = 1.98e-28,sigma_s = 7.9e-28;
    sigma_tr = sigma_s + sigma_a;
    v = 2.51;
    M_shkr = 10.146;
    k_bes = 3.09;
    p = 19740;
    R_kr = 0.05;
    Ed = 201;
    break;
    
    default:cout<<" Вы ошиблсь. Введите 233,235 или 239 в соответствии с вашим ДВ.";
    
}
long double NA = 6.0225 * 1e+26;  
long double const N0 = NA * M_dv / A;
long double R_sh = pow ((3 * M_dv) / (4 * Pi * p) ,1./3.);/* Радиус шара,м. */
long double R_shsm = R_sh * 100;// Радиус шара,см. 
long double Tau_sjatiya = R_shsm * 1e-6;// Время сжатия ДВ.
long double delta_v_max = 1.99;
long double delta_kp = sqrt(M_shkr / M_dv);
long double t_kp = (delta_kp - 1) * Tau_sjatiya / (delta_v_max - 1);// Момент перехода шара в надкритическое состояние.
long double t0 = 0.77 * Tau_sjatiya;// Начало инициировния.
long double delta_t = (Tau_sjatiya - t0) / 5;
long double SIGMA_a = (NA/A) * p * sigma_a,SIGMA_f = (NA/A) * p * sigma_f,SIGMA_s = (NA/A) * p * sigma_s,SIGMA_tr = SIGMA_a + SIGMA_s;
long double Tau_c = 1e-7; 
long double  k,ny,g,d;
long double delta_v_sj,delta_v_razl,delta_v,PER1;
long double Np[20],NP;
int i=0;
 
/*cout<<"sigma_tr= "<<sigma_tr<<endl;
cout<<"N0= "<<N0<<endl;
cout<<"delta_v_max= "<<delta_v_max<<endl;
cout<<"sigma_a_ock= "<<sigma_a_ock<<endl;
cout<<"Tau_sjatiya= "<<Tau_sjatiya<<endl;
cout<<"t0= "<<t0<<endl;
cout<<"delta_t= "<<delta_t<<endl;
cout<<"SIGMA_a= "<<SIGMA_a<<endl;
cout<<"SIGMA_f= "<<SIGMA_f<<endl;
cout<<"SIGMA_s= "<<SIGMA_s<<endl;
cout<<"SIGMA_tr= "<<SIGMA_tr<<endl;
cout<<"Tau_c= "<<Tau_c<<endl;*/
// РАСЧЁТ
// 1 цикл
    Np[i] = 209.;
    t = t0;
    delta_v_sj = (1 + (delta_v_max - 1) * t/Tau_sjatiya);
    delta_v_razl = 0;
    delta_v = delta_v_sj - delta_v_razl;
    
    d = 0.71 / (v * SIGMA_f - SIGMA_a + SIGMA_tr);
    
    k = 1 + (k_bes - 1) * (1 - pow( (R_kr + d) / (R_sh * pow(delta_v , 2./3.)+d ),2 ));
    int j = 1;
    Np[j] = Np[i] / (k - 1) * ( pow(k, (t / Tau_c + 1)) - 1);
/*cout<<"Np["<<j<<"] = "<< Np[j] <<endl;
cout<<"N0= "<<N0<<endl;*/
 
    ny = Np[j]/N0;
 
    g = 20 * ny * M_dv;
 
/*cout<<" Np["<<i<<"]= "<<Np[i]<<endl;
cout<<"k = "<<k <<endl;
cout<<"t / Tau_c + 1 = "<<t / Tau_c + 1<<endl;
cout<<" (k - 1) = "<<(k - 1)<<endl;
cout<<" pow(k, (t / Tau_c + 1)) = "<<pow(k, (t / Tau_c + 1))<<endl;
cout<<" ( pow(k, (t / Tau_c + 1)) - 1) = "<<( pow(k, (t / Tau_c + 1)) - 1)<<endl;*/
    cout.precision(4);
    cout<<" 1 цикл\n "<<endl;
    cout<<"Np[0] "<<Np[i]<<endl;
    cout<<" t = "<<t<<endl;
    cout<<" delta_v= "<<delta_v<<endl;
    cout<<"d  = "<<d <<endl;
    cout<<" k_bes = "<<k_bes<<endl;
    cout<<" R_kr = "<<R_kr<<endl;
    cout<<" d = "<<d<<endl;
    cout<<" R_sh = "<<R_sh<<endl;
    cout<<" delta_v = "<<delta_v<<endl;
    cout<<"k = "<<k <<endl;
    cout<<" Np"<<j<<"= "<<Np[j]<<endl;
    cout<<" ny = "<<ny*100<<"%"<<endl;
    cout<<"g = "<<g <<endl;
    cout<<"\n";
    
    do{
        i += 1;
        t = t0 + i * delta_t;
        if (t <= Tau_sjatiya) delta_v_sj = (1 + (delta_v_max-1) * t/Tau_sjatiya);
        else delta_v_sj = delta_v_max + (delta_v_max - 1) / (0.4 * Tau_sjatiya) * (t - Tau_sjatiya);
        delta_v_razl = pow( delta_v_sj*t / R_sh * sqrt(Ed * Np[i] / (2 * N0)) + 1 ,3.);
        delta_v = delta_v_sj - delta_v_razl;
        SIGMA_f = (delta_v * N0 - Np[i]) * sigma_f;
        SIGMA_a = (delta_v * N0 - Np[i]) * sigma_a + 2 * Np[i] * sigma_a_ock;
        SIGMA_tr = delta_v * N0 * sigma_tr;
        
        d = 0.71 / (v * SIGMA_f -SIGMA_a + SIGMA_tr);
        
        k = 1 + (k_bes - 1) * (1 - pow( (R_kr + d) / (R_sh * pow(delta_v * delta_v , 1./3.)+d ),2 ));
        
        Tau_c = 1 / (1.97 * SIGMA_a * 1e+7);
        j += 1;
        NP = Np[i] / (k - 1) * ( pow(k, t / Tau_c + 1) - 1);
        Np[j] = Np[i]+ NP;    
     
        ny = Np[j] / N0;
     
        g = 20 * ny * M_dv;
        cout<<" " <<j << "  цикл\n"<<endl;
        cout<<" t = "<<t<<endl;
        cout<<" Np0 = "<<Np[i]<<endl;
        cout<<" delta_v_sj = "<<delta_v_sj<<endl;
        cout<<" delta_v_razl = "<<delta_v_razl<<endl;
        cout<<" delta_v = "<<delta_v<<endl;
        cout<<" SIGMA_f = "<<SIGMA_f<<endl;
        cout<<" SIGMA_a = "<<SIGMA_a<<endl;
        cout<<" SIGMA_tr = "<<SIGMA_tr<<endl;
        cout<<" d = "<<d<<endl;
        cout<<" k = "<<k<<endl;
        cout<<" Tau_c = "<<Tau_c<<endl;
        cout<<" NP= "<<NP<<endl;
        cout<<" Tau_c = "<<Tau_c<<endl;
        cout<<" Np"<<j<<"= "<<Tau_c<<endl;
        cout<<" ny = "<<ny*100<<"%"<<endl;
        cout<<" g = "<<g<<endl;
        cout << "\n";
    } while (ny < 0.05e+00 && SIGMA_f > 0 && SIGMA_s > 0 && SIGMA_a > 0);
 
getch();
return 0;   
}