/src/enzo/po_solver.C

Large files files are truncated, but you can click here to view the full file

/* THIS FILE HAS BEEN AUTO-GENERATED.  DO NOT EDIT. */

/* This is C++ code to read HDF5 files for
   reaction rates, cooling rates, and initial
   conditions for the chemical network defined
   by the user.  In addition, this contains
   code for calculating temperature from the
   gas energy and computing the RHS and the
   Jacobian of the system of equations which
   will be fed into the solver.
*/


#include "po_solver.h"

po_data *po_setup_data(void) {

    po_data *data = (po_data *) malloc(sizeof(po_data));

    /* Temperature-related pieces */
    data->bounds[0] = 1.0;
    data->bounds[1] = 1e+12;
    data->nbins = 1024 - 1;
    data->dbin = (log(data->bounds[1]) - log(data->bounds[0])) / data->nbins;
    data->idbin = 1.0L / data->dbin;

    /* Redshift-related pieces */
    data->z_bounds[0] = 0.0;
    data->z_bounds[1] = 15.0;
    data->n_zbins = 100 - 1;
    data->d_zbin = (log(data->z_bounds[1] + 1.0) - log(data->z_bounds[0] + 1.0)) / data->n_zbins;
    data->id_zbin = 1.0L / data->d_zbin;
    
    po_read_rate_tables(data);
    fprintf(stderr, "Successfully read in rate tables.\n");

    po_read_cooling_tables(data);
    fprintf(stderr, "Successfully read in cooling rate tables.\n");

    return data;

}


