<commit_before>#include "burnupcalc.h"
#include <boost/math/special_functions/bessel.hpp>


using namespace std;

typedef unsigned long long timestamp_t;

static timestamp_t
get_timestamp ()
    {
      struct timeval now;
      gettimeofday (&now, NULL);
      return  now.tv_usec + (timestamp_t)now.tv_sec * 1000000;
    }

double intpol(double y0, double y1, double x0, double x1, double x) {
    // linear interpolation function
    double y = y0 + (y1 - y0)*(x - x0)/(x1 - x0);
    return y;
}



fuelBundle regionCollapse(fuelBundle fuel){
///add micro region flux effects
//struct effects accounted here
    //cout << "Begin regionCollapse" << endl;
    for(int i = 0; i < fuel.batch.size(); i++){
        //for(int j = 0; j < fuel.batch[i].iso.size(); j ++){
        fuel.batch[i].collapsed_iso = FuelBuilder(fuel.batch[i].iso);
        //}
        //builds total BU from BUd
        fuel.batch[i].collapsed_iso.BU.push_back(fuel.batch[i].collapsed_iso.BUd[0]);

        for(int j = 1; j < fuel.batch[i].collapsed_iso.BUd.size(); j++){
        //cout << "    1.5tst" << i+1 << "  " << j << endl;
            fuel.batch[i].collapsed_iso.BU.push_back(fuel.batch[i].collapsed_iso.BU[j-1]+fuel.batch[i].collapsed_iso.BUd[j]);
        }
        //test to see if the prod/dest vectors are the same length
        if(fuel.batch[i].collapsed_iso.neutron_prod.size() != fuel.batch[i].collapsed_iso.neutron_dest.size()){
            cout << "Error. Neutron production/destruction rate vector length mismatch." << endl;
        }
        if(fuel.batch[i].collapsed_iso.neutron_prod[0] == 0 || fuel.batch[i].collapsed_iso.neutron_dest[0] == 0){
            fuel.batch[i].collapsed_iso.neutron_prod.erase(fuel.batch[i].collapsed_iso.neutron_prod.begin());
            fuel.batch[i].collapsed_iso.neutron_dest.erase(fuel.batch[i].collapsed_iso.neutron_dest.begin());
            fuel.batch[i].collapsed_iso.BUd.erase(fuel.batch[i].collapsed_iso.BUd.begin());
            fuel.batch[i].collapsed_iso.BU.erase(fuel.batch[i].collapsed_iso.BU.begin());
            fuel.batch[i].collapsed_iso.fluence.erase(fuel.batch[i].collapsed_iso.fluence.begin());
            for(int j = 0; j < fuel.batch[i].collapsed_iso.iso_vector.size(); j++){
                fuel.batch[i].collapsed_iso.iso_vector[j].mass.erase(fuel.batch[i].collapsed_iso.iso_vector[j].mass.begin());
            }
        }
    }
    return fuel;
};


fuelBundle fast_region_collapse(fuelBundle fuel){
///add micro region flux effects
//struct effects accounted here
    //cout << "Begin regionCollapse" << endl;
    for(int i = 0; i < fuel.batch.size(); i++){
        //for(int j = 0; j < fuel.batch[i].iso.size(); j ++){
        fuel.batch[i].collapsed_iso = BurnupBuilder(fuel.batch[i].iso);
        //}
        //builds total BU from BUd
        fuel.batch[i].collapsed_iso.BU.push_back(fuel.batch[i].collapsed_iso.BUd[0]);

        for(int j = 1; j < fuel.batch[i].collapsed_iso.BUd.size(); j++){
        //cout << "    1.5tst" << i+1 << "  " << j << endl;
            fuel.batch[i].collapsed_iso.BU.push_back(fuel.batch[i].collapsed_iso.BU[j-1]+fuel.batch[i].collapsed_iso.BUd[j]);
        }
        //test to see if the prod/dest vectors are the same length
        if(fuel.batch[i].collapsed_iso.neutron_prod.size() != fuel.batch[i].collapsed_iso.neutron_dest.size()){
            cout << "Error. Neutron production/destruction rate vector length mismatch." << endl;
        }
        if(fuel.batch[i].collapsed_iso.neutron_prod[0] == 0 || fuel.batch[i].collapsed_iso.neutron_dest[0] == 0){
            fuel.batch[i].collapsed_iso.neutron_prod.erase(fuel.batch[i].collapsed_iso.neutron_prod.begin());
            fuel.batch[i].collapsed_iso.neutron_dest.erase(fuel.batch[i].collapsed_iso.neutron_dest.begin());
            fuel.batch[i].collapsed_iso.BUd.erase(fuel.batch[i].collapsed_iso.BUd.begin());
            fuel.batch[i].collapsed_iso.BU.erase(fuel.batch[i].collapsed_iso.BU.begin());
            fuel.batch[i].collapsed_iso.fluence.erase(fuel.batch[i].collapsed_iso.fluence.begin());
        }
    }
    return fuel;
};


map<int, double> tomass (int ti, double fluence, isoInformation isoinfo) {
    map<int, double> out = map<int, double>();
    double mass_i;
    int name_i;
    int nucid;
    for (int i = 0; i < isoinfo.iso_vector.size(); i++){
        name_i = isoinfo.iso_vector[i].name;
        nucid = pyne::nucname::zzaaam(name_i);
        mass_i = intpol(isoinfo.iso_vector[i].mass[ti-1],
                        isoinfo.iso_vector[i].mass[ti],
                        isoinfo.fluence[ti-1],
                        isoinfo.fluence[ti],
                        fluence);
        out[nucid] = mass_i;
    }
    return out;
}


fuelBundle phicalc_simple(fuelBundle core){
//updates the rflux of each batch in core.batch
//assumes the flux of a batch is proportional to the inverse neutron prod rate
    double maxphi = 0;
    //finds the inverse of neutron production at the batch_fluence
    //stores it in core.batch.rflux
    for(int i = 0; i < core.batch.size(); i++){
        int ii;
        if(core.batch[i].collapsed_iso.fluence.back() < core.batch[i].Fg){
            //cout << endl << "Maximum fluence error! Batch fluence exceeded max library fluence. (method2)" << endl;
            //cout << " Values on max fluence will be used. Do not trust results." << endl;
            for(int j = 0; j < core.batch.size(); j++){
                core.batch[j].rflux = 1;
            }
            return core;
        } else {
            for(ii = 0; core.batch[i].collapsed_iso.fluence[ii] < core.batch[i].Fg; ii++){}
        }
        if(ii == 0){
            core.batch[i].rflux = 1/core.batch[i].collapsed_iso.neutron_prod[0];
            if(maxphi < core.batch[i].rflux){
                maxphi = core.batch[i].rflux;
            }
        } else{
            core.batch[i].rflux =
                1/intpol(core.batch[i].collapsed_iso.neutron_prod[ii-1],core.batch[i].collapsed_iso.neutron_prod[ii], core.batch[i].collapsed_iso.fluence[ii-1], core.batch[i].collapsed_iso.fluence[ii], core.batch[i].Fg+1);
        }
        if(maxphi < core.batch[i].rflux){
            maxphi = core.batch[i].rflux;
        }
    }
    //normalizes all the flux values
    for(int i = 0; i < core.batch.size(); i++){
        core.batch[i].rflux /= maxphi;
        //cout << core.batch[i].rflux << endl;
    }
    //cout << endl;
    return core;
}

