/src/enzo/po_solver.C

/* 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++) {
        if (no_photo) {
            data->cs_OV_ph_OV_ph[i] = 0.0;
            data->dcs_OV_ph_OV_ph[i] = 0.0;
        } else {
            data->cs_OV_ph_OV_ph[i] = data->c_OV_ph_OV_ph[zbin_id] +
                data->zdef * (data->c_OV_ph_OV_ph[zbin_id+1] - data->c_OV_ph_OV_ph[zbin_id]);
            data->dcs_OV_ph_OV_ph[i] = (data->c_OV_ph_OV_ph[zbin_id+1] - data->c_OV_ph_OV_ph[zbin_id]);;
            data->dcs_OV_ph_OV_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OVI_c_OVI_c[i] = data->c_OVI_c_OVI_c[bin_id] +
            data->Tdef[i] * (data->c_OVI_c_OVI_c[bin_id+1] - data->c_OVI_c_OVI_c[bin_id]);
        data->dcs_OVI_c_OVI_c[i] = (data->c_OVI_c_OVI_c[bin_id+1] - data->c_OVI_c_OVI_c[bin_id]);;
        data->dcs_OVI_c_OVI_c[i] /= data->dT[i];
	data->dcs_OVI_c_OVI_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_OVI_ph_OVI_ph[i] = 0.0;
            data->dcs_OVI_ph_OVI_ph[i] = 0.0;
        } else {
            data->cs_OVI_ph_OVI_ph[i] = data->c_OVI_ph_OVI_ph[zbin_id] +
                data->zdef * (data->c_OVI_ph_OVI_ph[zbin_id+1] - data->c_OVI_ph_OVI_ph[zbin_id]);
            data->dcs_OVI_ph_OVI_ph[i] = (data->c_OVI_ph_OVI_ph[zbin_id+1] - data->c_OVI_ph_OVI_ph[zbin_id]);;
            data->dcs_OVI_ph_OVI_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OVII_c_OVII_c[i] = data->c_OVII_c_OVII_c[bin_id] +
            data->Tdef[i] * (data->c_OVII_c_OVII_c[bin_id+1] - data->c_OVII_c_OVII_c[bin_id]);
        data->dcs_OVII_c_OVII_c[i] = (data->c_OVII_c_OVII_c[bin_id+1] - data->c_OVII_c_OVII_c[bin_id]);;
        data->dcs_OVII_c_OVII_c[i] /= data->dT[i];
	data->dcs_OVII_c_OVII_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_OVII_ph_OVII_ph[i] = 0.0;
            data->dcs_OVII_ph_OVII_ph[i] = 0.0;
        } else {
            data->cs_OVII_ph_OVII_ph[i] = data->c_OVII_ph_OVII_ph[zbin_id] +
                data->zdef * (data->c_OVII_ph_OVII_ph[zbin_id+1] - data->c_OVII_ph_OVII_ph[zbin_id]);
            data->dcs_OVII_ph_OVII_ph[i] = (data->c_OVII_ph_OVII_ph[zbin_id+1] - data->c_OVII_ph_OVII_ph[zbin_id]);;
            data->dcs_OVII_ph_OVII_ph[i] /= data->dz;
        }
    }
    
    for (i = 0; i < nstrip; i++) {
        bin_id = data->bin_id[i];
        data->cs_OVIII_c_OVIII_c[i] = data->c_OVIII_c_OVIII_c[bin_id] +
            data->Tdef[i] * (data->c_OVIII_c_OVIII_c[bin_id+1] - data->c_OVIII_c_OVIII_c[bin_id]);
        data->dcs_OVIII_c_OVIII_c[i] = (data->c_OVIII_c_OVIII_c[bin_id+1] - data->c_OVIII_c_OVIII_c[bin_id]);;
        data->dcs_OVIII_c_OVIII_c[i] /= data->dT[i];
	data->dcs_OVIII_c_OVIII_c[i] *= data->invTs[i];
    }
    
    for (i = 0; i < nstrip; i++) {
        if (no_photo) {
            data->cs_OVIII_ph_OVIII_ph[i] = 0.0;
            data->dcs_OVIII_ph_OVIII_ph[i] = 0.0;
        } else {
            data->cs_OVIII_ph_OVIII_ph[i] = data->c_OVIII_ph_OVIII_ph[zbin_id] +
                data->zdef * (data->c_OVIII_ph_OVIII_ph[zbin_id+1] - data->c_OVIII_ph_OVIII_ph[zbin_id]);
            data->dcs_OVIII_ph_OVIII_ph[i] = (data->c_OVIII_ph_OVIII_ph[zbin_id+1] - data->c_OVIII_ph_OVIII_ph[zbin_id]);;
            data->dcs_OVIII_ph_OVIII_ph[i] /= data->dz;
        }
    }
    

}
 



int calculate_rhs_po(double *input, double *rhs, int nstrip,
                  int nchem, void *sdata)
{
    /* We iterate over all of the rates */
    /* Calculate temperature first */
    po_data *data = (po_data*)sdata;
    int i, j;
    po_calculate_temperature(data, input, nstrip, nchem);

    po_interpolate_rates(data, nstrip);

    /* Now we set up some temporaries */
    double *HeI_i = data->rs_HeI_i;
    double *HeI_pi = data->rs_HeI_pi;
    double *HeII_i = data->rs_HeII_i;
    double *HeII_pi = data->rs_HeII_pi;
    double *HeII_r = data->rs_HeII_r;
    double *HeIII_r = data->rs_HeIII_r;
    double *HI_i = data->rs_HI_i;
    double *HI_pi = data->rs_HI_pi;
    double *HII_r = data->rs_HII_r;
    double *OI_i = data->rs_OI_i;
    double *OI_pi = data->rs_OI_pi;
    double *OII_i = data->rs_OII_i;
    double *OII_pi = data->rs_OII_pi;
    double *OII_r = data->rs_OII_r;
    double *OIII_i = data->rs_OIII_i;
    double *OIII_pi = data->rs_OIII_pi;
    double *OIII_r = data->rs_OIII_r;
    double *OIV_i = data->rs_OIV_i;
    double *OIV_pi = data->rs_OIV_pi;
    double *OIV_r = data->rs_OIV_r;
    double *OIX_r = data->rs_OIX_r;
    double *OV_i = data->rs_OV_i;
    double *OV_pi = data->rs_OV_pi;
    double *OV_r = data->rs_OV_r;
    double *OVI_i = data->rs_OVI_i;
    double *OVI_pi = data->rs_OVI_pi;
    double *OVI_r = data->rs_OVI_r;
    double *OVII_i = data->rs_OVII_i;
    double *OVII_pi = data->rs_OVII_pi;
    double *OVII_r = data->rs_OVII_r;
    double *OVIII_i = data->rs_OVIII_i;
    double *OVIII_pi = data->rs_OVIII_pi;
    double *OVIII_r = data->rs_OVIII_r;
    double *compton_comp = data->cs_compton_comp;
    double *HeI_c_HeI_c = data->cs_HeI_c_HeI_c;
    double *HeI_ph_HeI_ph = data->cs_HeI_ph_HeI_ph;
    double *HeII_c_HeII_c = data->cs_HeII_c_HeII_c;
    double *HeII_ph_HeII_ph = data->cs_HeII_ph_HeII_ph;
    double *HeIII_c_HeIII_c = data->cs_HeIII_c_HeIII_c;
    double *HI_c_HI_c = data->cs_HI_c_HI_c;
    double *HI_ph_HI_ph = data->cs_HI_ph_HI_ph;
    double *HII_c_HII_c = data->cs_HII_c_HII_c;
    double *OI_c_OI_c = data->cs_OI_c_OI_c;
    double *OI_ph_OI_ph = data->cs_OI_ph_OI_ph;
    double *OII_c_OII_c = data->cs_OII_c_OII_c;
    double *OII_ph_OII_ph = data->cs_OII_ph_OII_ph;
    double *OIII_c_OIII_c = data->cs_OIII_c_OIII_c;
    double *OIII_ph_OIII_ph = data->cs_OIII_ph_OIII_ph;
    double *OIV_c_OIV_c = data->cs_OIV_c_OIV_c;
    double *OIV_ph_OIV_ph = data->cs_OIV_ph_OIV_ph;
    double *OIX_c_OIX_c = data->cs_OIX_c_OIX_c;
    double *OV_c_OV_c = data->cs_OV_c_OV_c;
    double *OV_ph_OV_ph = data->cs_OV_ph_OV_ph;
    double *OVI_c_OVI_c = data->cs_OVI_c_OVI_c;
    double *OVI_ph_OVI_ph = data->cs_OVI_ph_OVI_ph;
    double *OVII_c_OVII_c = data->cs_OVII_c_OVII_c;
    double *OVII_ph_OVII_ph = data->cs_OVII_ph_OVII_ph;
    double *OVIII_c_OVIII_c = data->cs_OVIII_c_OVIII_c;
    double *OVIII_ph_OVIII_ph = data->cs_OVIII_ph_OVIII_ph;
    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;
    double z;
    double T;

    double mh = 1.67e-24;
    double total, total_e, total_de, mdensity;
    for (i = 0; i<nstrip; i++) {
        j = i * nchem;
        total = total_e = total_de = mdensity = 0.0;
        T = data->Ts[i];
        z = data->current_z;
        HI = input[j];
        if (HI < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HI] = % 0.16g [%d]\n",
               i, HI, j); */
            return 1;
            //HI = 1e-20;
        }
        
        
          total+=HI * 1.00794;
        
        
        j++;
    
        HII = input[j];
        if (HII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HII] = % 0.16g [%d]\n",
               i, HII, j); */
            return 1;
            //HII = 1e-20;
        }
        
        
          total+=HII * 1.00794;
        
        
        j++;
    
        HeI = input[j];
        if (HeI < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HeI] = % 0.16g [%d]\n",
               i, HeI, j); */
            return 1;
            //HeI = 1e-20;
        }
        
        
          total+=HeI * 4.002602;
        
        
        j++;
    
        HeII = input[j];
        if (HeII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HeII] = % 0.16g [%d]\n",
               i, HeII, j); */
            return 1;
            //HeII = 1e-20;
        }
        
        
          total+=HeII * 4.002602;
        
        
        j++;
    
        HeIII = input[j];
        if (HeIII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HeIII] = % 0.16g [%d]\n",
               i, HeIII, j); */
            return 1;
            //HeIII = 1e-20;
        }
        
        
          total+=HeIII * 4.002602;
        
        
        j++;
    
        de = input[j];
        if (de < 0.0) {
            /* fprintf(stderr, "RNegative[%d][de] = % 0.16g [%d]\n",
               i, de, j); */
            return 1;
            //de = 1e-20;
        }
        
        
        
        j++;
    
        ge = input[j];
        if (ge < 0.0) {
            /* fprintf(stderr, "RNegative[%d][ge] = % 0.16g [%d]\n",
               i, ge, j); */
            return 1;
            //ge = 1e-20;
        }
        
        j++;
    
        OI = input[j];
        if (OI < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OI] = % 0.16g [%d]\n",
               i, OI, j); */
            return 1;
            //OI = 1e-20;
        }
        
        
          total+=OI * 15.9994;
        
        
        j++;
    
        OII = input[j];
        if (OII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OII] = % 0.16g [%d]\n",
               i, OII, j); */
            return 1;
            //OII = 1e-20;
        }
        
        
          total+=OII * 15.9994;
        
        
        j++;
    
        OIII = input[j];
        if (OIII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OIII] = % 0.16g [%d]\n",
               i, OIII, j); */
            return 1;
            //OIII = 1e-20;
        }
        
        
          total+=OIII * 15.9994;
        
        
        j++;
    
        OIV = input[j];
        if (OIV < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OIV] = % 0.16g [%d]\n",
               i, OIV, j); */
            return 1;
            //OIV = 1e-20;
        }
        
        
          total+=OIV * 15.9994;
        
        
        j++;
    
        OV = input[j];
        if (OV < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OV] = % 0.16g [%d]\n",
               i, OV, j); */
            return 1;
            //OV = 1e-20;
        }
        
        
          total+=OV * 15.9994;
        
        
        j++;
    
        OVI = input[j];
        if (OVI < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OVI] = % 0.16g [%d]\n",
               i, OVI, j); */
            return 1;
            //OVI = 1e-20;
        }
        
        
          total+=OVI * 15.9994;
        
        
        j++;
    
        OVII = input[j];
        if (OVII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OVII] = % 0.16g [%d]\n",
               i, OVII, j); */
            return 1;
            //OVII = 1e-20;
        }
        
        
          total+=OVII * 15.9994;
        
        
        j++;
    
        OVIII = input[j];
        if (OVIII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OVIII] = % 0.16g [%d]\n",
               i, OVIII, j); */
            return 1;
            //OVIII = 1e-20;
        }
        
        
          total+=OVIII * 15.9994;
        
        
        j++;
    
        OIX = input[j];
        if (OIX < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OIX] = % 0.16g [%d]\n",
               i, OIX, j); */
            return 1;
            //OIX = 1e-20;
        }
        
        
          total+=OIX * 15.9994;
        
        
        j++;
    
        mdensity = total * mh;
        total = 0.0;
        j = i * nchem;
        // 
        // Species: HI
        // 
        rhs[j] = HII_r[i]*HII*de - HI_i[i]*HI*de - HI_pi[i]*HI;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 1.00794;
            total_e += HI * 0.0;
        
        
            total_de += rhs[j] * 0.0;
        
        j++;
    
        // 
        // Species: HII
        // 
        rhs[j] = -HII_r[i]*HII*de + HI_i[i]*HI*de + HI_pi[i]*HI;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 1.00794;
            total_e += HII * 1.0;
        
        
            total_de += rhs[j] * 1.0;
        
        j++;
    
        // 
        // Species: HeI
        // 
        rhs[j] = HeII_r[i]*HeII*de - HeI_i[i]*HeI*de - HeI_pi[i]*HeI;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 4.002602;
            total_e += HeI * 0.0;
        
        
            total_de += rhs[j] * 0.0;
        
        j++;
    
        // 
        // Species: HeII
        // 
        rhs[j] = HeIII_r[i]*HeIII*de - HeII_i[i]*HeII*de - HeII_pi[i]*HeII - HeII_r[i]*HeII*de + HeI_i[i]*HeI*de + HeI_pi[i]*HeI;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 4.002602;
            total_e += HeII * 1.0;
        
        
            total_de += rhs[j] * 1.0;
        
        j++;
    
        // 
        // Species: HeIII
        // 
        rhs[j] = -HeIII_r[i]*HeIII*de + HeII_i[i]*HeII*de + HeII_pi[i]*HeII;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 4.002602;
            total_e += HeIII * 2.0;
        
        
            total_de += rhs[j] * 2.0;
        
        j++;
    
        // 
        // Species: de
        // 
        rhs[j] = -HII_r[i]*HII*de + HI_i[i]*HI*de + HI_pi[i]*HI - HeIII_r[i]*HeIII*de + HeII_i[i]*HeII*de + HeII_pi[i]*HeII - HeII_r[i]*HeII*de + HeI_i[i]*HeI*de + HeI_pi[i]*HeI + OIII_i[i]*OIII*de + OIII_pi[i]*OIII - OIII_r[i]*OIII*de + OII_i[i]*OII*de + OII_pi[i]*OII - OII_r[i]*OII*de + OIV_i[i]*OIV*de + OIV_pi[i]*OIV - OIV_r[i]*OIV*de - OIX_r[i]*OIX*de + OI_i[i]*OI*de + OI_pi[i]*OI + OVIII_i[i]*OVIII*de + OVIII_pi[i]*OVIII - OVIII_r[i]*OVIII*de + OVII_i[i]*OVII*de + OVII_pi[i]*OVII - OVII_r[i]*OVII*de + OVI_i[i]*OVI*de + OVI_pi[i]*OVI - OVI_r[i]*OVI*de + OV_i[i]*OV*de + OV_pi[i]*OV - OV_r[i]*OV*de;
        
        
            total_de += -rhs[j];
        
        j++;
    
        // 
        // Species: ge
        // 
        rhs[j] = -HI*HI_c_HI_c[i]*de + HI*HI_ph_HI_ph[i] - HII*HII_c_HII_c[i]*de - HeI*HeI_c_HeI_c[i]*de + HeI*HeI_ph_HeI_ph[i] - HeII*HeII_c_HeII_c[i]*de + HeII*HeII_ph_HeII_ph[i] - HeIII*HeIII_c_HeIII_c[i]*de - OI*OI_c_OI_c[i]*de + OI*OI_ph_OI_ph[i] - OII*OII_c_OII_c[i]*de + OII*OII_ph_OII_ph[i] - OIII*OIII_c_OIII_c[i]*de + OIII*OIII_ph_OIII_ph[i] - OIV*OIV_c_OIV_c[i]*de + OIV*OIV_ph_OIV_ph[i] - OIX*OIX_c_OIX_c[i]*de - OV*OV_c_OV_c[i]*de + OV*OV_ph_OV_ph[i] - OVI*OVI_c_OVI_c[i]*de + OVI*OVI_ph_OVI_ph[i] - OVII*OVII_c_OVII_c[i]*de + OVII*OVII_ph_OVII_ph[i] - OVIII*OVIII_c_OVIII_c[i]*de + OVIII*OVIII_ph_OVIII_ph[i] - compton_comp[i]*de*pow(z + 1.0, 4)*(T - 2.73*z - 2.73);
        
	    rhs[j] /= mdensity;
        
        
        j++;
    
        // 
        // Species: OI
        // 
        rhs[j] = OII_r[i]*OII*de - OI_i[i]*OI*de - OI_pi[i]*OI;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OI * 0;
        
        
            total_de += rhs[j] * 0;
        
        j++;
    
        // 
        // Species: OII
        // 
        rhs[j] = OIII_r[i]*OIII*de - OII_i[i]*OII*de - OII_pi[i]*OII - OII_r[i]*OII*de + OI_i[i]*OI*de + OI_pi[i]*OI;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OII * 1;
        
        
            total_de += rhs[j] * 1;
        
        j++;
    
        // 
        // Species: OIII
        // 
        rhs[j] = -OIII_i[i]*OIII*de - OIII_pi[i]*OIII - OIII_r[i]*OIII*de + OII_i[i]*OII*de + OII_pi[i]*OII + OIV_r[i]*OIV*de;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OIII * 2;
        
        
            total_de += rhs[j] * 2;
        
        j++;
    
        // 
        // Species: OIV
        // 
        rhs[j] = OIII_i[i]*OIII*de + OIII_pi[i]*OIII - OIV_i[i]*OIV*de - OIV_pi[i]*OIV - OIV_r[i]*OIV*de + OV_r[i]*OV*de;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OIV * 3;
        
        
            total_de += rhs[j] * 3;
        
        j++;
    
        // 
        // Species: OV
        // 
        rhs[j] = OIV_i[i]*OIV*de + OIV_pi[i]*OIV + OVI_r[i]*OVI*de - OV_i[i]*OV*de - OV_pi[i]*OV - OV_r[i]*OV*de;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OV * 4;
        
        
            total_de += rhs[j] * 4;
        
        j++;
    
        // 
        // Species: OVI
        // 
        rhs[j] = OVII_r[i]*OVII*de - OVI_i[i]*OVI*de - OVI_pi[i]*OVI - OVI_r[i]*OVI*de + OV_i[i]*OV*de + OV_pi[i]*OV;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OVI * 5;
        
        
            total_de += rhs[j] * 5;
        
        j++;
    
        // 
        // Species: OVII
        // 
        rhs[j] = OVIII_r[i]*OVIII*de - OVII_i[i]*OVII*de - OVII_pi[i]*OVII - OVII_r[i]*OVII*de + OVI_i[i]*OVI*de + OVI_pi[i]*OVI;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OVII * 6;
        
        
            total_de += rhs[j] * 6;
        
        j++;
    
        // 
        // Species: OVIII
        // 
        rhs[j] = OIX_r[i]*OIX*de - OVIII_i[i]*OVIII*de - OVIII_pi[i]*OVIII - OVIII_r[i]*OVIII*de + OVII_i[i]*OVII*de + OVII_pi[i]*OVII;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OVIII * 7;
        
        
            total_de += rhs[j] * 7;
        
        j++;
    
        // 
        // Species: OIX
        // 
        rhs[j] = -OIX_r[i]*OIX*de + OVIII_i[i]*OVIII*de + OVIII_pi[i]*OVIII;
        
            /* Already in number density, not mass density */
            total += rhs[j] * 15.9994;
            total_e += OIX * 8;
        
        
            total_de += rhs[j] * 8;
        
        j++;
    
    }  
    return 0;
}