int po_main(int argc, char** argv)
{
    po_data *data = po_setup_data();

    /* Initial conditions */

    hid_t file_id = H5Fopen("po_initial_conditions.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
    if (file_id < 0) {fprintf(stderr, "Failed to open "
        "po_initial_conditions.h5 so dying.\n");
        return(1);}

    /* Allocate the correct number of cells */
    hsize_t dims; /* We have flat versus number of species */

    /* Check gas energy to get the number of cells */
    fprintf(stderr, "Getting dimensionality from ge:\n");
    herr_t status = H5LTget_dataset_info(file_id, "/ge", &dims, NULL, NULL);
    if(status == -1) {
        fprintf(stderr, "Error opening initial conditions file.\n");
        return 1;
    }
    fprintf(stderr, "  ncells = % 3i\n", (int) dims);
    data->ncells = dims;

    int N = 16;

    double *atol, *rtol;
    atol = (double *) alloca(N * dims * sizeof(double));
    rtol = (double *) alloca(N * dims * sizeof(double));

    double *tics = (double *) alloca(dims * sizeof(double));
    double *ics = (double *) alloca(dims * N * sizeof(double));
    double *input = (double *) alloca(dims * N * sizeof(double));
    
    unsigned int i = 0, j;
    
    fprintf(stderr, "Reading I.C. for /HI\n");
    H5LTread_dataset_double(file_id, "/HI", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 1.00794; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "HI[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /HII\n");
    H5LTread_dataset_double(file_id, "/HII", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 1.00794; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "HII[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /HeI\n");
    H5LTread_dataset_double(file_id, "/HeI", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 4.002602; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "HeI[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /HeII\n");
    H5LTread_dataset_double(file_id, "/HeII", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 4.002602; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "HeII[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /HeIII\n");
    H5LTread_dataset_double(file_id, "/HeIII", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 4.002602; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "HeIII[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /de\n");
    H5LTread_dataset_double(file_id, "/de", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 1.0; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "de[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /ge\n");
    H5LTread_dataset_double(file_id, "/ge", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 1.0; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "ge[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OI\n");
    H5LTread_dataset_double(file_id, "/OI", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OI[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OII\n");
    H5LTread_dataset_double(file_id, "/OII", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OII[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OIII\n");
    H5LTread_dataset_double(file_id, "/OIII", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OIII[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OIV\n");
    H5LTread_dataset_double(file_id, "/OIV", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OIV[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OV\n");
    H5LTread_dataset_double(file_id, "/OV", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OV[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OVI\n");
    H5LTread_dataset_double(file_id, "/OVI", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OVI[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OVII\n");
    H5LTread_dataset_double(file_id, "/OVII", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OVII[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OVIII\n");
    H5LTread_dataset_double(file_id, "/OVIII", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OVIII[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    
    fprintf(stderr, "Reading I.C. for /OIX\n");
    H5LTread_dataset_double(file_id, "/OIX", tics);
    for (j = 0; j < dims; j++) {
        ics[j * N + i] = tics[j] / 15.9994; /* Convert to number density */
        atol[j * N + i] = tics[j] * 1e-09;
        rtol[j * N + i] = 1e-09;
        if(j==0) {
            fprintf(stderr, "OIX[0] = %0.3g, atol => % 0.16g\n",
                    tics[j], atol[j]);
        }
    }
    i++;
    

    H5Fclose(file_id);

    double dtf = 3.1557e13;
    double dt = -1.0;
    double z = -1.0;
    for (i = 0; i < dims * N; i++) input[i] = ics[i];
    double ttot;
    ttot = dengo_evolve_po(dtf, dt, z, input, rtol, atol, dims, data);

    /* Write results to HDF5 file */
    file_id = H5Fcreate("po_solution.h5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    hsize_t dimsarr[1];
    dimsarr[0] = dims;
    i = 0;
    
    double HI[dims];
    for (j = 0; j < dims; j++) {
        HI[j] = input[j * N + i] * 1.00794;
    }
    fprintf(stderr, "Writing solution for /HI\n");
    H5LTmake_dataset_double(file_id, "/HI", 1, dimsarr, HI);
    i++;
    
    double HII[dims];
    for (j = 0; j < dims; j++) {
        HII[j] = input[j * N + i] * 1.00794;
    }
    fprintf(stderr, "Writing solution for /HII\n");
    H5LTmake_dataset_double(file_id, "/HII", 1, dimsarr, HII);
    i++;
    
    double HeI[dims];
    for (j = 0; j < dims; j++) {
        HeI[j] = input[j * N + i] * 4.002602;
    }
    fprintf(stderr, "Writing solution for /HeI\n");
    H5LTmake_dataset_double(file_id, "/HeI", 1, dimsarr, HeI);
    i++;
    
    double HeII[dims];
    for (j = 0; j < dims; j++) {
        HeII[j] = input[j * N + i] * 4.002602;
    }
    fprintf(stderr, "Writing solution for /HeII\n");
    H5LTmake_dataset_double(file_id, "/HeII", 1, dimsarr, HeII);
    i++;
    
    double HeIII[dims];
    for (j = 0; j < dims; j++) {
        HeIII[j] = input[j * N + i] * 4.002602;
    }
    fprintf(stderr, "Writing solution for /HeIII\n");
    H5LTmake_dataset_double(file_id, "/HeIII", 1, dimsarr, HeIII);
    i++;
    
    double de[dims];
    for (j = 0; j < dims; j++) {
        de[j] = input[j * N + i] * 1.0;
    }
    fprintf(stderr, "Writing solution for /de\n");
    H5LTmake_dataset_double(file_id, "/de", 1, dimsarr, de);
    i++;
    
    double ge[dims];
    for (j = 0; j < dims; j++) {
        ge[j] = input[j * N + i] * 1.0;
    }
    fprintf(stderr, "Writing solution for /ge\n");
    H5LTmake_dataset_double(file_id, "/ge", 1, dimsarr, ge);
    i++;
    
    double OI[dims];
    for (j = 0; j < dims; j++) {
        OI[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OI\n");
    H5LTmake_dataset_double(file_id, "/OI", 1, dimsarr, OI);
    i++;
    
    double OII[dims];
    for (j = 0; j < dims; j++) {
        OII[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OII\n");
    H5LTmake_dataset_double(file_id, "/OII", 1, dimsarr, OII);
    i++;
    
    double OIII[dims];
    for (j = 0; j < dims; j++) {
        OIII[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OIII\n");
    H5LTmake_dataset_double(file_id, "/OIII", 1, dimsarr, OIII);
    i++;
    
    double OIV[dims];
    for (j = 0; j < dims; j++) {
        OIV[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OIV\n");
    H5LTmake_dataset_double(file_id, "/OIV", 1, dimsarr, OIV);
    i++;
    
    double OV[dims];
    for (j = 0; j < dims; j++) {
        OV[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OV\n");
    H5LTmake_dataset_double(file_id, "/OV", 1, dimsarr, OV);
    i++;
    
    double OVI[dims];
    for (j = 0; j < dims; j++) {
        OVI[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OVI\n");
    H5LTmake_dataset_double(file_id, "/OVI", 1, dimsarr, OVI);
    i++;
    
    double OVII[dims];
    for (j = 0; j < dims; j++) {
        OVII[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OVII\n");
    H5LTmake_dataset_double(file_id, "/OVII", 1, dimsarr, OVII);
    i++;
    
    double OVIII[dims];
    for (j = 0; j < dims; j++) {
        OVIII[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OVIII\n");
    H5LTmake_dataset_double(file_id, "/OVIII", 1, dimsarr, OVIII);
    i++;
    
    double OIX[dims];
    for (j = 0; j < dims; j++) {
        OIX[j] = input[j * N + i] * 15.9994;
    }
    fprintf(stderr, "Writing solution for /OIX\n");
    H5LTmake_dataset_double(file_id, "/OIX", 1, dimsarr, OIX);
    i++;
    
    double temperature[dims];
    for (j = 0; j < dims; j++) {
    	temperature[j] = data->Ts[j];
    }
    H5LTmake_dataset_double(file_id, "/T", 1, dimsarr, temperature);
    double time[1];
    time[0] = ttot;
    double timestep[1];
    timestep[0] = dt;
    H5LTset_attribute_double(file_id, "/", "time", time, 1); 
    H5LTset_attribute_double(file_id, "/", "timestep", timestep, 1);
    H5Fclose(file_id);
    
    return 0;
}
 



double dengo_evolve_po (double dtf, double &dt, double z, double *input,
            double *rtol, double *atol, long long dims, po_data *data) {
    int i, j;
    hid_t file_id;
    /* fprintf(stderr, "  ncells = % 3i\n", (int) dims); */

    int N = 16;
    rhs_f f = calculate_rhs_po;
    jac_f jf = calculate_jacobian_po;
    if (dt < 0) dt = dtf; // / 1e3;
    data->current_z = z;
    int niter = 0;
    int siter = 0;
    double ttot = 0;
    double *scale = (double *) alloca(dims * N * sizeof(double));
    double *prev = (double *) alloca(dims * N * sizeof(double));
    for (i = 0; i < dims * N; i++) scale[i] = input[i];
    for (i = 0; i < dims * N; i++) prev[i] = input[i];
    double *u0 = (double *) alloca(N*dims*sizeof(double));
    double *s  = (double *) alloca(N*sizeof(double));
    double *gu = (double *) alloca(N*dims*sizeof(double));
    double *Ju = (double *) alloca(N*N*dims*sizeof(double));
    double floor_value = 1e-25;
    while (ttot < dtf) {
        int rv = BE_chem_solve(f, jf, input, dt, rtol, atol, dims, N,
                               scale, (void *) data, u0, s, gu, Ju);
        /*
        fprintf(stderr, "Return value [%d]: %i.  %0.5g / %0.5g = %0.5g (%0.5g)\n",
                niter, rv, ttot, dtf, ttot/dtf, dt);
        fprintf(stderr, "Value[0] = %0.5g %0.5g\n",
                input[0], prev[0]);
        */
        for (i = 0; i < dims * N; i++) {
            if (input[i] < 0) {
                rv = 1;
                break;
            }
        }
        if (rv == 0) {
	    if (siter == 500000) break;
	    siter++;
            if (siter % 100000 == 0) {
                fprintf(stderr, "Successful Iteration[%d]: (%0.4g) %0.16g / %0.16g\n",
                        siter, dt, ttot, dtf);
            }
            ttot += dt;
	    dt = DMIN(dt * 1.1, dtf - ttot);
	    
	    for (i = 0; i < dims * N; i++) prev[i] = input[i];
            for (i = 0; i < dims * N; i++) {     
                if (input[i] < floor_value) {
                  input[i] = floor_value;
                }
                scale[i] = input[i];
            }
        } else {
            dt /= 2.0;
            for (i = 0; i < dims * N; i++) {
                input[i] = prev[i];
                scale[i] = input[i];
            }
            if (dt/dtf < 1e-15)  {
                fprintf(stderr, "Dying! dt/dtf = %0.5g\n", dt/dtf);
                break;
            }
        }
        /* if ((dt/dtf < 1e-3) && ((dtf - ttot) != 0.0)) {
           fprintf(stderr, "Relative dt is small! (dt = %0.5g, dtf = %0.5g, dt/dtf = %0.5g)\n", dt, dtf, dt/dtf);
           } */
        niter++;
    }
    /* fprintf(stderr, "End: %0.5g / %0.5g (%0.5g)\n",
       ttot, dtf, dtf-ttot); */

    return ttot;
}
 


void po_read_rate_tables(po_data *data)
{
    hid_t file_id = H5Fopen("po_tables.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
    /* Allocate the correct number of rate tables */
    H5LTread_dataset_double(file_id, "/HeI_i", data->r_HeI_i);
    H5LTread_dataset_double(file_id, "/HeI_pi", data->r_HeI_pi);
    H5LTread_dataset_double(file_id, "/HeII_i", data->r_HeII_i);
    H5LTread_dataset_double(file_id, "/HeII_pi", data->r_HeII_pi);
    H5LTread_dataset_double(file_id, "/HeII_r", data->r_HeII_r);
    H5LTread_dataset_double(file_id, "/HeIII_r", data->r_HeIII_r);
    H5LTread_dataset_double(file_id, "/HI_i", data->r_HI_i);
    H5LTread_dataset_double(file_id, "/HI_pi", data->r_HI_pi);
    H5LTread_dataset_double(file_id, "/HII_r", data->r_HII_r);
    H5LTread_dataset_double(file_id, "/OI_i", data->r_OI_i);
    H5LTread_dataset_double(file_id, "/OI_pi", data->r_OI_pi);
    H5LTread_dataset_double(file_id, "/OII_i", data->r_OII_i);
    H5LTread_dataset_double(file_id, "/OII_pi", data->r_OII_pi);
    H5LTread_dataset_double(file_id, "/OII_r", data->r_OII_r);
    H5LTread_dataset_double(file_id, "/OIII_i", data->r_OIII_i);
    H5LTread_dataset_double(file_id, "/OIII_pi", data->r_OIII_pi);
    H5LTread_dataset_double(file_id, "/OIII_r", data->r_OIII_r);
    H5LTread_dataset_double(file_id, "/OIV_i", data->r_OIV_i);
    H5LTread_dataset_double(file_id, "/OIV_pi", data->r_OIV_pi);
    H5LTread_dataset_double(file_id, "/OIV_r", data->r_OIV_r);
    H5LTread_dataset_double(file_id, "/OIX_r", data->r_OIX_r);
    H5LTread_dataset_double(file_id, "/OV_i", data->r_OV_i);
    H5LTread_dataset_double(file_id, "/OV_pi", data->r_OV_pi);
    H5LTread_dataset_double(file_id, "/OV_r", data->r_OV_r);
    H5LTread_dataset_double(file_id, "/OVI_i", data->r_OVI_i);
    H5LTread_dataset_double(file_id, "/OVI_pi", data->r_OVI_pi);
    H5LTread_dataset_double(file_id, "/OVI_r", data->r_OVI_r);
    H5LTread_dataset_double(file_id, "/OVII_i", data->r_OVII_i);
    H5LTread_dataset_double(file_id, "/OVII_pi", data->r_OVII_pi);
    H5LTread_dataset_double(file_id, "/OVII_r", data->r_OVII_r);
    H5LTread_dataset_double(file_id, "/OVIII_i", data->r_OVIII_i);
    H5LTread_dataset_double(file_id, "/OVIII_pi", data->r_OVIII_pi);
    H5LTread_dataset_double(file_id, "/OVIII_r", data->r_OVIII_r);

    H5Fclose(file_id);
}

void po_read_cooling_tables(po_data *data)
{

    hid_t file_id = H5Fopen("po_tables.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
    /* Allocate the correct number of rate tables */
    H5LTread_dataset_double(file_id, "/compton_comp",
                            data->c_compton_comp);
    H5LTread_dataset_double(file_id, "/HeI_c_HeI_c",
                            data->c_HeI_c_HeI_c);
    H5LTread_dataset_double(file_id, "/HeI_ph_HeI_ph",
                            data->c_HeI_ph_HeI_ph);
    H5LTread_dataset_double(file_id, "/HeII_c_HeII_c",
                            data->c_HeII_c_HeII_c);
    H5LTread_dataset_double(file_id, "/HeII_ph_HeII_ph",
                            data->c_HeII_ph_HeII_ph);
    H5LTread_dataset_double(file_id, "/HeIII_c_HeIII_c",
                            data->c_HeIII_c_HeIII_c);
    H5LTread_dataset_double(file_id, "/HI_c_HI_c",
                            data->c_HI_c_HI_c);
    H5LTread_dataset_double(file_id, "/HI_ph_HI_ph",
                            data->c_HI_ph_HI_ph);
    H5LTread_dataset_double(file_id, "/HII_c_HII_c",
                            data->c_HII_c_HII_c);
    H5LTread_dataset_double(file_id, "/OI_c_OI_c",
                            data->c_OI_c_OI_c);
    H5LTread_dataset_double(file_id, "/OI_ph_OI_ph",
                            data->c_OI_ph_OI_ph);
    H5LTread_dataset_double(file_id, "/OII_c_OII_c",
                            data->c_OII_c_OII_c);
    H5LTread_dataset_double(file_id, "/OII_ph_OII_ph",
                            data->c_OII_ph_OII_ph);
    H5LTread_dataset_double(file_id, "/OIII_c_OIII_c",
                            data->c_OIII_c_OIII_c);
    H5LTread_dataset_double(file_id, "/OIII_ph_OIII_ph",
                            data->c_OIII_ph_OIII_ph);
    H5LTread_dataset_double(file_id, "/OIV_c_OIV_c",
                            data->c_OIV_c_OIV_c);
    H5LTread_dataset_double(file_id, "/OIV_ph_OIV_ph",
                            data->c_OIV_ph_OIV_ph);
    H5LTread_dataset_double(file_id, "/OIX_c_OIX_c",
                            data->c_OIX_c_OIX_c);
    H5LTread_dataset_double(file_id, "/OV_c_OV_c",
                            data->c_OV_c_OV_c);
    H5LTread_dataset_double(file_id, "/OV_ph_OV_ph",
                            data->c_OV_ph_OV_ph);
    H5LTread_dataset_double(file_id, "/OVI_c_OVI_c",
                            data->c_OVI_c_OVI_c);
    H5LTread_dataset_double(file_id, "/OVI_ph_OVI_ph",
                            data->c_OVI_ph_OVI_ph);
    H5LTread_dataset_double(file_id, "/OVII_c_OVII_c",
                            data->c_OVII_c_OVII_c);
    H5LTread_dataset_double(file_id, "/OVII_ph_OVII_ph",
                            data->c_OVII_ph_OVII_ph);
    H5LTread_dataset_double(file_id, "/OVIII_c_OVIII_c",
                            data->c_OVIII_c_OVIII_c);
    H5LTread_dataset_double(file_id, "/OVIII_ph_OVIII_ph",
                            data->c_OVIII_ph_OVIII_ph);

    H5Fclose(file_id);
}

 


void po_calculate_temperature(po_data *data,
                        double *input, int nstrip, int nchem)
{
    int i, j;
    double density;
    double kb = 1.3806504e-16; // Boltzmann constant [erg/K]
    double mh = 1.67e-24;
    double gamma = 5.e0/3.e0;
    
    /* Calculate total density */
    double HI;
    double HII;
    double HeI;
    double HeII;
    double HeIII;
    double de;
    double ge;
    double OI;
    double OII;
    double OIII;
    double OIV;
    double OV;
    double OVI;
    double OVII;
    double OVIII;
    double OIX;

    for (i = 0; i<nstrip; i++) {
        j = i * nchem;
        HI = input[j];
        /*fprintf(stderr, "HI[%d] = % 0.16g\n",
                i, HI);*/
        j++;
    
        HII = input[j];
        /*fprintf(stderr, "HII[%d] = % 0.16g\n",
                i, HII);*/
        j++;
    
        HeI = input[j];
        /*fprintf(stderr, "HeI[%d] = % 0.16g\n",
                i, HeI);*/
        j++;
    
        HeII = input[j];
        /*fprintf(stderr, "HeII[%d] = % 0.16g\n",
                i, HeII);*/
        j++;
    
        HeIII = input[j];
        /*fprintf(stderr, "HeIII[%d] = % 0.16g\n",
                i, HeIII);*/
        j++;
    
        de = input[j];
        /*fprintf(stderr, "de[%d] = % 0.16g\n",
                i, de);*/
        j++;
    
        ge = input[j];
        /*fprintf(stderr, "ge[%d] = % 0.16g\n",
                i, ge);*/
        j++;
    
        OI = input[j];
        /*fprintf(stderr, "OI[%d] = % 0.16g\n",
                i, OI);*/
        j++;
    
        OII = input[j];
        /*fprintf(stderr, "OII[%d] = % 0.16g\n",
                i, OII);*/
        j++;
    
        OIII = input[j];
        /*fprintf(stderr, "OIII[%d] = % 0.16g\n",
                i, OIII);*/
        j++;
    
        OIV = input[j];
        /*fprintf(stderr, "OIV[%d] = % 0.16g\n",
                i, OIV);*/
        j++;
    
        OV = input[j];
        /*fprintf(stderr, "OV[%d] = % 0.16g\n",
                i, OV);*/
        j++;
    
        OVI = input[j];
        /*fprintf(stderr, "OVI[%d] = % 0.16g\n",
                i, OVI);*/
        j++;
    
        OVII = input[j];
        /*fprintf(stderr, "OVII[%d] = % 0.16g\n",
                i, OVII);*/
        j++;
    
        OVIII = input[j];
        /*fprintf(stderr, "OVIII[%d] = % 0.16g\n",
                i, OVIII);*/
        j++;
    
        OIX = input[j];
        /*fprintf(stderr, "OIX[%d] = % 0.16g\n",
                i, OIX);*/
        j++;
    
        density = 1.00794*HI + 1.00794*HII + 4.002602*HeI + 4.002602*HeII + 4.002602*HeIII + 15.9994*OI + 15.9994*OII + 15.9994*OIII + 15.9994*OIV + 15.9994*OIX + 15.9994*OV + 15.9994*OVI + 15.9994*OVII + 15.9994*OVIII;
        data->Ts[i] = density*ge*mh/(kb*(HI/(gamma - 1.0) + HII/(gamma - 1.0) + HeI/(gamma - 1.0) + HeII/(gamma - 1.0) + HeIII/(gamma - 1.0) + OI/(gamma - 1.0) + OII/(gamma - 1.0) + OIII/(gamma - 1.0) + OIV/(gamma - 1.0) + OIX/(gamma - 1.0) + OV/(gamma - 1.0) + OVI/(gamma - 1.0) + OVII/(gamma - 1.0) + OVIII/(gamma - 1.0) + de/(gamma - 1.0)));
        if (data->Ts[i] < data->bounds[0]) {
            data->Ts[i] = data->bounds[0];
        } else if (data->Ts[i] > data->bounds[1]) {
            data->Ts[i] = data->bounds[1];
        }
        data->logTs[i] = log(data->Ts[i]);
        data->invTs[i] = 1.0 / data->Ts[i];
	data->dTs_ge[i] = 
        density*mh/(kb*(HI/(gamma - 1.0) + HII/(gamma - 1.0) + HeI/(gamma - 1.0) + HeII/(gamma - 1.0) + HeIII/(gamma - 1.0) + OI/(gamma - 1.0) + OII/(gamma - 1.0) + OIII/(gamma - 1.0) + OIV/(gamma - 1.0) + OIX/(gamma - 1.0) + OV/(gamma - 1.0) + OVI/(gamma - 1.0) + OVII/(gamma - 1.0) + OVIII/(gamma - 1.0) + de/(gamma - 1.0)));
        /*fprintf(stderr, "T[%d] = % 0.16g, density = % 0.16g\n",
                i, data->Ts[i], density);*/
    }
         
}
 


/*
   This setup may be different than the user may anticipate, as a result
   of the lockstep timestep we use for a pencil beam through the grid.
   As such, it accepts the number of things to interpolate and makes
   assumptions about the sizes of the rates.
*/

/* This also requires no templating other than for the solver name...*/
void po_interpolate_rates(po_data *data,
                    int nstrip)
{
    int i, bin_id, zbin_id;
    double lb, t1, t2;
    double lbz, z1, z2;
    int no_photo = 0;
    lb = log(data->bounds[0]);
    lbz = log(data->z_bounds[0] + 1.0);
    /*fprintf(stderr, "lb = % 0.16g, ub = % 0.16g\n", lb, ub);*/
    for (i = 0; i < nstrip; i++) {
        data->bin_id[i] = bin_id = (int) (data->idbin * (data->logTs[i] - lb));
        if (data->bin_id[i] <= 0) {
            data->bin_id[i] = 0;
        } else if (data->bin_id[i] >= data->nbins) {
            data->bin_id[i] = data->nbins - 1;
        }
        t1 = (lb + (bin_id    ) * data->dbin);
        t2 = (lb + (bin_id + 1) * data->dbin);
        data->Tdef[i] = (data->logTs[i] - t1)/(t2 - t1);
        data->dT[i] = (t2 - t1);
        /*fprintf(stderr, "INTERP: %d, bin_id = %d, dT = % 0.16g, T = % 0.16g, logT = % 0.16g\n",
                i, data->bin_id[i], data->dT[i], data->Ts[i],
                data->logTs[i]);*/
    }
    
    if ((data->current_z >= data->z_bounds[0]) && (data->current_z < data->z_bounds[1])) {
        zbin_id = (int) (data->id_zbin * (log(data->current_z + 1.0) - lbz));
        if (zbin_id <= 0) {
            zbin_id = 0;
        } else if (zbin_id >= data->n_zbins) {
            zbin_id = data->n_zbins - 1;
        }
        z1 = (lbz + (zbin_id    ) * data->d_zbin);
        z2 = (lbz + (zbin_id + 1) * data->d_zbin);
        data->zdef = (log(data->current_z + 1.0) - z1)/(z2 - z1);
        data->dz = (exp(z2) - exp(z1)); //note: given this, we don't have to divide rate of change by z
    } else {
        no_photo = 1;
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_HeI_i[i] = data->r_HeI_i[bin_id] +
            data->Tdef[i] * (data->r_HeI_i[bin_id+1] - data->r_HeI_i[bin_id]);
        data->drs_HeI_i[i] = (data->r_HeI_i[bin_id+1] - data->r_HeI_i[bin_id]);
        data->drs_HeI_i[i] /= data->dT[i];
	data->drs_HeI_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_HeI_pi[i] = 0.0;
            data->drs_HeI_pi[i] = 0.0;
        } else {
            data->rs_HeI_pi[i] = data->r_HeI_pi[zbin_id] +
                data->zdef * (data->r_HeI_pi[zbin_id+1] - data->r_HeI_pi[zbin_id]);
            data->drs_HeI_pi[i] = (data->r_HeI_pi[zbin_id+1] - data->r_HeI_pi[zbin_id]);
            data->drs_HeI_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_HeII_i[i] = data->r_HeII_i[bin_id] +
            data->Tdef[i] * (data->r_HeII_i[bin_id+1] - data->r_HeII_i[bin_id]);
        data->drs_HeII_i[i] = (data->r_HeII_i[bin_id+1] - data->r_HeII_i[bin_id]);
        data->drs_HeII_i[i] /= data->dT[i];
	data->drs_HeII_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_HeII_pi[i] = 0.0;
            data->drs_HeII_pi[i] = 0.0;
        } else {
            data->rs_HeII_pi[i] = data->r_HeII_pi[zbin_id] +
                data->zdef * (data->r_HeII_pi[zbin_id+1] - data->r_HeII_pi[zbin_id]);
            data->drs_HeII_pi[i] = (data->r_HeII_pi[zbin_id+1] - data->r_HeII_pi[zbin_id]);
            data->drs_HeII_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_HeII_r[i] = data->r_HeII_r[bin_id] +
            data->Tdef[i] * (data->r_HeII_r[bin_id+1] - data->r_HeII_r[bin_id]);
        data->drs_HeII_r[i] = (data->r_HeII_r[bin_id+1] - data->r_HeII_r[bin_id]);
        data->drs_HeII_r[i] /= data->dT[i];
	data->drs_HeII_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_HeIII_r[i] = data->r_HeIII_r[bin_id] +
            data->Tdef[i] * (data->r_HeIII_r[bin_id+1] - data->r_HeIII_r[bin_id]);
        data->drs_HeIII_r[i] = (data->r_HeIII_r[bin_id+1] - data->r_HeIII_r[bin_id]);
        data->drs_HeIII_r[i] /= data->dT[i];
	data->drs_HeIII_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_HI_i[i] = data->r_HI_i[bin_id] +
            data->Tdef[i] * (data->r_HI_i[bin_id+1] - data->r_HI_i[bin_id]);
        data->drs_HI_i[i] = (data->r_HI_i[bin_id+1] - data->r_HI_i[bin_id]);
        data->drs_HI_i[i] /= data->dT[i];
	data->drs_HI_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_HI_pi[i] = 0.0;
            data->drs_HI_pi[i] = 0.0;
        } else {
            data->rs_HI_pi[i] = data->r_HI_pi[zbin_id] +
                data->zdef * (data->r_HI_pi[zbin_id+1] - data->r_HI_pi[zbin_id]);
            data->drs_HI_pi[i] = (data->r_HI_pi[zbin_id+1] - data->r_HI_pi[zbin_id]);
            data->drs_HI_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_HII_r[i] = data->r_HII_r[bin_id] +
            data->Tdef[i] * (data->r_HII_r[bin_id+1] - data->r_HII_r[bin_id]);
        data->drs_HII_r[i] = (data->r_HII_r[bin_id+1] - data->r_HII_r[bin_id]);
        data->drs_HII_r[i] /= data->dT[i];
	data->drs_HII_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OI_i[i] = data->r_OI_i[bin_id] +
            data->Tdef[i] * (data->r_OI_i[bin_id+1] - data->r_OI_i[bin_id]);
        data->drs_OI_i[i] = (data->r_OI_i[bin_id+1] - data->r_OI_i[bin_id]);
        data->drs_OI_i[i] /= data->dT[i];
	data->drs_OI_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_OI_pi[i] = 0.0;
            data->drs_OI_pi[i] = 0.0;
        } else {
            data->rs_OI_pi[i] = data->r_OI_pi[zbin_id] +
                data->zdef * (data->r_OI_pi[zbin_id+1] - data->r_OI_pi[zbin_id]);
            data->drs_OI_pi[i] = (data->r_OI_pi[zbin_id+1] - data->r_OI_pi[zbin_id]);
            data->drs_OI_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OII_i[i] = data->r_OII_i[bin_id] +
            data->Tdef[i] * (data->r_OII_i[bin_id+1] - data->r_OII_i[bin_id]);
        data->drs_OII_i[i] = (data->r_OII_i[bin_id+1] - data->r_OII_i[bin_id]);
        data->drs_OII_i[i] /= data->dT[i];
	data->drs_OII_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_OII_pi[i] = 0.0;
            data->drs_OII_pi[i] = 0.0;
        } else {
            data->rs_OII_pi[i] = data->r_OII_pi[zbin_id] +
                data->zdef * (data->r_OII_pi[zbin_id+1] - data->r_OII_pi[zbin_id]);
            data->drs_OII_pi[i] = (data->r_OII_pi[zbin_id+1] - data->r_OII_pi[zbin_id]);
            data->drs_OII_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OII_r[i] = data->r_OII_r[bin_id] +
            data->Tdef[i] * (data->r_OII_r[bin_id+1] - data->r_OII_r[bin_id]);
        data->drs_OII_r[i] = (data->r_OII_r[bin_id+1] - data->r_OII_r[bin_id]);
        data->drs_OII_r[i] /= data->dT[i];
	data->drs_OII_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OIII_i[i] = data->r_OIII_i[bin_id] +
            data->Tdef[i] * (data->r_OIII_i[bin_id+1] - data->r_OIII_i[bin_id]);
        data->drs_OIII_i[i] = (data->r_OIII_i[bin_id+1] - data->r_OIII_i[bin_id]);
        data->drs_OIII_i[i] /= data->dT[i];
	data->drs_OIII_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_OIII_pi[i] = 0.0;
            data->drs_OIII_pi[i] = 0.0;
        } else {
            data->rs_OIII_pi[i] = data->r_OIII_pi[zbin_id] +
                data->zdef * (data->r_OIII_pi[zbin_id+1] - data->r_OIII_pi[zbin_id]);
            data->drs_OIII_pi[i] = (data->r_OIII_pi[zbin_id+1] - data->r_OIII_pi[zbin_id]);
            data->drs_OIII_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OIII_r[i] = data->r_OIII_r[bin_id] +
            data->Tdef[i] * (data->r_OIII_r[bin_id+1] - data->r_OIII_r[bin_id]);
        data->drs_OIII_r[i] = (data->r_OIII_r[bin_id+1] - data->r_OIII_r[bin_id]);
        data->drs_OIII_r[i] /= data->dT[i];
	data->drs_OIII_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OIV_i[i] = data->r_OIV_i[bin_id] +
            data->Tdef[i] * (data->r_OIV_i[bin_id+1] - data->r_OIV_i[bin_id]);
        data->drs_OIV_i[i] = (data->r_OIV_i[bin_id+1] - data->r_OIV_i[bin_id]);
        data->drs_OIV_i[i] /= data->dT[i];
	data->drs_OIV_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_OIV_pi[i] = 0.0;
            data->drs_OIV_pi[i] = 0.0;
        } else {
            data->rs_OIV_pi[i] = data->r_OIV_pi[zbin_id] +
                data->zdef * (data->r_OIV_pi[zbin_id+1] - data->r_OIV_pi[zbin_id]);
            data->drs_OIV_pi[i] = (data->r_OIV_pi[zbin_id+1] - data->r_OIV_pi[zbin_id]);
            data->drs_OIV_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OIV_r[i] = data->r_OIV_r[bin_id] +
            data->Tdef[i] * (data->r_OIV_r[bin_id+1] - data->r_OIV_r[bin_id]);
        data->drs_OIV_r[i] = (data->r_OIV_r[bin_id+1] - data->r_OIV_r[bin_id]);
        data->drs_OIV_r[i] /= data->dT[i];
	data->drs_OIV_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OIX_r[i] = data->r_OIX_r[bin_id] +
            data->Tdef[i] * (data->r_OIX_r[bin_id+1] - data->r_OIX_r[bin_id]);
        data->drs_OIX_r[i] = (data->r_OIX_r[bin_id+1] - data->r_OIX_r[bin_id]);
        data->drs_OIX_r[i] /= data->dT[i];
	data->drs_OIX_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OV_i[i] = data->r_OV_i[bin_id] +
            data->Tdef[i] * (data->r_OV_i[bin_id+1] - data->r_OV_i[bin_id]);
        data->drs_OV_i[i] = (data->r_OV_i[bin_id+1] - data->r_OV_i[bin_id]);
        data->drs_OV_i[i] /= data->dT[i];
	data->drs_OV_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_OV_pi[i] = 0.0;
            data->drs_OV_pi[i] = 0.0;
        } else {
            data->rs_OV_pi[i] = data->r_OV_pi[zbin_id] +
                data->zdef * (data->r_OV_pi[zbin_id+1] - data->r_OV_pi[zbin_id]);
            data->drs_OV_pi[i] = (data->r_OV_pi[zbin_id+1] - data->r_OV_pi[zbin_id]);
            data->drs_OV_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OV_r[i] = data->r_OV_r[bin_id] +
            data->Tdef[i] * (data->r_OV_r[bin_id+1] - data->r_OV_r[bin_id]);
        data->drs_OV_r[i] = (data->r_OV_r[bin_id+1] - data->r_OV_r[bin_id]);
        data->drs_OV_r[i] /= data->dT[i];
	data->drs_OV_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OVI_i[i] = data->r_OVI_i[bin_id] +
            data->Tdef[i] * (data->r_OVI_i[bin_id+1] - data->r_OVI_i[bin_id]);
        data->drs_OVI_i[i] = (data->r_OVI_i[bin_id+1] - data->r_OVI_i[bin_id]);
        data->drs_OVI_i[i] /= data->dT[i];
	data->drs_OVI_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_OVI_pi[i] = 0.0;
            data->drs_OVI_pi[i] = 0.0;
        } else {
            data->rs_OVI_pi[i] = data->r_OVI_pi[zbin_id] +
                data->zdef * (data->r_OVI_pi[zbin_id+1] - data->r_OVI_pi[zbin_id]);
            data->drs_OVI_pi[i] = (data->r_OVI_pi[zbin_id+1] - data->r_OVI_pi[zbin_id]);
            data->drs_OVI_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OVI_r[i] = data->r_OVI_r[bin_id] +
            data->Tdef[i] * (data->r_OVI_r[bin_id+1] - data->r_OVI_r[bin_id]);
        data->drs_OVI_r[i] = (data->r_OVI_r[bin_id+1] - data->r_OVI_r[bin_id]);
        data->drs_OVI_r[i] /= data->dT[i];
	data->drs_OVI_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OVII_i[i] = data->r_OVII_i[bin_id] +
            data->Tdef[i] * (data->r_OVII_i[bin_id+1] - data->r_OVII_i[bin_id]);
        data->drs_OVII_i[i] = (data->r_OVII_i[bin_id+1] - data->r_OVII_i[bin_id]);
        data->drs_OVII_i[i] /= data->dT[i];
	data->drs_OVII_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_OVII_pi[i] = 0.0;
            data->drs_OVII_pi[i] = 0.0;
        } else {
            data->rs_OVII_pi[i] = data->r_OVII_pi[zbin_id] +
                data->zdef * (data->r_OVII_pi[zbin_id+1] - data->r_OVII_pi[zbin_id]);
            data->drs_OVII_pi[i] = (data->r_OVII_pi[zbin_id+1] - data->r_OVII_pi[zbin_id]);
            data->drs_OVII_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OVII_r[i] = data->r_OVII_r[bin_id] +
            data->Tdef[i] * (data->r_OVII_r[bin_id+1] - data->r_OVII_r[bin_id]);
        data->drs_OVII_r[i] = (data->r_OVII_r[bin_id+1] - data->r_OVII_r[bin_id]);
        data->drs_OVII_r[i] /= data->dT[i];
	data->drs_OVII_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OVIII_i[i] = data->r_OVIII_i[bin_id] +
            data->Tdef[i] * (data->r_OVIII_i[bin_id+1] - data->r_OVIII_i[bin_id]);
        data->drs_OVIII_i[i] = (data->r_OVIII_i[bin_id+1] - data->r_OVIII_i[bin_id]);
        data->drs_OVIII_i[i] /= data->dT[i];
	data->drs_OVIII_i[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->rs_OVIII_pi[i] = 0.0;
            data->drs_OVIII_pi[i] = 0.0;
        } else {
            data->rs_OVIII_pi[i] = data->r_OVIII_pi[zbin_id] +
                data->zdef * (data->r_OVIII_pi[zbin_id+1] - data->r_OVIII_pi[zbin_id]);
            data->drs_OVIII_pi[i] = (data->r_OVIII_pi[zbin_id+1] - data->r_OVIII_pi[zbin_id]);
            data->drs_OVIII_pi[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->rs_OVIII_r[i] = data->r_OVIII_r[bin_id] +
            data->Tdef[i] * (data->r_OVIII_r[bin_id+1] - data->r_OVIII_r[bin_id]);
        data->drs_OVIII_r[i] = (data->r_OVIII_r[bin_id+1] - data->r_OVIII_r[bin_id]);
        data->drs_OVIII_r[i] /= data->dT[i];
	data->drs_OVIII_r[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_compton_comp[i] = data->c_compton_comp[bin_id] +
            data->Tdef[i] * (data->c_compton_comp[bin_id+1] - data->c_compton_comp[bin_id]);
        data->dcs_compton_comp[i] = (data->c_compton_comp[bin_id+1] - data->c_compton_comp[bin_id]);;
        data->dcs_compton_comp[i] /= data->dT[i];
	data->dcs_compton_comp[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_HeI_c_HeI_c[i] = data->c_HeI_c_HeI_c[bin_id] +
            data->Tdef[i] * (data->c_HeI_c_HeI_c[bin_id+1] - data->c_HeI_c_HeI_c[bin_id]);
        data->dcs_HeI_c_HeI_c[i] = (data->c_HeI_c_HeI_c[bin_id+1] - data->c_HeI_c_HeI_c[bin_id]);;
        data->dcs_HeI_c_HeI_c[i] /= data->dT[i];
	data->dcs_HeI_c_HeI_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_HeI_ph_HeI_ph[i] = 0.0;
            data->dcs_HeI_ph_HeI_ph[i] = 0.0;
        } else {
            data->cs_HeI_ph_HeI_ph[i] = data->c_HeI_ph_HeI_ph[zbin_id] +
                data->zdef * (data->c_HeI_ph_HeI_ph[zbin_id+1] - data->c_HeI_ph_HeI_ph[zbin_id]);
            data->dcs_HeI_ph_HeI_ph[i] = (data->c_HeI_ph_HeI_ph[zbin_id+1] - data->c_HeI_ph_HeI_ph[zbin_id]);;
            data->dcs_HeI_ph_HeI_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_HeII_c_HeII_c[i] = data->c_HeII_c_HeII_c[bin_id] +
            data->Tdef[i] * (data->c_HeII_c_HeII_c[bin_id+1] - data->c_HeII_c_HeII_c[bin_id]);
        data->dcs_HeII_c_HeII_c[i] = (data->c_HeII_c_HeII_c[bin_id+1] - data->c_HeII_c_HeII_c[bin_id]);;
        data->dcs_HeII_c_HeII_c[i] /= data->dT[i];
	data->dcs_HeII_c_HeII_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_HeII_ph_HeII_ph[i] = 0.0;
            data->dcs_HeII_ph_HeII_ph[i] = 0.0;
        } else {
            data->cs_HeII_ph_HeII_ph[i] = data->c_HeII_ph_HeII_ph[zbin_id] +
                data->zdef * (data->c_HeII_ph_HeII_ph[zbin_id+1] - data->c_HeII_ph_HeII_ph[zbin_id]);
            data->dcs_HeII_ph_HeII_ph[i] = (data->c_HeII_ph_HeII_ph[zbin_id+1] - data->c_HeII_ph_HeII_ph[zbin_id]);;
            data->dcs_HeII_ph_HeII_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_HeIII_c_HeIII_c[i] = data->c_HeIII_c_HeIII_c[bin_id] +
            data->Tdef[i] * (data->c_HeIII_c_HeIII_c[bin_id+1] - data->c_HeIII_c_HeIII_c[bin_id]);
        data->dcs_HeIII_c_HeIII_c[i] = (data->c_HeIII_c_HeIII_c[bin_id+1] - data->c_HeIII_c_HeIII_c[bin_id]);;
        data->dcs_HeIII_c_HeIII_c[i] /= data->dT[i];
	data->dcs_HeIII_c_HeIII_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_HI_c_HI_c[i] = data->c_HI_c_HI_c[bin_id] +
            data->Tdef[i] * (data->c_HI_c_HI_c[bin_id+1] - data->c_HI_c_HI_c[bin_id]);
        data->dcs_HI_c_HI_c[i] = (data->c_HI_c_HI_c[bin_id+1] - data->c_HI_c_HI_c[bin_id]);;
        data->dcs_HI_c_HI_c[i] /= data->dT[i];
	data->dcs_HI_c_HI_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_HI_ph_HI_ph[i] = 0.0;
            data->dcs_HI_ph_HI_ph[i] = 0.0;
        } else {
            data->cs_HI_ph_HI_ph[i] = data->c_HI_ph_HI_ph[zbin_id] +
                data->zdef * (data->c_HI_ph_HI_ph[zbin_id+1] - data->c_HI_ph_HI_ph[zbin_id]);
            data->dcs_HI_ph_HI_ph[i] = (data->c_HI_ph_HI_ph[zbin_id+1] - data->c_HI_ph_HI_ph[zbin_id]);;
            data->dcs_HI_ph_HI_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_HII_c_HII_c[i] = data->c_HII_c_HII_c[bin_id] +
            data->Tdef[i] * (data->c_HII_c_HII_c[bin_id+1] - data->c_HII_c_HII_c[bin_id]);
        data->dcs_HII_c_HII_c[i] = (data->c_HII_c_HII_c[bin_id+1] - data->c_HII_c_HII_c[bin_id]);;
        data->dcs_HII_c_HII_c[i] /= data->dT[i];
	data->dcs_HII_c_HII_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OI_c_OI_c[i] = data->c_OI_c_OI_c[bin_id] +
            data->Tdef[i] * (data->c_OI_c_OI_c[bin_id+1] - data->c_OI_c_OI_c[bin_id]);
        data->dcs_OI_c_OI_c[i] = (data->c_OI_c_OI_c[bin_id+1] - data->c_OI_c_OI_c[bin_id]);;
        data->dcs_OI_c_OI_c[i] /= data->dT[i];
	data->dcs_OI_c_OI_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_OI_ph_OI_ph[i] = 0.0;
            data->dcs_OI_ph_OI_ph[i] = 0.0;
        } else {
            data->cs_OI_ph_OI_ph[i] = data->c_OI_ph_OI_ph[zbin_id] +
                data->zdef * (data->c_OI_ph_OI_ph[zbin_id+1] - data->c_OI_ph_OI_ph[zbin_id]);
            data->dcs_OI_ph_OI_ph[i] = (data->c_OI_ph_OI_ph[zbin_id+1] - data->c_OI_ph_OI_ph[zbin_id]);;
            data->dcs_OI_ph_OI_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OII_c_OII_c[i] = data->c_OII_c_OII_c[bin_id] +
            data->Tdef[i] * (data->c_OII_c_OII_c[bin_id+1] - data->c_OII_c_OII_c[bin_id]);
        data->dcs_OII_c_OII_c[i] = (data->c_OII_c_OII_c[bin_id+1] - data->c_OII_c_OII_c[bin_id]);;
        data->dcs_OII_c_OII_c[i] /= data->dT[i];
	data->dcs_OII_c_OII_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_OII_ph_OII_ph[i] = 0.0;
            data->dcs_OII_ph_OII_ph[i] = 0.0;
        } else {
            data->cs_OII_ph_OII_ph[i] = data->c_OII_ph_OII_ph[zbin_id] +
                data->zdef * (data->c_OII_ph_OII_ph[zbin_id+1] - data->c_OII_ph_OII_ph[zbin_id]);
            data->dcs_OII_ph_OII_ph[i] = (data->c_OII_ph_OII_ph[zbin_id+1] - data->c_OII_ph_OII_ph[zbin_id]);;
            data->dcs_OII_ph_OII_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OIII_c_OIII_c[i] = data->c_OIII_c_OIII_c[bin_id] +
            data->Tdef[i] * (data->c_OIII_c_OIII_c[bin_id+1] - data->c_OIII_c_OIII_c[bin_id]);
        data->dcs_OIII_c_OIII_c[i] = (data->c_OIII_c_OIII_c[bin_id+1] - data->c_OIII_c_OIII_c[bin_id]);;
        data->dcs_OIII_c_OIII_c[i] /= data->dT[i];
	data->dcs_OIII_c_OIII_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_OIII_ph_OIII_ph[i] = 0.0;
            data->dcs_OIII_ph_OIII_ph[i] = 0.0;
        } else {
            data->cs_OIII_ph_OIII_ph[i] = data->c_OIII_ph_OIII_ph[zbin_id] +
                data->zdef * (data->c_OIII_ph_OIII_ph[zbin_id+1] - data->c_OIII_ph_OIII_ph[zbin_id]);
            data->dcs_OIII_ph_OIII_ph[i] = (data->c_OIII_ph_OIII_ph[zbin_id+1] - data->c_OIII_ph_OIII_ph[zbin_id]);;
            data->dcs_OIII_ph_OIII_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OIV_c_OIV_c[i] = data->c_OIV_c_OIV_c[bin_id] +
            data->Tdef[i] * (data->c_OIV_c_OIV_c[bin_id+1] - data->c_OIV_c_OIV_c[bin_id]);
        data->dcs_OIV_c_OIV_c[i] = (data->c_OIV_c_OIV_c[bin_id+1] - data->c_OIV_c_OIV_c[bin_id]);;
        data->dcs_OIV_c_OIV_c[i] /= data->dT[i];
	data->dcs_OIV_c_OIV_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_OIV_ph_OIV_ph[i] = 0.0;
            data->dcs_OIV_ph_OIV_ph[i] = 0.0;
        } else {
            data->cs_OIV_ph_OIV_ph[i] = data->c_OIV_ph_OIV_ph[zbin_id] +
                data->zdef * (data->c_OIV_ph_OIV_ph[zbin_id+1] - data->c_OIV_ph_OIV_ph[zbin_id]);
            data->dcs_OIV_ph_OIV_ph[i] = (data->c_OIV_ph_OIV_ph[zbin_id+1] - data->c_OIV_ph_OIV_ph[zbin_id]);;
            data->dcs_OIV_ph_OIV_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OIX_c_OIX_c[i] = data->c_OIX_c_OIX_c[bin_id] +
            data->Tdef[i] * (data->c_OIX_c_OIX_c[bin_id+1] - data->c_OIX_c_OIX_c[bin_id]);
        data->dcs_OIX_c_OIX_c[i] = (data->c_OIX_c_OIX_c[bin_id+1] - data->c_OIX_c_OIX_c[bin_id]);;
        data->dcs_OIX_c_OIX_c[i] /= data->dT[i];
	data->dcs_OIX_c_OIX_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OV_c_OV_c[i] = data->c_OV_c_OV_c[bin_id] +
            data->Tdef[i] * (data->c_OV_c_OV_c[bin_id+1] - data->c_OV_c_OV_c[bin_id]);
        data->dcs_OV_c_OV_c[i] = (data->c_OV_c_OV_c[bin_id+1] - data->c_OV_c_OV_c[bin_id]);;
        data->dcs_OV_c_OV_c[i] /= data->dT[i];
	data->dcs_OV_c_OV_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
     …