fuelBundle phicalc_eqpow(fuelBundle core){
    //updates the rflux of each batch in core.batch
    //assumes the flux of a batch is proportional to the inverse neutron prod rate
    double maxphi = 0;
    //finds the inverse of neutron production at the batch_fluence
    //stores it in core.batch.rflux
    for(int i = 0; i < core.batch.size(); i++){
        int ii;


        if(core.batch[i].collapsed_iso.fluence.back() < core.batch[i].Fg){
            //cout << endl << "Maximum fluence error! Batch fluence exceeded max library fluence. (method2)" << endl;
            //cout << "  Values on max fluence will be used. Do not trust results." << endl;
            for(int j = 0; j < core.batch.size(); j++){
                core.batch[j].rflux = 1;
            }
            return core;

        } else {
            for(ii = 0; core.batch[i].collapsed_iso.fluence[ii] < core.batch[i].Fg; ii++){}
        }

        if(ii == 0){
            //core.batch[i].rflux = 1/core.batch[i].collapsed_iso.neutron_prod[0];
            core.batch[i].rflux = (core.batch[i].collapsed_iso.fluence[ii+1]-core.batch[i].collapsed_iso.fluence[ii])
                /(core.batch[i].collapsed_iso.BU[ii+1]-core.batch[i].collapsed_iso.BU[ii]);

        } else{

            core.batch[i].rflux = (core.batch[i].collapsed_iso.fluence[ii]-core.batch[i].collapsed_iso.fluence[ii-1])
                /(core.batch[i].collapsed_iso.BU[ii]-core.batch[i].collapsed_iso.BU[ii-1]);
        }

        if(maxphi < core.batch[i].rflux){
            maxphi = core.batch[i].rflux;
        }
    }

    //normalizes all the flux values
    //cout << "EqPow Flux: ";
    for(int i = 0; i < core.batch.size(); i++){
        core.batch[i].rflux /= maxphi;
        //cout << core.batch[i].rflux << "   ";
    }
    //cout << endl;

    return core;
}

double nusigf_finder(batch_info batch){
    double Nusig_f;
    int ii;
    double prod = 0;
    int A = 235;

    //find the index to interpolate on
    for(ii = 0; batch.collapsed_iso.fluence[ii] < batch.Fg; ii++);
    if(batch.collapsed_iso.fluence.back() < batch.Fg){
        ii = batch.collapsed_iso.fluence.size() - 1;
    }

    prod = intpol(batch.collapsed_iso.neutron_prod[ii-1], batch.collapsed_iso.neutron_prod[ii], batch.collapsed_iso.fluence[ii-1], batch.collapsed_iso.fluence[ii], batch.Fg);


    if(ii == 0){
        prod = batch.collapsed_iso.neutron_prod[0];
    }

    //Nusig_f[cm-1] = prod rate / Avagadros number * mass * barn
    //Nusig_f = prod/(6.0221413E+23)*A*1E24/1000;
    Nusig_f = prod*0.01097;

    return Nusig_f;
}

double siga_finder(batch_info batch){
//returns the macroscopic cross section at fluence Fg
    double sig_a;
    int ii;
    double dest = 0;
    int A = 235;

    //find the index to interpolate on
    for(ii = 0; batch.collapsed_iso.fluence[ii] < batch.Fg; ii++);
    if(batch.collapsed_iso.fluence.back() < batch.Fg){
        ii = batch.collapsed_iso.fluence.size() - 1;
    }

    if(ii == 0){
        dest = batch.collapsed_iso.neutron_dest[0];
    }


    dest = intpol(batch.collapsed_iso.neutron_dest[ii-1], batch.collapsed_iso.neutron_dest[ii], batch.collapsed_iso.fluence[ii-1], batch.collapsed_iso.fluence[ii], batch.Fg);

    //sig_a[barn] = dest rate [cm2/kg] / Avagadros number * mass * barn
    //sig_a = dest/(6.0221413E+23)*A*1E24/1000;
    sig_a = dest*0.01097;

    return sig_a;
}