int calculate_jacobian_po(double *input, double *Joutput,
        int nstrip, int nchem, void *sdata)
{
    /* We iterate over all of the rates */
    /* Calculate temperature first */
    po_data *data = (po_data*)sdata;

    int i, j;
    //po_calculate_temperature(data, input, nstrip, nchem);

    //po_interpolate_rates(data, nstrip);

    /* Now we set up some temporaries */
    double *Tge = data->dTs_ge;
    double *HeI_i = data->rs_HeI_i;
    double *rHeI_i = data->drs_HeI_i;
    double *HeI_pi = data->rs_HeI_pi;
    double *rHeI_pi = data->drs_HeI_pi;
    double *HeII_i = data->rs_HeII_i;
    double *rHeII_i = data->drs_HeII_i;
    double *HeII_pi = data->rs_HeII_pi;
    double *rHeII_pi = data->drs_HeII_pi;
    double *HeII_r = data->rs_HeII_r;
    double *rHeII_r = data->drs_HeII_r;
    double *HeIII_r = data->rs_HeIII_r;
    double *rHeIII_r = data->drs_HeIII_r;
    double *HI_i = data->rs_HI_i;
    double *rHI_i = data->drs_HI_i;
    double *HI_pi = data->rs_HI_pi;
    double *rHI_pi = data->drs_HI_pi;
    double *HII_r = data->rs_HII_r;
    double *rHII_r = data->drs_HII_r;
    double *OI_i = data->rs_OI_i;
    double *rOI_i = data->drs_OI_i;
    double *OI_pi = data->rs_OI_pi;
    double *rOI_pi = data->drs_OI_pi;
    double *OII_i = data->rs_OII_i;
    double *rOII_i = data->drs_OII_i;
    double *OII_pi = data->rs_OII_pi;
    double *rOII_pi = data->drs_OII_pi;
    double *OII_r = data->rs_OII_r;
    double *rOII_r = data->drs_OII_r;
    double *OIII_i = data->rs_OIII_i;
    double *rOIII_i = data->drs_OIII_i;
    double *OIII_pi = data->rs_OIII_pi;
    double *rOIII_pi = data->drs_OIII_pi;
    double *OIII_r = data->rs_OIII_r;
    double *rOIII_r = data->drs_OIII_r;
    double *OIV_i = data->rs_OIV_i;
    double *rOIV_i = data->drs_OIV_i;
    double *OIV_pi = data->rs_OIV_pi;
    double *rOIV_pi = data->drs_OIV_pi;
    double *OIV_r = data->rs_OIV_r;
    double *rOIV_r = data->drs_OIV_r;
    double *OIX_r = data->rs_OIX_r;
    double *rOIX_r = data->drs_OIX_r;
    double *OV_i = data->rs_OV_i;
    double *rOV_i = data->drs_OV_i;
    double *OV_pi = data->rs_OV_pi;
    double *rOV_pi = data->drs_OV_pi;
    double *OV_r = data->rs_OV_r;
    double *rOV_r = data->drs_OV_r;
    double *OVI_i = data->rs_OVI_i;
    double *rOVI_i = data->drs_OVI_i;
    double *OVI_pi = data->rs_OVI_pi;
    double *rOVI_pi = data->drs_OVI_pi;
    double *OVI_r = data->rs_OVI_r;
    double *rOVI_r = data->drs_OVI_r;
    double *OVII_i = data->rs_OVII_i;
    double *rOVII_i = data->drs_OVII_i;
    double *OVII_pi = data->rs_OVII_pi;
    double *rOVII_pi = data->drs_OVII_pi;
    double *OVII_r = data->rs_OVII_r;
    double *rOVII_r = data->drs_OVII_r;
    double *OVIII_i = data->rs_OVIII_i;
    double *rOVIII_i = data->drs_OVIII_i;
    double *OVIII_pi = data->rs_OVIII_pi;
    double *rOVIII_pi = data->drs_OVIII_pi;
    double *OVIII_r = data->rs_OVIII_r;
    double *rOVIII_r = data->drs_OVIII_r;
    double *compton_comp = data->cs_compton_comp;
    double *rcompton_comp = data->dcs_compton_comp;
    double *HeI_c_HeI_c = data->cs_HeI_c_HeI_c;
    double *rHeI_c_HeI_c = data->dcs_HeI_c_HeI_c;
    double *HeI_ph_HeI_ph = data->cs_HeI_ph_HeI_ph;
    double *rHeI_ph_HeI_ph = data->dcs_HeI_ph_HeI_ph;
    double *HeII_c_HeII_c = data->cs_HeII_c_HeII_c;
    double *rHeII_c_HeII_c = data->dcs_HeII_c_HeII_c;
    double *HeII_ph_HeII_ph = data->cs_HeII_ph_HeII_ph;
    double *rHeII_ph_HeII_ph = data->dcs_HeII_ph_HeII_ph;
    double *HeIII_c_HeIII_c = data->cs_HeIII_c_HeIII_c;
    double *rHeIII_c_HeIII_c = data->dcs_HeIII_c_HeIII_c;
    double *HI_c_HI_c = data->cs_HI_c_HI_c;
    double *rHI_c_HI_c = data->dcs_HI_c_HI_c;
    double *HI_ph_HI_ph = data->cs_HI_ph_HI_ph;
    double *rHI_ph_HI_ph = data->dcs_HI_ph_HI_ph;
    double *HII_c_HII_c = data->cs_HII_c_HII_c;
    double *rHII_c_HII_c = data->dcs_HII_c_HII_c;
    double *OI_c_OI_c = data->cs_OI_c_OI_c;
    double *rOI_c_OI_c = data->dcs_OI_c_OI_c;
    double *OI_ph_OI_ph = data->cs_OI_ph_OI_ph;
    double *rOI_ph_OI_ph = data->dcs_OI_ph_OI_ph;
    double *OII_c_OII_c = data->cs_OII_c_OII_c;
    double *rOII_c_OII_c = data->dcs_OII_c_OII_c;
    double *OII_ph_OII_ph = data->cs_OII_ph_OII_ph;
    double *rOII_ph_OII_ph = data->dcs_OII_ph_OII_ph;
    double *OIII_c_OIII_c = data->cs_OIII_c_OIII_c;
    double *rOIII_c_OIII_c = data->dcs_OIII_c_OIII_c;
    double *OIII_ph_OIII_ph = data->cs_OIII_ph_OIII_ph;
    double *rOIII_ph_OIII_ph = data->dcs_OIII_ph_OIII_ph;
    double *OIV_c_OIV_c = data->cs_OIV_c_OIV_c;
    double *rOIV_c_OIV_c = data->dcs_OIV_c_OIV_c;
    double *OIV_ph_OIV_ph = data->cs_OIV_ph_OIV_ph;
    double *rOIV_ph_OIV_ph = data->dcs_OIV_ph_OIV_ph;
    double *OIX_c_OIX_c = data->cs_OIX_c_OIX_c;
    double *rOIX_c_OIX_c = data->dcs_OIX_c_OIX_c;
    double *OV_c_OV_c = data->cs_OV_c_OV_c;
    double *rOV_c_OV_c = data->dcs_OV_c_OV_c;
    double *OV_ph_OV_ph = data->cs_OV_ph_OV_ph;
    double *rOV_ph_OV_ph = data->dcs_OV_ph_OV_ph;
    double *OVI_c_OVI_c = data->cs_OVI_c_OVI_c;
    double *rOVI_c_OVI_c = data->dcs_OVI_c_OVI_c;
    double *OVI_ph_OVI_ph = data->cs_OVI_ph_OVI_ph;
    double *rOVI_ph_OVI_ph = data->dcs_OVI_ph_OVI_ph;
    double *OVII_c_OVII_c = data->cs_OVII_c_OVII_c;
    double *rOVII_c_OVII_c = data->dcs_OVII_c_OVII_c;
    double *OVII_ph_OVII_ph = data->cs_OVII_ph_OVII_ph;
    double *rOVII_ph_OVII_ph = data->dcs_OVII_ph_OVII_ph;
    double *OVIII_c_OVIII_c = data->cs_OVIII_c_OVIII_c;
    double *rOVIII_c_OVIII_c = data->dcs_OVIII_c_OVIII_c;
    double *OVIII_ph_OVIII_ph = data->cs_OVIII_ph_OVIII_ph;
    double *rOVIII_ph_OVIII_ph = data->dcs_OVIII_ph_OVIII_ph;
    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;
    double z;
    double T;

    double mh = 1.67e-24;
    double total, mdensity;
    for (i = 0; i<nstrip; i++) {
        j = i * nchem;
	total = mdensity = 0.0;
        T = data->Ts[i];
        z = data->current_z;
	    HI = input[j];
        if (HI < 0.0) {
          fprintf(stderr, "JNegative[%d][HI] = % 0.16g [%d]\n",
            i, HI, j);
          /*HI = 0.0;*/
          //HI = 1e-20;
          return 1;
        }
	
        
          total+=HI * 1.00794;
        
	
        j++;
    
	    HII = input[j];
        if (HII < 0.0) {
          fprintf(stderr, "JNegative[%d][HII] = % 0.16g [%d]\n",
            i, HII, j);
          /*HII = 0.0;*/
          //HII = 1e-20;
          return 1;
        }
	
        
          total+=HII * 1.00794;
        
	
        j++;
    
	    HeI = input[j];
        if (HeI < 0.0) {
          fprintf(stderr, "JNegative[%d][HeI] = % 0.16g [%d]\n",
            i, HeI, j);
          /*HeI = 0.0;*/
          //HeI = 1e-20;
          return 1;
        }
	
        
          total+=HeI * 4.002602;
        
	
        j++;
    
	    HeII = input[j];
        if (HeII < 0.0) {
          fprintf(stderr, "JNegative[%d][HeII] = % 0.16g [%d]\n",
            i, HeII, j);
          /*HeII = 0.0;*/
          //HeII = 1e-20;
          return 1;
        }
	
        
          total+=HeII * 4.002602;
        
	
        j++;
    
	    HeIII = input[j];
        if (HeIII < 0.0) {
          fprintf(stderr, "JNegative[%d][HeIII] = % 0.16g [%d]\n",
            i, HeIII, j);
          /*HeIII = 0.0;*/
          //HeIII = 1e-20;
          return 1;
        }
	
        
          total+=HeIII * 4.002602;
        
	
        j++;
    
	    de = input[j];
        if (de < 0.0) {
          fprintf(stderr, "JNegative[%d][de] = % 0.16g [%d]\n",
            i, de, j);
          /*de = 0.0;*/
          //de = 1e-20;
          return 1;
        }
	
        
	
        j++;
    
	    ge = input[j];
        if (ge < 0.0) {
          fprintf(stderr, "JNegative[%d][ge] = % 0.16g [%d]\n",
            i, ge, j);
          /*ge = 0.0;*/
          //ge = 1e-20;
          return 1;
        }
	
        j++;
    
	    OI = input[j];
        if (OI < 0.0) {
          fprintf(stderr, "JNegative[%d][OI] = % 0.16g [%d]\n",
            i, OI, j);
          /*OI = 0.0;*/
          //OI = 1e-20;
          return 1;
        }
	
        
          total+=OI * 15.9994;
        
	
        j++;
    
	    OII = input[j];
        if (OII < 0.0) {
          fprintf(stderr, "JNegative[%d][OII] = % 0.16g [%d]\n",
            i, OII, j);
          /*OII = 0.0;*/
          //OII = 1e-20;
          return 1;
        }
	
        
          total+=OII * 15.9994;
        
	
        j++;
    
	    OIII = input[j];
        if (OIII < 0.0) {
          fprintf(stderr, "JNegative[%d][OIII] = % 0.16g [%d]\n",
            i, OIII, j);
          /*OIII = 0.0;*/
          //OIII = 1e-20;
          return 1;
        }
	
        
          total+=OIII * 15.9994;
        
	
        j++;
    
	    OIV = input[j];
        if (OIV < 0.0) {
          fprintf(stderr, "JNegative[%d][OIV] = % 0.16g [%d]\n",
            i, OIV, j);
          /*OIV = 0.0;*/
          //OIV = 1e-20;
          return 1;
        }
	
        
          total+=OIV * 15.9994;
        
	
        j++;
    
	    OV = input[j];
        if (OV < 0.0) {
          fprintf(stderr, "JNegative[%d][OV] = % 0.16g [%d]\n",
            i, OV, j);
          /*OV = 0.0;*/
          //OV = 1e-20;
          return 1;
        }
	
        
          total+=OV * 15.9994;
        
	
        j++;
    
	    OVI = input[j];
        if (OVI < 0.0) {
          fprintf(stderr, "JNegative[%d][OVI] = % 0.16g [%d]\n",
            i, OVI, j);
          /*OVI = 0.0;*/
          //OVI = 1e-20;
          return 1;
        }
	
        
          total+=OVI * 15.9994;
        
	
        j++;
    
	    OVII = input[j];
        if (OVII < 0.0) {
          fprintf(stderr, "JNegative[%d][OVII] = % 0.16g [%d]\n",
            i, OVII, j);
          /*OVII = 0.0;*/
          //OVII = 1e-20;
          return 1;
        }
	
        
          total+=OVII * 15.9994;
        
	
        j++;
    
	    OVIII = input[j];
        if (OVIII < 0.0) {
          fprintf(stderr, "JNegative[%d][OVIII] = % 0.16g [%d]\n",
            i, OVIII, j);
          /*OVIII = 0.0;*/
          //OVIII = 1e-20;
          return 1;
        }
	
        
          total+=OVIII * 15.9994;
        
	
        j++;
    
	    OIX = input[j];
        if (OIX < 0.0) {
          fprintf(stderr, "JNegative[%d][OIX] = % 0.16g [%d]\n",
            i, OIX, j);
          /*OIX = 0.0;*/
          //OIX = 1e-20;
          return 1;
        }
	
        
          total+=OIX * 15.9994;
        
	
        j++;
    
        mdensity = total * mh;
        
        j = i * nchem * nchem;
        // 
        // Species: HI
        //
            // HI by HI
            Joutput[j] = -HI_i[i]*de - HI_pi[i];
	    
	    
            j++;
            // HII by HI
            Joutput[j] = HI_i[i]*de + HI_pi[i];
	    
	    
            j++;
            // HeI by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // de by HI
            Joutput[j] = HI_i[i]*de + HI_pi[i];
	    
	    
            j++;
            // ge by HI
            Joutput[j] = -HI_c_HI_c[i]*de + HI_ph_HI_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by HI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by HI
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: HII
        //
            // HI by HII
            Joutput[j] = HII_r[i]*de;
	    
	    
            j++;
            // HII by HII
            Joutput[j] = -HII_r[i]*de;
	    
	    
            j++;
            // HeI by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // de by HII
            Joutput[j] = -HII_r[i]*de;
	    
	    
            j++;
            // ge by HII
            Joutput[j] = -HII_c_HII_c[i]*de;
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by HII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by HII
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: HeI
        //
            // HI by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by HeI
            Joutput[j] = -HeI_i[i]*de - HeI_pi[i];
	    
	    
            j++;
            // HeII by HeI
            Joutput[j] = HeI_i[i]*de + HeI_pi[i];
	    
	    
            j++;
            // HeIII by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // de by HeI
            Joutput[j] = HeI_i[i]*de + HeI_pi[i];
	    
	    
            j++;
            // ge by HeI
            Joutput[j] = -HeI_c_HeI_c[i]*de + HeI_ph_HeI_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by HeI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by HeI
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: HeII
        //
            // HI by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by HeII
            Joutput[j] = HeII_r[i]*de;
	    
	    
            j++;
            // HeII by HeII
            Joutput[j] = -HeII_i[i]*de - HeII_pi[i] - HeII_r[i]*de;
	    
	    
            j++;
            // HeIII by HeII
            Joutput[j] = HeII_i[i]*de + HeII_pi[i];
	    
	    
            j++;
            // de by HeII
            Joutput[j] = HeII_i[i]*de + HeII_pi[i] - HeII_r[i]*de;
	    
	    
            j++;
            // ge by HeII
            Joutput[j] = -HeII_c_HeII_c[i]*de + HeII_ph_HeII_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by HeII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by HeII
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: HeIII
        //
            // HI by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by HeIII
            Joutput[j] = HeIII_r[i]*de;
	    
	    
            j++;
            // HeIII by HeIII
            Joutput[j] = -HeIII_r[i]*de;
	    
	    
            j++;
            // de by HeIII
            Joutput[j] = -HeIII_r[i]*de;
	    
	    
            j++;
            // ge by HeIII
            Joutput[j] = -HeIII_c_HeIII_c[i]*de;
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by HeIII
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: de
        //
            // HI by de
            Joutput[j] = HII_r[i]*HII - HI_i[i]*HI;
	    
	    
            j++;
            // HII by de
            Joutput[j] = -HII_r[i]*HII + HI_i[i]*HI;
	    
	    
            j++;
            // HeI by de
            Joutput[j] = HeII_r[i]*HeII - HeI_i[i]*HeI;
	    
	    
            j++;
            // HeII by de
            Joutput[j] = HeIII_r[i]*HeIII - HeII_i[i]*HeII - HeII_r[i]*HeII + HeI_i[i]*HeI;
	    
	    
            j++;
            // HeIII by de
            Joutput[j] = -HeIII_r[i]*HeIII + HeII_i[i]*HeII;
	    
	    
            j++;
            // de by de
            Joutput[j] = -HII_r[i]*HII + HI_i[i]*HI - HeIII_r[i]*HeIII + HeII_i[i]*HeII - HeII_r[i]*HeII + HeI_i[i]*HeI + OIII_i[i]*OIII - OIII_r[i]*OIII + OII_i[i]*OII - OII_r[i]*OII + OIV_i[i]*OIV - OIV_r[i]*OIV - OIX_r[i]*OIX + OI_i[i]*OI + OVIII_i[i]*OVIII - OVIII_r[i]*OVIII + OVII_i[i]*OVII - OVII_r[i]*OVII + OVI_i[i]*OVI - OVI_r[i]*OVI + OV_i[i]*OV - OV_r[i]*OV;
	    
	    
            j++;
            // ge by de
            Joutput[j] = -HI*HI_c_HI_c[i] - HII*HII_c_HII_c[i] - HeI*HeI_c_HeI_c[i] - HeII*HeII_c_HeII_c[i] - HeIII*HeIII_c_HeIII_c[i] - OI*OI_c_OI_c[i] - OII*OII_c_OII_c[i] - OIII*OIII_c_OIII_c[i] - OIV*OIV_c_OIV_c[i] - OIX*OIX_c_OIX_c[i] - OV*OV_c_OV_c[i] - OVI*OVI_c_OVI_c[i] - OVII*OVII_c_OVII_c[i] - OVIII*OVIII_c_OVIII_c[i] - compton_comp[i]*pow(z + 1.0, 4)*(T - 2.73*z - 2.73);
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by de
            Joutput[j] = OII_r[i]*OII - OI_i[i]*OI;
	    
	    
            j++;
            // OII by de
            Joutput[j] = OIII_r[i]*OIII - OII_i[i]*OII - OII_r[i]*OII + OI_i[i]*OI;
	    
	    
            j++;
            // OIII by de
            Joutput[j] = -OIII_i[i]*OIII - OIII_r[i]*OIII + OII_i[i]*OII + OIV_r[i]*OIV;
	    
	    
            j++;
            // OIV by de
            Joutput[j] = OIII_i[i]*OIII - OIV_i[i]*OIV - OIV_r[i]*OIV + OV_r[i]*OV;
	    
	    
            j++;
            // OV by de
            Joutput[j] = OIV_i[i]*OIV + OVI_r[i]*OVI - OV_i[i]*OV - OV_r[i]*OV;
	    
	    
            j++;
            // OVI by de
            Joutput[j] = OVII_r[i]*OVII - OVI_i[i]*OVI - OVI_r[i]*OVI + OV_i[i]*OV;
	    
	    
            j++;
            // OVII by de
            Joutput[j] = OVIII_r[i]*OVIII - OVII_i[i]*OVII - OVII_r[i]*OVII + OVI_i[i]*OVI;
	    
	    
            j++;
            // OVIII by de
            Joutput[j] = OIX_r[i]*OIX - OVIII_i[i]*OVIII - OVIII_r[i]*OVIII + OVII_i[i]*OVII;
	    
	    
            j++;
            // OIX by de
            Joutput[j] = -OIX_r[i]*OIX + OVIII_i[i]*OVIII;
	    
	    
            j++;
    
        // 
        // Species: ge
        //
            // HI by ge
            Joutput[j] = -HI*de*rHI_i[i] - HI*rHI_pi[i] + HII*de*rHII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // HII by ge
            Joutput[j] = HI*de*rHI_i[i] + HI*rHI_pi[i] - HII*de*rHII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // HeI by ge
            Joutput[j] = -HeI*de*rHeI_i[i] - HeI*rHeI_pi[i] + HeII*de*rHeII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // HeII by ge
            Joutput[j] = HeI*de*rHeI_i[i] + HeI*rHeI_pi[i] - HeII*de*rHeII_i[i] - HeII*de*rHeII_r[i] - HeII*rHeII_pi[i] + HeIII*de*rHeIII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // HeIII by ge
            Joutput[j] = HeII*de*rHeII_i[i] + HeII*rHeII_pi[i] - HeIII*de*rHeIII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // de by ge
            Joutput[j] = HI*de*rHI_i[i] + HI*rHI_pi[i] - HII*de*rHII_r[i] + HeI*de*rHeI_i[i] + HeI*rHeI_pi[i] + HeII*de*rHeII_i[i] - HeII*de*rHeII_r[i] + HeII*rHeII_pi[i] - HeIII*de*rHeIII_r[i] + OI*de*rOI_i[i] + OI*rOI_pi[i] + OII*de*rOII_i[i] - OII*de*rOII_r[i] + OII*rOII_pi[i] + OIII*de*rOIII_i[i] - OIII*de*rOIII_r[i] + OIII*rOIII_pi[i] + OIV*de*rOIV_i[i] - OIV*de*rOIV_r[i] + OIV*rOIV_pi[i] - OIX*de*rOIX_r[i] + OV*de*rOV_i[i] - OV*de*rOV_r[i] + OV*rOV_pi[i] + OVI*de*rOVI_i[i] - OVI*de*rOVI_r[i] + OVI*rOVI_pi[i] + OVII*de*rOVII_i[i] - OVII*de*rOVII_r[i] + OVII*rOVII_pi[i] + OVIII*de*rOVIII_i[i] - OVIII*de*rOVIII_r[i] + OVIII*rOVIII_pi[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // ge by ge
            Joutput[j] = 0;
	    
	    Joutput[j] /= mdensity;
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OI by ge
            Joutput[j] = -OI*de*rOI_i[i] - OI*rOI_pi[i] + OII*de*rOII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OII by ge
            Joutput[j] = OI*de*rOI_i[i] + OI*rOI_pi[i] - OII*de*rOII_i[i] - OII*de*rOII_r[i] - OII*rOII_pi[i] + OIII*de*rOIII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OIII by ge
            Joutput[j] = OII*de*rOII_i[i] + OII*rOII_pi[i] - OIII*de*rOIII_i[i] - OIII*de*rOIII_r[i] - OIII*rOIII_pi[i] + OIV*de*rOIV_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OIV by ge
            Joutput[j] = OIII*de*rOIII_i[i] + OIII*rOIII_pi[i] - OIV*de*rOIV_i[i] - OIV*de*rOIV_r[i] - OIV*rOIV_pi[i] + OV*de*rOV_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OV by ge
            Joutput[j] = OIV*de*rOIV_i[i] + OIV*rOIV_pi[i] - OV*de*rOV_i[i] - OV*de*rOV_r[i] - OV*rOV_pi[i] + OVI*de*rOVI_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OVI by ge
            Joutput[j] = OV*de*rOV_i[i] + OV*rOV_pi[i] - OVI*de*rOVI_i[i] - OVI*de*rOVI_r[i] - OVI*rOVI_pi[i] + OVII*de*rOVII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OVII by ge
            Joutput[j] = OVI*de*rOVI_i[i] + OVI*rOVI_pi[i] - OVII*de*rOVII_i[i] - OVII*de*rOVII_r[i] - OVII*rOVII_pi[i] + OVIII*de*rOVIII_r[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OVIII by ge
            Joutput[j] = OIX*de*rOIX_r[i] + OVII*de*rOVII_i[i] + OVII*rOVII_pi[i] - OVIII*de*rOVIII_i[i] - OVIII*de*rOVIII_r[i] - OVIII*rOVIII_pi[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
            // OIX by ge
            Joutput[j] = -OIX*de*rOIX_r[i] + OVIII*de*rOVIII_i[i] + OVIII*rOVIII_pi[i];
	    
	    
            Joutput[j] *= Tge[i];
            
            j++;
    
        // 
        // Species: OI
        //
            // HI by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OI
            Joutput[j] = OI_i[i]*de + OI_pi[i];
	    
	    
            j++;
            // ge by OI
            Joutput[j] = -OI_c_OI_c[i]*de + OI_ph_OI_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OI
            Joutput[j] = -OI_i[i]*de - OI_pi[i];
	    
	    
            j++;
            // OII by OI
            Joutput[j] = OI_i[i]*de + OI_pi[i];
	    
	    
            j++;
            // OIII by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by OI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by OI
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: OII
        //
            // HI by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OII
            Joutput[j] = OII_i[i]*de + OII_pi[i] - OII_r[i]*de;
	    
	    
            j++;
            // ge by OII
            Joutput[j] = -OII_c_OII_c[i]*de + OII_ph_OII_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OII
            Joutput[j] = OII_r[i]*de;
	    
	    
            j++;
            // OII by OII
            Joutput[j] = -OII_i[i]*de - OII_pi[i] - OII_r[i]*de;
	    
	    
            j++;
            // OIII by OII
            Joutput[j] = OII_i[i]*de + OII_pi[i];
	    
	    
            j++;
            // OIV by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by OII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by OII
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: OIII
        //
            // HI by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OIII
            Joutput[j] = OIII_i[i]*de + OIII_pi[i] - OIII_r[i]*de;
	    
	    
            j++;
            // ge by OIII
            Joutput[j] = -OIII_c_OIII_c[i]*de + OIII_ph_OIII_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by OIII
            Joutput[j] = OIII_r[i]*de;
	    
	    
            j++;
            // OIII by OIII
            Joutput[j] = -OIII_i[i]*de - OIII_pi[i] - OIII_r[i]*de;
	    
	    
            j++;
            // OIV by OIII
            Joutput[j] = OIII_i[i]*de + OIII_pi[i];
	    
	    
            j++;
            // OV by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by OIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by OIII
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: OIV
        //
            // HI by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OIV
            Joutput[j] = OIV_i[i]*de + OIV_pi[i] - OIV_r[i]*de;
	    
	    
            j++;
            // ge by OIV
            Joutput[j] = -OIV_c_OIV_c[i]*de + OIV_ph_OIV_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by OIV
            Joutput[j] = OIV_r[i]*de;
	    
	    
            j++;
            // OIV by OIV
            Joutput[j] = -OIV_i[i]*de - OIV_pi[i] - OIV_r[i]*de;
	    
	    
            j++;
            // OV by OIV
            Joutput[j] = OIV_i[i]*de + OIV_pi[i];
	    
	    
            j++;
            // OVI by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by OIV
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by OIV
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: OV
        //
            // HI by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OV
            Joutput[j] = OV_i[i]*de + OV_pi[i] - OV_r[i]*de;
	    
	    
            j++;
            // ge by OV
            Joutput[j] = -OV_c_OV_c[i]*de + OV_ph_OV_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by OV
            Joutput[j] = OV_r[i]*de;
	    
	    
            j++;
            // OV by OV
            Joutput[j] = -OV_i[i]*de - OV_pi[i] - OV_r[i]*de;
	    
	    
            j++;
            // OVI by OV
            Joutput[j] = OV_i[i]*de + OV_pi[i];
	    
	    
            j++;
            // OVII by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by OV
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by OV
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: OVI
        //
            // HI by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OVI
            Joutput[j] = OVI_i[i]*de + OVI_pi[i] - OVI_r[i]*de;
	    
	    
            j++;
            // ge by OVI
            Joutput[j] = -OVI_c_OVI_c[i]*de + OVI_ph_OVI_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by OVI
            Joutput[j] = OVI_r[i]*de;
	    
	    
            j++;
            // OVI by OVI
            Joutput[j] = -OVI_i[i]*de - OVI_pi[i] - OVI_r[i]*de;
	    
	    
            j++;
            // OVII by OVI
            Joutput[j] = OVI_i[i]*de + OVI_pi[i];
	    
	    
            j++;
            // OVIII by OVI
            Joutput[j] = 0;
	    
	    
            j++;
            // OIX by OVI
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: OVII
        //
            // HI by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OVII
            Joutput[j] = OVII_i[i]*de + OVII_pi[i] - OVII_r[i]*de;
	    
	    
            j++;
            // ge by OVII
            Joutput[j] = -OVII_c_OVII_c[i]*de + OVII_ph_OVII_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by OVII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by OVII
            Joutput[j] = OVII_r[i]*de;
	    
	    
            j++;
            // OVII by OVII
            Joutput[j] = -OVII_i[i]*de - OVII_pi[i] - OVII_r[i]*de;
	    
	    
            j++;
            // OVIII by OVII
            Joutput[j] = OVII_i[i]*de + OVII_pi[i];
	    
	    
            j++;
            // OIX by OVII
            Joutput[j] = 0;
	    
	    
            j++;
    
        // 
        // Species: OVIII
        //
            // HI by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OVIII
            Joutput[j] = OVIII_i[i]*de + OVIII_pi[i] - OVIII_r[i]*de;
	    
	    
            j++;
            // ge by OVIII
            Joutput[j] = -OVIII_c_OVIII_c[i]*de + OVIII_ph_OVIII_ph[i];
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by OVIII
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by OVIII
            Joutput[j] = OVIII_r[i]*de;
	    
	    
            j++;
            // OVIII by OVIII
            Joutput[j] = -OVIII_i[i]*de - OVIII_pi[i] - OVIII_r[i]*de;
	    
	    
            j++;
            // OIX by OVIII
            Joutput[j] = OVIII_i[i]*de + OVIII_pi[i];
	    
	    
            j++;
    
        // 
        // Species: OIX
        //
            // HI by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // HII by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // HeI by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // HeII by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // HeIII by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // de by OIX
            Joutput[j] = -OIX_r[i]*de;
	    
	    
            j++;
            // ge by OIX
            Joutput[j] = -OIX_c_OIX_c[i]*de;
	    
	    Joutput[j] /= mdensity;
	    
	    
            j++;
            // OI by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // OII by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // OIII by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // OIV by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // OV by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // OVI by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // OVII by OIX
            Joutput[j] = 0;
	    
	    
            j++;
            // OVIII by OIX
            Joutput[j] = OIX_r[i]*de;
	    
	    
            j++;
            // OIX by OIX
            Joutput[j] = -OIX_r[i]*de;
	    
	    
            j++;
    
    }

    return 0;
    
}


int calculate_cooling_rate_po(double *input, int nstrip, int nchem,
                              double z_current, double *strip_edot, po_data *data)
{
    int i, j;

    data->current_z = z_current;

    po_calculate_temperature(data, input, nstrip, nchem);

    po_interpolate_rates(data, nstrip);

    /* Now we set up some temporaries */
    double *HeI_i = data->rs_HeI_i;
    double *HeI_pi = data->rs_HeI_pi;
    double *HeII_i = data->rs_HeII_i;
    double *HeII_pi = data->rs_HeII_pi;
    double *HeII_r = data->rs_HeII_r;
    double *HeIII_r = data->rs_HeIII_r;
    double *HI_i = data->rs_HI_i;
    double *HI_pi = data->rs_HI_pi;
    double *HII_r = data->rs_HII_r;
    double *OI_i = data->rs_OI_i;
    double *OI_pi = data->rs_OI_pi;
    double *OII_i = data->rs_OII_i;
    double *OII_pi = data->rs_OII_pi;
    double *OII_r = data->rs_OII_r;
    double *OIII_i = data->rs_OIII_i;
    double *OIII_pi = data->rs_OIII_pi;
    double *OIII_r = data->rs_OIII_r;
    double *OIV_i = data->rs_OIV_i;
    double *OIV_pi = data->rs_OIV_pi;
    double *OIV_r = data->rs_OIV_r;
    double *OIX_r = data->rs_OIX_r;
    double *OV_i = data->rs_OV_i;
    double *OV_pi = data->rs_OV_pi;
    double *OV_r = data->rs_OV_r;
    double *OVI_i = data->rs_OVI_i;
    double *OVI_pi = data->rs_OVI_pi;
    double *OVI_r = data->rs_OVI_r;
    double *OVII_i = data->rs_OVII_i;
    double *OVII_pi = data->rs_OVII_pi;
    double *OVII_r = data->rs_OVII_r;
    double *OVIII_i = data->rs_OVIII_i;
    double *OVIII_pi = data->rs_OVIII_pi;
    double *OVIII_r = data->rs_OVIII_r;
    double *compton_comp = data->cs_compton_comp;
    double *HeI_c_HeI_c = data->cs_HeI_c_HeI_c;
    double *HeI_ph_HeI_ph = data->cs_HeI_ph_HeI_ph;
    double *HeII_c_HeII_c = data->cs_HeII_c_HeII_c;
    double *HeII_ph_HeII_ph = data->cs_HeII_ph_HeII_ph;
    double *HeIII_c_HeIII_c = data->cs_HeIII_c_HeIII_c;
    double *HI_c_HI_c = data->cs_HI_c_HI_c;
    double *HI_ph_HI_ph = data->cs_HI_ph_HI_ph;
    double *HII_c_HII_c = data->cs_HII_c_HII_c;
    double *OI_c_OI_c = data->cs_OI_c_OI_c;
    double *OI_ph_OI_ph = data->cs_OI_ph_OI_ph;
    double *OII_c_OII_c = data->cs_OII_c_OII_c;
    double *OII_ph_OII_ph = data->cs_OII_ph_OII_ph;
    double *OIII_c_OIII_c = data->cs_OIII_c_OIII_c;
    double *OIII_ph_OIII_ph = data->cs_OIII_ph_OIII_ph;
    double *OIV_c_OIV_c = data->cs_OIV_c_OIV_c;
    double *OIV_ph_OIV_ph = data->cs_OIV_ph_OIV_ph;
    double *OIX_c_OIX_c = data->cs_OIX_c_OIX_c;
    double *OV_c_OV_c = data->cs_OV_c_OV_c;
    double *OV_ph_OV_ph = data->cs_OV_ph_OV_ph;
    double *OVI_c_OVI_c = data->cs_OVI_c_OVI_c;
    double *OVI_ph_OVI_ph = data->cs_OVI_ph_OVI_ph;
    double *OVII_c_OVII_c = data->cs_OVII_c_OVII_c;
    double *OVII_ph_OVII_ph = data->cs_OVII_ph_OVII_ph;
    double *OVIII_c_OVIII_c = data->cs_OVIII_c_OVIII_c;
    double *OVIII_ph_OVIII_ph = data->cs_OVIII_ph_OVIII_ph;
    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;
    double z;
    double T;

    double mh = 1.67e-24;
    double total, total_e, total_de, mdensity;
    for (i = 0; i<nstrip; i++) {
        j = i * nchem;
        total = total_e = total_de = mdensity = 0.0;
        T = data->Ts[i];
        z = data->current_z;
        HI = input[j];
        if (HI < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HI] = % 0.16g [%d]\n",
               i, HI, j); */
            return 1;
            //HI = 1e-20;
        }        
        
        total+=HI * 1.00794;
        
        j++;
    
        HII = input[j];
        if (HII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HII] = % 0.16g [%d]\n",
               i, HII, j); */
            return 1;
            //HII = 1e-20;
        }        
        
        total+=HII * 1.00794;     
        
        j++;
    
        HeI = input[j];
        if (HeI < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HeI] = % 0.16g [%d]\n",
               i, HeI, j); */
            return 1;
            //HeI = 1e-20;
        }       
        
        total+=HeI * 4.002602;       
        
        j++;
    
        HeII = input[j];
        if (HeII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HeII] = % 0.16g [%d]\n",
               i, HeII, j); */
            return 1;
            //HeII = 1e-20;
        }        
        
        total+=HeII * 4.002602;        
        
        j++;
    
        HeIII = input[j];
        if (HeIII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][HeIII] = % 0.16g [%d]\n",
               i, HeIII, j); */
            return 1;
            //HeIII = 1e-20;
        }        
        
        total+=HeIII * 4.002602;        
        
        j++;
    
        de = input[j];
        if (de < 0.0) {
            /* fprintf(stderr, "RNegative[%d][de] = % 0.16g [%d]\n",
               i, de, j); */
            return 1;
            //de = 1e-20;
        }                
        
        j++;
    
        ge = input[j];
        if (ge < 0.0) {
            /* fprintf(stderr, "RNegative[%d][ge] = % 0.16g [%d]\n",
               i, ge, j); */
            return 1;
            //ge = 1e-20;
        }
        
        j++;
    
        OI = input[j];
        if (OI < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OI] = % 0.16g [%d]\n",
               i, OI, j); */
            return 1;
            //OI = 1e-20;
        }        
        
        total+=OI * 15.9994;        
        
        j++;
    
        OII = input[j];
        if (OII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OII] = % 0.16g [%d]\n",
               i, OII, j); */
            return 1;
            //OII = 1e-20;
        }        
        
        total+=OII * 15.9994;        
        
        j++;
    
        OIII = input[j];
        if (OIII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OIII] = % 0.16g [%d]\n",
               i, OIII, j); */
            return 1;
            //OIII = 1e-20;
        }        
        
        total+=OIII * 15.9994;        
        
        j++;
    
        OIV = input[j];
        if (OIV < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OIV] = % 0.16g [%d]\n",
               i, OIV, j); */
            return 1;
            //OIV = 1e-20;
        }        
        
        total+=OIV * 15.9994;        
        
        j++;
    
        OV = input[j];
        if (OV < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OV] = % 0.16g [%d]\n",
               i, OV, j); */
            return 1;
            //OV = 1e-20;
        }        
        
        total+=OV * 15.9994;        
        
        j++;
    
        OVI = input[j];
        if (OVI < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OVI] = % 0.16g [%d]\n",
               i, OVI, j); */
            return 1;
            //OVI = 1e-20;
        }        
        
        total+=OVI * 15.9994;        
        
        j++;
    
        OVII = input[j];
        if (OVII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OVII] = % 0.16g [%d]\n",
               i, OVII, j); */
            return 1;
            //OVII = 1e-20;
        }        
        
        total+=OVII * 15.9994;        
        
        j++;
    
        OVIII = input[j];
        if (OVIII < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OVIII] = % 0.16g [%d]\n",
               i, OVIII, j); */
            return 1;
            //OVIII = 1e-20;
        }        
        
        total+=OVIII * 15.9994;        
        
        j++;
    
        OIX = input[j];
        if (OIX < 0.0) {
            /* fprintf(stderr, "RNegative[%d][OIX] = % 0.16g [%d]\n",
               i, OIX, j); */
            return 1;
            //OIX = 1e-20;
        }
        
        total+=OIX * 15.9994;
        
        j++;
    
        mdensity = total * mh;

        strip_edot[i] = -HI*HI_c_HI_c[i]*de + HI*HI_ph_HI_ph[i] - HII*HII_c_HII_c[i]*de - HeI*HeI_c_HeI_c[i]*de + HeI*HeI_ph_HeI_ph[i] - HeII*HeII_c_HeII_c[i]*de + HeII*HeII_ph_HeII_ph[i] - HeIII*HeIII_c_HeIII_c[i]*de - OI*OI_c_OI_c[i]*de + OI*OI_ph_OI_ph[i] - OII*OII_c_OII_c[i]*de + OII*OII_ph_OII_ph[i] - OIII*OIII_c_OIII_c[i]*de + OIII*OIII_ph_OIII_ph[i] - OIV*OIV_c_OIV_c[i]*de + OIV*OIV_ph_OIV_ph[i] - OIX*OIX_c_OIX_c[i]*de - OV*OV_c_OV_c[i]*de + OV*OV_ph_OV_ph[i] - OVI*OVI_c_OVI_c[i]*de + OVI*OVI_ph_OVI_ph[i] - OVII*OVII_c_OVII_c[i]*de + OVII*OVII_ph_OVII_ph[i] - OVIII*OVIII_c_OVIII_c[i]*de + OVIII*OVIII_ph_OVIII_ph[i] - compton_comp[i]*de*pow(z + 1.0, 4)*(T - 2.73*z - 2.73);
        
	strip_edot[i] /= mdensity;
    
    }  
    return 0;
}