////////////////////////////////////////////////////////////////////////////////////////////
fuelBundle phicalc_cylindrical(fuelBundle core){
//cout << "cylindrical" << endl<<endl<<endl<<endl<<endl;

//the word 'region' is interchanged with the word 'batch' in comments, water layer is the outermost region

    int region = core.batch.size(); // number of regions (batches) not counting the outer region
    double delta = core.cylindrical_delta;
    double R[region+1]; //radial thickness of each region
    int N[region+1]; //number of mesh points in each region
    int NC[region+1]; //cumulative N
    int NTotal; //total number of mesh points
    double dd2[region+1]; // D/delta^2
    double Sigma_a[region+1]; //mac. abs. cs of each region
    double NuSigma_f[region+1]; //nu sigma f of each region
    double Sigma_tr[region+1]; //mac. transport cs of each region
    double D[region+1]; //diff coef. for each region
    double LSquared[region+1];
    double k = 1;
    double k_prev = 0.9;
    double prod, prod_prev;
    double flux[region+1], maxflux=0;
    double sum = 0;

    //set the radial thickness of each region
    R[0] = sqrt(core.fuel_area/region/3.141592);
    for(int i = 1; i < region; i++){
        R[i] = sqrt(core.fuel_area/region/3.141592*(i+1));
    }
    R[region] = R[region-1] + core.mod_thickness; //this is the moderator region


/*
    Sigma_a[0] = 0.0230;
    Sigma_a[1] = 0.0246;
    Sigma_a[2] = 0.0324;
    NuSigma_f[0] = 0.0184;
    NuSigma_f[1] = 0.0217;
    NuSigma_f[2] = 0.0382;*/
    //assign fuel cross sections
    for(int i = 0; i < region; i++){

        NuSigma_f[i] = nusigf_finder(core.batch[i]);

        Sigma_a[i] = siga_finder(core.batch[i]);

        Sigma_tr[i] = core.fuel_Sig_tr;
        D[i] = 1/(Sigma_tr[i]*3.);
        LSquared[i] = D[i]/Sigma_a[i];

    }

    //assign moderator cross sections
    Sigma_a[region] = 0.0066;// core.mod_Sig_a;
    Sigma_tr[region] = core.mod_Sig_tr;
    D[region] = 1/(Sigma_tr[region]*3.);
    LSquared[region] = D[region]/Sigma_a[region];
    NuSigma_f[region] = core.mod_Sig_f;

    //populate dd2
    for(int i = 0; i < region+1; i++){
        dd2[i] = D[i]/delta/delta;
    }

    /*
    cout << "Cross sections in cylindrical calc: " << endl << "Sig_a      Sig_f      R" << endl;
    for(int r = 0; r < region+1; r++){
        cout << Sigma_a[r] << "  " << NuSigma_f[r] << "  " << R[r] << endl;
    }*/

    //populate N, number of mesh points in each region
    N[0] = R[0]/delta;
    NC[0] = N[0];
    NTotal = N[0];

    for(int i = 1; i < region+1; i++){
        N[i] = ceil((R[i] - R[i-1]) / delta);
        NC[i] = NC[i-1] + N[i];
        NTotal += N[i];
    }
    NC[region] += 1;
    NTotal += 1;/*
    Eigen::MatrixXf A(NTotal, NTotal);
    Eigen::MatrixXf F(NTotal, 1);
    Eigen::MatrixXf phi(NTotal, 1);
    Eigen::MatrixXf phi_prev(NTotal, 1);
    Eigen::MatrixXf S(NTotal, 1);
    Eigen::MatrixXf S_prev(NTotal, 1);

    A.setZero();
    F.setZero();

    int jprev = 0;
    int j;

    int r = 0; //region index
    for(int i = 1; i < NTotal-1; i++){
        A(i, i-1) = (-1.)*dd2[r]*(2*i-1)/(2*i);
        A(i,i) = dd2[r]*2. + Sigma_a[r];
        A(i, i+1) = (-1.)*dd2[r]*(2*i+1)/(2*i);
        if(i == NC[r]){
            r += 1;
        }
        if(core.batch[r].rflux < 1.1 && core.batch[r].rflux > 0){
            phi_prev(i) = core.batch[r].rflux; //uses last runs results if available
        } else{
            phi_prev(i) = 1;
        }

    }
    A(0,0) = 1;
    A(0,1) = -1;
    A(NTotal-1,NTotal-1) = 1;
    phi_prev(NTotal-1) = 0;

    for(r = 0; r < region; r++){
        A(NC[r],NC[r]-1) = D[r];
        A(NC[r],NC[r]) = -D[r]-D[r+1];
        A(NC[r],NC[r]+1) = D[r+1];
    }

    r = 0;
    for(int i = 1; i < NTotal; i++){
        if(i != NC[r]){
            F(i) = NuSigma_f[r];
        }
        if( i == NC[r]){
            r += 1;
        }
        S_prev(i) = F(i)*phi_prev(i);
    }


    for(int iter = 0; iter < 10; iter++){
        phi = A.colPivHouseholderQr().solve(S_prev)/k_prev;


        if(!phi.allFinite()){
            phi = phi_prev;
        }

        S = F.array() * phi.array();

        prod = (0.25)*3.141592*NuSigma_f[0]*phi(0)*delta*delta;
        prod_prev = (0.25)*3.141592*NuSigma_f[0]*phi_prev(0)*delta*delta;
        r = 0;

        for(int i = 0; i < NTotal; i++){
            if(i == NC[r]){
                prod += 3.141592*NuSigma_f[r]*phi(i)*(i-0.25)*delta*delta;
                prod_prev += 3.141592*NuSigma_f[r]*phi_prev(i)*(i-0.25)*delta*delta;
                r += 1;
                prod += 3.141592*NuSigma_f[r]*phi(i)*(i-0.25)*delta*delta;
                prod_prev += 3.141592*NuSigma_f[r]*phi_prev(i)*(i-0.25)*delta*delta;
            } else {
                prod += 2.*3.141592*NuSigma_f[r]*phi(i)*i*delta*delta;
                prod_prev += 2.*3.141592*NuSigma_f[r]*phi_prev(i)*i*delta*delta;
            }
        }

        if(abs((k_prev-k)/k) < 0.01 && iter > 3){
            break;
        }
        //cout << "prod: " << prod << "  prod_prev: " << prod_prev << "  k: " << prod/prod_prev << endl;
        k = prod/prod_prev*k_prev;
        phi_prev = phi;
        k_prev = k;
        S_prev = S;


    }

    //NO NEED TO normalize phi
    //phi = phi.array()/phi.maxCoeff();

    //find area weighted average phi per batch
    r = 0;
    flux[0] = 0;
    for(int i = 0; i < NTotal; i++){
        flux[r] += phi(i)*(2*(i+1)-1);
        sum += (2*(i+1)-1);

        if(i == NC[r] || i == NTotal-1){
            flux[r] /= sum;
            sum = 0;
            r += 1;
            flux[r] = 0;
        }
    }
    for(r = 0; r < region+1; r++){
        if(flux[r] > maxflux){
            maxflux = flux[r];
        }
    }

    //cout << "flux: ";
    //cout << core.fuel_area << " " << delta << " | " << Sigma_a[0] << " " << Sigma_a[1] << " " << Sigma_a[2] << " | "
    //<< NuSigma_f[0] << " " << NuSigma_f[1] << " " << NuSigma_f[2] << " | ";
    //normalize the fluxes
    for(r = 0; r < region+1; r++){
        flux[r] /= maxflux;
        //cout << flux[r] << " ";
    }
    //ĺcout << endl;




    //NO NEED TO normalize phi
    //phi = phi.array()/phi.maxCoeff();

    //find area weighted average phi per batch
    r = 0;
    flux[0] = 0;
    for(int i = 0; i < NTotal; i++){
        flux[r] += phi(i)*(2*(i+1)-1);
        sum += (2*(i+1)-1);

        if(i == NC[r] || i == NTotal-1){
            flux[r] /= sum;
            sum = 0;
            r += 1;
            flux[r] = 0;
        }
    }
    for(r = 0; r < region+1; r++){
        if(flux[r] > maxflux){
            maxflux = flux[r];
        }
    }

    //cout << "flux: ";
    //cout << core.fuel_area << " " << delta << " | " << Sigma_a[0] << " " << Sigma_a[1] << " " << Sigma_a[2] << " | "
    //<< NuSigma_f[0] << " " << NuSigma_f[1] << " " << NuSigma_f[2] << " | ";
    //normalize the fluxes
    for(r = 0; r < region+1; r++){
        flux[r] /= maxflux;
        //cout << flux[r] << " ";
    }
    //cout << endl;


    for(int i = 0; i < core.batch.size(); i++){
        core.batch[i].rflux = flux[i];
    }
*/
    return core;
}
/////////////////////////////////////////////////////////////////////////////////////////////

double kcalc(fuelBundle core){
    //uses the fluence values in core.batch.Fg to calculate the core criticality
    //boost::timer t;
    int N = core.batch.size();
    double prod_tot = 0;
    double dest_tot = 0;
    double pnl = core.pnl;


    //finds the index j where fluence is just under target, then interpolates
    //  to find the prod and dest values for each batch
    int j;
    for(int i = 0; i < N; i++){

        if(core.batch[i].collapsed_iso.fluence.back() < core.batch[i].Fg){
            //cout << endl << "Maximum fluence error! Batch fluence exceeded max library fluence. (kcalc)" << endl;
            //cout << "  Values on max fluence will be used. Do not trust results." << endl;
            j = core.batch[i].collapsed_iso.fluence.size() - 1;
        } else {
            //find dicrete point j to interpolate on
            for(j = 0; core.batch[i].collapsed_iso.fluence[j] < core.batch[i].Fg; j++){}
        }

        //add the production rate of batch i to total production
        prod_tot += intpol(core.batch[i].collapsed_iso.neutron_prod[j-1], core.batch[i].collapsed_iso.neutron_prod[j], core.batch[i].collapsed_iso.fluence[j-1], core.batch[i].collapsed_iso.fluence[j], core.batch[i].Fg);

        //add structural material production of this batch, scaled up using disadvantage factor
        //prod_tot += core.struct_prod * core.batch[i].DA;

        //add the destruction rate of batch i to total destruction
        dest_tot += intpol(core.batch[i].collapsed_iso.neutron_dest[j-1], core.batch[i].collapsed_iso.neutron_dest[j], core.batch[i].collapsed_iso.fluence[j-1], core.batch[i].collapsed_iso.fluence[j], core.batch[i].Fg);

        //add structural material destruction of this batch, scaled up using disadvantage factor
        //dest_tot += core.struct_dest * core.batch[i].DA;
        //cout << "diff: " << prod_tot - dest_tot << endl;
    }
    if(pnl < 0 || pnl > 1.5){
        cout << endl << "Error in core nonleakage! Assumed new value: 1.000" << endl;
        pnl = 1.000;
        }

    //std::cout << "k calc time " << t.elapsed() << std::endl;
    return prod_tot * pnl / dest_tot;
}

double CR_finder(fuelBundle core){
    //
    //i is the batch number starting from zero
    //boost::timer t;
    double FP = 0, FP0 = 0, FP1 = 0;
    double fissile = 0, fissile0 = 0, fissile1 = 0;
    double ini_fissile = 0;
    double CR;
    int ii, ZZ;

    //cout << "\nCR_finder starting." << endl << "  Upper: " << core.CR_upper << "  lower: " << core.CR_lower << endl;


    for(int i = 0; i < core.batch.size(); i++){
        for(ii = 0; core.batch[i].collapsed_iso.fluence[ii] < core.batch[i].Fg; ii++){}
        if(ii == 0){
            ii = 1;
        }

        for(int j = 0; j < core.batch[i].collapsed_iso.iso_vector.size(); j++){
            //convert name to mass number
            ZZ = core.batch[i].collapsed_iso.iso_vector[j].name;
            ZZ = ZZ % 10000;
            ZZ /= 10;

            //add up the FP
            if(ZZ < core.CR_upper && ZZ > core.CR_lower){
                //cout << "    Batch" << i+1 << " ZZ: " << ZZ << endl;
                //interpolation will be done at the end
                FP0 += core.batch[i].collapsed_iso.iso_vector[j].mass[ii-1];
                FP1 += core.batch[i].collapsed_iso.iso_vector[j].mass[ii];
            }

            //add up fissiles
            for(int fis = 0; fis < core.CR_fissile.size(); fis++){
                if(core.batch[i].collapsed_iso.iso_vector[j].name == core.CR_fissile[fis]){
                    fissile0 += core.batch[i].collapsed_iso.iso_vector[j].mass[ii-1];
                    fissile1 += core.batch[i].collapsed_iso.iso_vector[j].mass[ii];

                    ini_fissile += core.batch[i].collapsed_iso.iso_vector[j].mass[0]; //mass at fluence zero
                }
            }


        }


        // recycling variable FP0 here to check greater than zero
        FP0 = intpol(FP0, FP1, core.batch[i].collapsed_iso.fluence[ii-1], core.batch[i].collapsed_iso.fluence[ii], core.batch[i].Fg);
        if(FP0 > 0){
            FP += FP0;
        }

        fissile += intpol(fissile0, fissile1, core.batch[i].collapsed_iso.fluence[ii-1], core.batch[i].collapsed_iso.fluence[ii], core.batch[i].Fg);

        FP0 = 0;
        FP1 = 0;
        fissile0 = 0;
        fissile1 = 0;
    }

    //cout << "  FP: " << FP << "  fiss: " << fissile << "  ini_fiss: " << ini_fissile << "      CR: " << (FP+fissile-ini_fissile)/FP << endl;

    if(FP > 0){
        CR = (FP+fissile-ini_fissile)/FP;
    } else {
        CR  = 0;
    }
    //std::cout << "CR Calc " << t.elapsed() << std::endl;
    return CR;
}


fuelBundle burnupcalc(fuelBundle core, int mode, int DA_mode, double delta) {
timestamp_t t0 = get_timestamp();
    //this function only uses the COLLAPSED_ISO of each BATCH in the structure CORE
    //all factors that contribute to a change in neutron prod/dest rates have to be factored
    //      before calling this function
    //cout << endl << "Burnupcalc" << endl;
    boost::timer t;

    int N = core.batch.size(); //number of batches
    double dt = delta*24*60*60; //days to [s]
    double kcore, kcore_prev;
    double y0, y1;
    double burnup = 0;

/*
    cout << "burnupcalc starting neutron prods: " << core.batch[0].collapsed_iso.neutron_prod[0] << " " << core.batch[1].collapsed_iso.neutron_prod[0] << " " << core.batch[2].collapsed_iso.neutron_prod[0]  << endl;
    cout << "   " << core.batch[0].collapsed_iso.BU[1] << " " << core.batch[1].collapsed_iso.BU[1] << " " << core.batch[2].collapsed_iso.BU[1]  << endl;
    cout << "   " << core.batch[0].collapsed_iso.fluence[1] << " " << core.batch[1].collapsed_iso.fluence[1] << " " << core.batch[2].collapsed_iso.fluence[1]  << endl;
    cout << "     " << dt << endl;
    cout << " " << core.batch[0].batch_fluence << " " << core.batch[1].batch_fluence << " " << core.batch[2].batch_fluence  << endl;
*/

    //assign the batch_fluence to Fg
    for(int i = 0; i < N; i++){
        core.batch[i].Fg = core.batch[i].batch_fluence;
    }

    kcore = 3.141592;

    //more forward in time until kcore drops under 1
    while(kcore > 1){
        kcore_prev = kcore;
        //if(kcore != 3.141592){cout << "kcore: " << kcore << endl;}

        //find the normalized relative flux of each batch
        if(mode == 1){
            //simplest case, all batches get the same flux
            for(int i = 0; i < N; i++){
                core.batch[i].rflux = 1;
            }
        }else if(mode == 2){
            //inverse-production flux calculation
            core = phicalc_simple(core);
        }else if(mode == 3){
            core = phicalc_cylindrical(core);
        }else if(mode == 0){
            core = phicalc_eqpow(core);
        }else{
            cout << endl << "Error in mode input for batch-level flux calculation." << endl;
            return core;
        }

        //disadvantage calculation
        if(DA_mode == 1){
            core = DA_calc(core);
        }

        //CR calculation
        //cout << CR_numerator(core, 5)/CR_denominator(core, 5) << " ";
        core.CR = CR_finder(core);
        //update fluences
        for(int i = 0; i < N; i++){
            //cout << "flux: " << core.batch[i].rflux << endl;
            core.batch[i].Fg += core.batch[i].rflux * core.base_flux * dt;

        }

        kcore = kcalc(core);

        /*if(core.CR_target != 0 && std::abs(core.CR - core.CR_target)/core.CR < 0.1){
            break;
        }*/
    }
    core.CR = CR_finder(core);
    //cout << " CR: " << core.CR << endl;


    //update core fluences and CR
    for(int i = 0; i < N; i++){
        core.batch[i].batch_fluence = intpol(core.batch[i].Fg - (core.batch[i].rflux * core.base_flux * dt), core.batch[i].Fg, kcore_prev, kcore, 1);
        //cout << "batch " << i+1 << " CR: " << core.batch[i].CR << endl;
        //cout << "  fluence end of burnupcalc: " << core.batch[i].batch_fluence << endl;

    }<commit_msg>CR only calculated at end of cycle<commit_after>#include "burnupcalc.h"
#include <boost/math/special_functions/bessel.hpp>


using namespace std;

typedef unsigned long long timestamp_t;

static timestamp_t
get_timestamp ()
    {
      struct timeval now;
      gettimeofday (&now, NULL);
      return  now.tv_usec + (timestamp_t)now.tv_sec * 1000000;
    }

double intpol(double y0, double y1, double x0, double x1, double x) {
    // linear interpolation function
    double y = y0 + (y1 - y0)*(x - x0)/(x1 - x0);
    return y;
}



fuelBundle regionCollapse(fuelBundle fuel){
///add micro region flux effects
//struct effects accounted here
    //cout << "Begin regionCollapse" << endl;
    for(int i = 0; i < fuel.batch.size(); i++){
        //for(int j = 0; j < fuel.batch[i].iso.size(); j ++){
        fuel.batch[i].collapsed_iso = FuelBuilder(fuel.batch[i].iso);
        //}
        //builds total BU from BUd
        fuel.batch[i].collapsed_iso.BU.push_back(fuel.batch[i].collapsed_iso.BUd[0]);

        for(int j = 1; j < fuel.batch[i].collapsed_iso.BUd.size(); j++){
        //cout << "    1.5tst" << i+1 << "  " << j << endl;
            fuel.batch[i].collapsed_iso.BU.push_back(fuel.batch[i].collapsed_iso.BU[j-1]+fuel.batch[i].collapsed_iso.BUd[j]);
        }
        //test to see if the prod/dest vectors are the same length
        if(fuel.batch[i].collapsed_iso.neutron_prod.size() != fuel.batch[i].collapsed_iso.neutron_dest.size()){
            cout << "Error. Neutron production/destruction rate vector length mismatch." << endl;
        }
        if(fuel.batch[i].collapsed_iso.neutron_prod[0] == 0 || fuel.batch[i].collapsed_iso.neutron_dest[0] == 0){
            fuel.batch[i].collapsed_iso.neutron_prod.erase(fuel.batch[i].collapsed_iso.neutron_prod.begin());
            fuel.batch[i].collapsed_iso.neutron_dest.erase(fuel.batch[i].collapsed_iso.neutron_dest.begin());
            fuel.batch[i].collapsed_iso.BUd.erase(fuel.batch[i].collapsed_iso.BUd.begin());
            fuel.batch[i].collapsed_iso.BU.erase(fuel.batch[i].collapsed_iso.BU.begin());
            fuel.batch[i].collapsed_iso.fluence.erase(fuel.batch[i].collapsed_iso.fluence.begin());
            for(int j = 0; j < fuel.batch[i].collapsed_iso.iso_vector.size(); j++){
                fuel.batch[i].collapsed_iso.iso_vector[j].mass.erase(fuel.batch[i].collapsed_iso.iso_vector[j].mass.begin());
            }
        }
    }
    return fuel;
};


fuelBundle fast_region_collapse(fuelBundle fuel){
///add micro region flux effects
//struct effects accounted here
    //cout << "Begin regionCollapse" << endl;
    for(int i = 0; i < fuel.batch.size(); i++){
        //for(int j = 0; j < fuel.batch[i].iso.size(); j ++){
        fuel.batch[i].collapsed_iso = BurnupBuilder(fuel.batch[i].iso);
        //}
        //builds total BU from BUd
        fuel.batch[i].collapsed_iso.BU.push_back(fuel.batch[i].collapsed_iso.BUd[0]);

        for(int j = 1; j < fuel.batch[i].collapsed_iso.BUd.size(); j++){
        //cout << "    1.5tst" << i+1 << "  " << j << endl;
            fuel.batch[i].collapsed_iso.BU.push_back(fuel.batch[i].collapsed_iso.BU[j-1]+fuel.batch[i].collapsed_iso.BUd[j]);
        }
        //test to see if the prod/dest vectors are the same length
        if(fuel.batch[i].collapsed_iso.neutron_prod.size() != fuel.batch[i].collapsed_iso.neutron_dest.size()){
            cout << "Error. Neutron production/destruction rate vector length mismatch." << endl;
        }
        if(fuel.batch[i].collapsed_iso.neutron_prod[0] == 0 || fuel.batch[i].collapsed_iso.neutron_dest[0] == 0){
            fuel.batch[i].collapsed_iso.neutron_prod.erase(fuel.batch[i].collapsed_iso.neutron_prod.begin());
            fuel.batch[i].collapsed_iso.neutron_dest.erase(fuel.batch[i].collapsed_iso.neutron_dest.begin());
            fuel.batch[i].collapsed_iso.BUd.erase(fuel.batch[i].collapsed_iso.BUd.begin());
            fuel.batch[i].collapsed_iso.BU.erase(fuel.batch[i].collapsed_iso.BU.begin());
            fuel.batch[i].collapsed_iso.fluence.erase(fuel.batch[i].collapsed_iso.fluence.begin());
        }
    }
    return fuel;
};


map<int, double> tomass (int ti, double fluence, isoInformation isoinfo) {
    map<int, double> out = map<int, double>();
    double mass_i;
    int name_i;
    int nucid;
    for (int i = 0; i < isoinfo.iso_vector.size(); i++){
        name_i = isoinfo.iso_vector[i].name;
        nucid = pyne::nucname::zzaaam(name_i);
        mass_i = intpol(isoinfo.iso_vector[i].mass[ti-1],
                        isoinfo.iso_vector[i].mass[ti],
                        isoinfo.fluence[ti-1],
                        isoinfo.fluence[ti],
                        fluence);
        out[nucid] = mass_i;
    }
    return out;
}


fuelBundle phicalc_simple(fuelBundle core){
//updates the rflux of each batch in core.batch
//assumes the flux of a batch is proportional to the inverse neutron prod rate
    double maxphi = 0;
    //finds the inverse of neutron production at the batch_fluence
    //stores it in core.batch.rflux
    for(int i = 0; i < core.batch.size(); i++){
        int ii;
        if(core.batch[i].collapsed_iso.fluence.back() < core.batch[i].Fg){
            //cout << endl << "Maximum fluence error! Batch fluence exceeded max library fluence. (method2)" << endl;
            //cout << " Values on max fluence will be used. Do not trust results." << endl;
            for(int j = 0; j < core.batch.size(); j++){
                core.batch[j].rflux = 1;
            }
            return core;
        } else {
            for(ii = 0; core.batch[i].collapsed_iso.fluence[ii] < core.batch[i].Fg; ii++){}
        }
        if(ii == 0){
            core.batch[i].rflux = 1/core.batch[i].collapsed_iso.neutron_prod[0];
            if(maxphi < core.batch[i].rflux){
                maxphi = core.batch[i].rflux;
            }
        } else{
            core.batch[i].rflux =
                1/intpol(core.batch[i].collapsed_iso.neutron_prod[ii-1],core.batch[i].collapsed_iso.neutron_prod[ii], core.batch[i].collapsed_iso.fluence[ii-1], core.batch[i].collapsed_iso.fluence[ii], core.batch[i].Fg+1);
        }
        if(maxphi < core.batch[i].rflux){
            maxphi = core.batch[i].rflux;
        }
    }
    //normalizes all the flux values
    for(int i = 0; i < core.batch.size(); i++){
        core.batch[i].rflux /= maxphi;
        //cout << core.batch[i].rflux << endl;
    }
    //cout << endl;
    return core;
}

fuelBundle phicalc_eqpow(fuelBundle core){
    //updates the rflux of each batch in core.batch
    //assumes the flux of a batch is proportional to the inverse neutron prod rate
    double maxphi = 0;
    //finds the inverse of neutron production at the batch_fluence
    //stores it in core.batch.rflux
    for(int i = 0; i < core.batch.size(); i++){
        int ii;


        if(core.batch[i].collapsed_iso.fluence.back() < core.batch[i].Fg){
            //cout << endl << "Maximum fluence error! Batch fluence exceeded max library fluence. (method2)" << endl;
            //cout << "  Values on max fluence will be used. Do not trust results." << endl;
            for(int j = 0; j < core.batch.size(); j++){
                core.batch[j].rflux = 1;
            }
            return core;

        } else {
            for(ii = 0; core.batch[i].collapsed_iso.fluence[ii] < core.batch[i].Fg; ii++){}
        }

        if(ii == 0){
            //core.batch[i].rflux = 1/core.batch[i].collapsed_iso.neutron_prod[0];
            core.batch[i].rflux = (core.batch[i].collapsed_iso.fluence[ii+1]-core.batch[i].collapsed_iso.fluence[ii])
                /(core.batch[i].collapsed_iso.BU[ii+1]-core.batch[i].collapsed_iso.BU[ii]);

        } else{

            core.batch[i].rflux = (core.batch[i].collapsed_iso.fluence[ii]-core.batch[i].collapsed_iso.fluence[ii-1])
                /(core.batch[i].collapsed_iso.BU[ii]-core.batch[i].collapsed_iso.BU[ii-1]);
        }

        if(maxphi < core.batch[i].rflux){
            maxphi = core.batch[i].rflux;
        }
    }

    //normalizes all the flux values
    //cout << "EqPow Flux: ";
    for(int i = 0; i < core.batch.size(); i++){
        core.batch[i].rflux /= maxphi;
        //cout << core.batch[i].rflux << "   ";
    }
    //cout << endl;

    return core;
}

double nusigf_finder(batch_info batch){
    double Nusig_f;
    int ii;
    double prod = 0;
    int A = 235;

    //find the index to interpolate on
    for(ii = 0; batch.collapsed_iso.fluence[ii] < batch.Fg; ii++);
    if(batch.collapsed_iso.fluence.back() < batch.Fg){
        ii = batch.collapsed_iso.fluence.size() - 1;
    }

    prod = intpol(batch.collapsed_iso.neutron_prod[ii-1], batch.collapsed_iso.neutron_prod[ii], batch.collapsed_iso.fluence[ii-1], batch.collapsed_iso.fluence[ii], batch.Fg);


    if(ii == 0){
        prod = batch.collapsed_iso.neutron_prod[0];
    }

    //Nusig_f[cm-1] = prod rate / Avagadros number * mass * barn
    //Nusig_f = prod/(6.0221413E+23)*A*1E24/1000;
    Nusig_f = prod*0.01097;

    return Nusig_f;
}

double siga_finder(batch_info batch){
//returns the macroscopic cross section at fluence Fg
    double sig_a;
    int ii;
    double dest = 0;
    int A = 235;

    //find the index to interpolate on
    for(ii = 0; batch.collapsed_iso.fluence[ii] < batch.Fg; ii++);
    if(batch.collapsed_iso.fluence.back() < batch.Fg){
        ii = batch.collapsed_iso.fluence.size() - 1;
    }

    if(ii == 0){
        dest = batch.collapsed_iso.neutron_dest[0];
    }


    dest = intpol(batch.collapsed_iso.neutron_dest[ii-1], batch.collapsed_iso.neutron_dest[ii], batch.collapsed_iso.fluence[ii-1], batch.collapsed_iso.fluence[ii], batch.Fg);

    //sig_a[barn] = dest rate [cm2/kg] / Avagadros number * mass * barn
    //sig_a = dest/(6.0221413E+23)*A*1E24/1000;
    sig_a = dest*0.01097;

    return sig_a;
}


////////////////////////////////////////////////////////////////////////////////////////////
fuelBundle phicalc_cylindrical(fuelBundle core){
//cout << "cylindrical" << endl<<endl<<endl<<endl<<endl;

//the word 'region' is interchanged with the word 'batch' in comments, water layer is the outermost region

    int region = core.batch.size(); // number of regions (batches) not counting the outer region
    double delta = core.cylindrical_delta;
    double R[region+1]; //radial thickness of each region
    int N[region+1]; //number of mesh points in each region
    int NC[region+1]; //cumulative N
    int NTotal; //total number of mesh points
    double dd2[region+1]; // D/delta^2
    double Sigma_a[region+1]; //mac. abs. cs of each region
    double NuSigma_f[region+1]; //nu sigma f of each region
    double Sigma_tr[region+1]; //mac. transport cs of each region
    double D[region+1]; //diff coef. for each region
    double LSquared[region+1];
    double k = 1;
    double k_prev = 0.9;
    double prod, prod_prev;
    double flux[region+1], maxflux=0;
    double sum = 0;

    //set the radial thickness of each region
    R[0] = sqrt(core.fuel_area/region/3.141592);
    for(int i = 1; i < region; i++){
        R[i] = sqrt(core.fuel_area/region/3.141592*(i+1));
    }
    R[region] = R[region-1] + core.mod_thickness; //this is the moderator region


/*
    Sigma_a[0] = 0.0230;
    Sigma_a[1] = 0.0246;
    Sigma_a[2] = 0.0324;
    NuSigma_f[0] = 0.0184;
    NuSigma_f[1] = 0.0217;
    NuSigma_f[2] = 0.0382;*/
    //assign fuel cross sections
    for(int i = 0; i < region; i++){

        NuSigma_f[i] = nusigf_finder(core.batch[i]);

        Sigma_a[i] = siga_finder(core.batch[i]);

        Sigma_tr[i] = core.fuel_Sig_tr;
        D[i] = 1/(Sigma_tr[i]*3.);
        LSquared[i] = D[i]/Sigma_a[i];

    }

    //assign moderator cross sections
    Sigma_a[region] = 0.0066;// core.mod_Sig_a;
    Sigma_tr[region] = core.mod_Sig_tr;
    D[region] = 1/(Sigma_tr[region]*3.);
    LSquared[region] = D[region]/Sigma_a[region];
    NuSigma_f[region] = core.mod_Sig_f;

    //populate dd2
    for(int i = 0; i < region+1; i++){
        dd2[i] = D[i]/delta/delta;
    }

    /*
    cout << "Cross sections in cylindrical calc: " << endl << "Sig_a      Sig_f      R" << endl;
    for(int r = 0; r < region+1; r++){
        cout << Sigma_a[r] << "  " << NuSigma_f[r] << "  " << R[r] << endl;
    }*/

    //populate N, number of mesh points in each region
    N[0] = R[0]/delta;
    NC[0] = N[0];
    NTotal = N[0];

    for(int i = 1; i < region+1; i++){
        N[i] = ceil((R[i] - R[i-1]) / delta);
        NC[i] = NC[i-1] + N[i];
        NTotal += N[i];
    }
    NC[region] += 1;
    NTotal += 1;/*
    Eigen::MatrixXf A(NTotal, NTotal);
    Eigen::MatrixXf F(NTotal, 1);
    Eigen::MatrixXf phi(NTotal, 1);
    Eigen::MatrixXf phi_prev(NTotal, 1);
    Eigen::MatrixXf S(NTotal, 1);
    Eigen::MatrixXf S_prev(NTotal, 1);

    A.setZero();
    F.setZero();

    int jprev = 0;
    int j;

    int r = 0; //region index
    for(int i = 1; i < NTotal-1; i++){
        A(i, i-1) = (-1.)*dd2[r]*(2*i-1)/(2*i);
        A(i,i) = dd2[r]*2. + Sigma_a[r];
        A(i, i+1) = (-1.)*dd2[r]*(2*i+1)/(2*i);
        if(i == NC[r]){
            r += 1;
        }
        if(core.batch[r].rflux < 1.1 && core.batch[r].rflux > 0){
            phi_prev(i) = core.batch[r].rflux; //uses last runs results if available
        } else{
            phi_prev(i) = 1;
        }

    }
    A(0,0) = 1;
    A(0,1) = -1;
    A(NTotal-1,NTotal-1) = 1;
    phi_prev(NTotal-1) = 0;

    for(r = 0; r < region; r++){
        A(NC[r],NC[r]-1) = D[r];
        A(NC[r],NC[r]) = -D[r]-D[r+1];
        A(NC[r],NC[r]+1) = D[r+1];
    }

    r = 0;
    for(int i = 1; i < NTotal; i++){
        if(i != NC[r]){
            F(i) = NuSigma_f[r];
        }
        if( i == NC[r]){
            r += 1;
        }
        S_prev(i) = F(i)*phi_prev(i);
    }


    for(int iter = 0; iter < 10; iter++){
        phi = A.colPivHouseholderQr().solve(S_prev)/k_prev;


        if(!phi.allFinite()){
            phi = phi_prev;
        }

        S = F.array() * phi.array();

        prod = (0.25)*3.141592*NuSigma_f[0]*phi(0)*delta*delta;
        prod_prev = (0.25)*3.141592*NuSigma_f[0]*phi_prev(0)*delta*delta;
        r = 0;

        for(int i = 0; i < NTotal; i++){
            if(i == NC[r]){
                prod += 3.141592*NuSigma_f[r]*phi(i)*(i-0.25)*delta*delta;
                prod_prev += 3.141592*NuSigma_f[r]*phi_prev(i)*(i-0.25)*delta*delta;
                r += 1;
                prod += 3.141592*NuSigma_f[r]*phi(i)*(i-0.25)*delta*delta;
                prod_prev += 3.141592*NuSigma_f[r]*phi_prev(i)*(i-0.25)*delta*delta;
            } else {
                prod += 2.*3.141592*NuSigma_f[r]*phi(i)*i*delta*delta;
                prod_prev += 2.*3.141592*NuSigma_f[r]*phi_prev(i)*i*delta*delta;
            }
        }

        if(abs((k_prev-k)/k) < 0.01 && iter > 3){
            break;
        }
        //cout << "prod: " << prod << "  prod_prev: " << prod_prev << "  k: " << prod/prod_prev << endl;
        k = prod/prod_prev*k_prev;
        phi_prev = phi;
        k_prev = k;
        S_prev = S;


    }

    //NO NEED TO normalize phi
    //phi = phi.array()/phi.maxCoeff();

    //find area weighted average phi per batch
    r = 0;
    flux[0] = 0;
    for(int i = 0; i < NTotal; i++){
        flux[r] += phi(i)*(2*(i+1)-1);
        sum += (2*(i+1)-1);

        if(i == NC[r] || i == NTotal-1){
            flux[r] /= sum;
            sum = 0;
            r += 1;
            flux[r] = 0;
        }
    }
    for(r = 0; r < region+1; r++){
        if(flux[r] > maxflux){
            maxflux = flux[r];
        }
    }

    //cout << "flux: ";
    //cout << core.fuel_area << " " << delta << " | " << Sigma_a[0] << " " << Sigma_a[1] << " " << Sigma_a[2] << " | "
    //<< NuSigma_f[0] << " " << NuSigma_f[1] << " " << NuSigma_f[2] << " | ";
    //normalize the fluxes
    for(r = 0; r < region+1; r++){
        flux[r] /= maxflux;
        //cout << flux[r] << " ";
    }
    //ĺcout << endl;




    //NO NEED TO normalize phi
    //phi = phi.array()/phi.maxCoeff();

    //find area weighted average phi per batch
    r = 0;
    flux[0] = 0;
    for(int i = 0; i < NTotal; i++){
        flux[r] += phi(i)*(2*(i+1)-1);
        sum += (2*(i+1)-1);

        if(i == NC[r] || i == NTotal-1){
            flux[r] /= sum;
            sum = 0;
            r += 1;
            flux[r] = 0;
        }
    }
    for(r = 0; r < region+1; r++){
        if(flux[r] > maxflux){
            maxflux = flux[r];
        }
    }

    //cout << "flux: ";
    //cout << core.fuel_area << " " << delta << " | " << Sigma_a[0] << " " << Sigma_a[1] << " " << Sigma_a[2] << " | "
    //<< NuSigma_f[0] << " " << NuSigma_f[1] << " " << NuSigma_f[2] << " | ";
    //normalize the fluxes
    for(r = 0; r < region+1; r++){
        flux[r] /= maxflux;
        //cout << flux[r] << " ";
    }
    //cout << endl;


    for(int i = 0; i < core.batch.size(); i++){
        core.batch[i].rflux = flux[i];
    }
*/
    return core;
}
/////////////////////////////////////////////////////////////////////////////////////////////

double kcalc(fuelBundle core){
    //uses the fluence values in core.batch.Fg to calculate the core criticality
    //boost::timer t;
    int N = core.batch.size();
    double prod_tot = 0;
    double dest_tot = 0;
    double pnl = core.pnl;


    //finds the index j where fluence is just under target, then interpolates
    //  to find the prod and dest values for each batch
    int j;
    for(int i = 0; i < N; i++){

        if(core.batch[i].collapsed_iso.fluence.back() < core.batch[i].Fg){
            //cout << endl << "Maximum fluence error! Batch fluence exceeded max library fluence. (kcalc)" << endl;
            //cout << "  Values on max fluence will be used. Do not trust results." << endl;
            j = core.batch[i].collapsed_iso.fluence.size() - 1;
        } else {
            //find dicrete point j to interpolate on
            for(j = 0; core.batch[i].collapsed_iso.fluence[j] < core.batch[i].Fg; j++){}
        }

        //add the production rate of batch i to total production
        prod_tot += intpol(core.batch[i].collapsed_iso.neutron_prod[j-1], core.batch[i].collapsed_iso.neutron_prod[j], core.batch[i].collapsed_iso.fluence[j-1], core.batch[i].collapsed_iso.fluence[j], core.batch[i].Fg);

        //add structural material production of this batch, scaled up using disadvantage factor
        //prod_tot += core.struct_prod * core.batch[i].DA;

        //add the destruction rate of batch i to total destruction
        dest_tot += intpol(core.batch[i].collapsed_iso.neutron_dest[j-1], core.batch[i].collapsed_iso.neutron_dest[j], core.batch[i].collapsed_iso.fluence[j-1], core.batch[i].collapsed_iso.fluence[j], core.batch[i].Fg);

        //add structural material destruction of this batch, scaled up using disadvantage factor
        //dest_tot += core.struct_dest * core.batch[i].DA;
        //cout << "diff: " << prod_tot - dest_tot << endl;
    }
    if(pnl < 0 || pnl > 1.5){
        cout << endl << "Error in core nonleakage! Assumed new value: 1.000" << endl;
        pnl = 1.000;
        }

    //std::cout << "k calc time " << t.elapsed() << std::endl;
    return prod_tot * pnl / dest_tot;
}

double CR_finder(fuelBundle core){
    //
    //i is the batch number starting from zero
    //boost::timer t;
    double FP = 0, FP0 = 0, FP1 = 0;
    double fissile = 0, fissile0 = 0, fissile1 = 0;
    double ini_fissile = 0;
    double CR;
    int ii, ZZ;

    //cout << "\nCR_finder starting." << endl << "  Upper: " << core.CR_upper << "  lower: " << core.CR_lower << endl;


    for(int i = 0; i < core.batch.size(); i++){
        for(ii = 0; core.batch[i].collapsed_iso.fluence[ii] < core.batch[i].Fg; ii++){}
        if(ii == 0){
            ii = 1;
        }

        for(int j = 0; j < core.batch[i].collapsed_iso.iso_vector.size(); j++){
            //convert name to mass number
            ZZ = core.batch[i].collapsed_iso.iso_vector[j].name;
            ZZ = ZZ % 10000;
            ZZ /= 10;

            //add up the FP
            if(ZZ < core.CR_upper && ZZ > core.CR_lower){
                //cout << "    Batch" << i+1 << " ZZ: " << ZZ << endl;
                //interpolation will be done at the end
                FP0 += core.batch[i].collapsed_iso.iso_vector[j].mass[ii-1];
                FP1 += core.batch[i].collapsed_iso.iso_vector[j].mass[ii];
            }

            //add up fissiles
            for(int fis = 0; fis < core.CR_fissile.size(); fis++){
                if(core.batch[i].collapsed_iso.iso_vector[j].name == core.CR_fissile[fis]){
                    fissile0 += core.batch[i].collapsed_iso.iso_vector[j].mass[ii-1];
                    fissile1 += core.batch[i].collapsed_iso.iso_vector[j].mass[ii];

                    ini_fissile += core.batch[i].collapsed_iso.iso_vector[j].mass[0]; //mass at fluence zero
                }
            }


        }


        // recycling variable FP0 here to check greater than zero
        FP0 = intpol(FP0, FP1, core.batch[i].collapsed_iso.fluence[ii-1], core.batch[i].collapsed_iso.fluence[ii], core.batch[i].Fg);
        if(FP0 > 0){
            FP += FP0;
        }

        fissile += intpol(fissile0, fissile1, core.batch[i].collapsed_iso.fluence[ii-1], core.batch[i].collapsed_iso.fluence[ii], core.batch[i].Fg);

        FP0 = 0;
        FP1 = 0;
        fissile0 = 0;
        fissile1 = 0;
    }

    //cout << "  FP: " << FP << "  fiss: " << fissile << "  ini_fiss: " << ini_fissile << "      CR: " << (FP+fissile-ini_fissile)/FP << endl;

    if(FP > 0){
        CR = (FP+fissile-ini_fissile)/FP;
    } else {
        CR  = 0;
    }
    //std::cout << "CR Calc " << t.elapsed() << std::endl;
    return CR;
}


fuelBundle burnupcalc(fuelBundle core, int mode, int DA_mode, double delta) {
timestamp_t t0 = get_timestamp();
    //this function only uses the COLLAPSED_ISO of each BATCH in the structure CORE
    //all factors that contribute to a change in neutron prod/dest rates have to be factored
    //      before calling this function
    //cout << endl << "Burnupcalc" << endl;
    boost::timer t;

    int N = core.batch.size(); //number of batches
    double dt = delta*24*60*60; //days to [s]
    double kcore, kcore_prev;
    double y0, y1;
    double burnup = 0;

/*
    cout << "burnupcalc starting neutron prods: " << core.batch[0].collapsed_iso.neutron_prod[0] << " " << core.batch[1].collapsed_iso.neutron_prod[0] << " " << core.batch[2].collapsed_iso.neutron_prod[0]  << endl;
    cout << "   " << core.batch[0].collapsed_iso.BU[1] << " " << core.batch[1].collapsed_iso.BU[1] << " " << core.batch[2].collapsed_iso.BU[1]  << endl;
    cout << "   " << core.batch[0].collapsed_iso.fluence[1] << " " << core.batch[1].collapsed_iso.fluence[1] << " " << core.batch[2].collapsed_iso.fluence[1]  << endl;
    cout << "     " << dt << endl;
    cout << " " << core.batch[0].batch_fluence << " " << core.batch[1].batch_fluence << " " << core.batch[2].batch_fluence  << endl;
*/

    //assign the batch_fluence to Fg
    for(int i = 0; i < N; i++){
        core.batch[i].Fg = core.batch[i].batch_fluence;
    }

    kcore = 3.141592;

    //more forward in time until kcore drops under 1
    while(kcore > 1){
        kcore_prev = kcore;
        //if(kcore != 3.141592){cout << "kcore: " << kcore << endl;}

        //find the normalized relative flux of each batch
        if(mode == 1){
            //simplest case, all batches get the same flux
            for(int i = 0; i < N; i++){
                core.batch[i].rflux = 1;
            }
        }else if(mode == 2){
            //inverse-production flux calculation
            core = phicalc_simple(core);
        }else if(mode == 3){
            core = phicalc_cylindrical(core);
        }else if(mode == 0){
            core = phicalc_eqpow(core);
        }else{
            cout << endl << "Error in mode input for batch-level flux calculation." << endl;
            return core;
        }

        //disadvantage calculation
        if(DA_mode == 1){
            core = DA_calc(core);
        }

        //CR calculation
        //cout << CR_numerator(core, 5)/CR_denominator(core, 5) << " ";

        //update fluences
        for(int i = 0; i < N; i++){
            //cout << "flux: " << core.batch[i].rflux << endl;
            core.batch[i].Fg += core.batch[i].rflux * core.base_flux * dt;

        }
        //core.CR = CR_finder(core);
        kcore = kcalc(core);

        /*if(core.CR_target != 0 && std::abs(core.CR - core.CR_target)/core.CR < 0.1){
            break;
        }*/
    }
    core.CR = CR_finder(core);
    //cout << " CR: " << core.CR << endl;


    //update core fluences and CR
    for(int i = 0; i < N; i++){
        core.batch[i].batch_fluence = intpol(core.batch[i].Fg - (core.batch[i].rflux * core.base_flux * dt), core.batch[i].Fg, kcore_prev, kcore, 1);
        //cout << "batch " << i+1 << " CR: " << core.batch[i].CR << endl;
        //cout << "  fluence end of burnupcalc: " << core.batch[i].batch_fluence << endl;

    }
{'fraction_non_alphanumeric': 0.17712901223443664, 'fraction_numerical': 0.022506557738130994, 'mean_word_length': 2.1384402689624835, 'pattern_counts': {'":': 0, '<': 895, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2180113', 'n_tokens_mistral': 19154, 'n_tokens_neox': 17794, 'n_words': 3710}